]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/MEDLoader/SauvMedConvertor.cxx
Salome HOME
Merge from V6_main (04/10/2012)
[tools/medcoupling.git] / src / MEDLoader / SauvMedConvertor.cxx
1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // File      : SauvMedConvertor.cxx
20 // Created   : Tue Aug 16 14:43:20 2011
21 // Author    : Edward AGAPOV (eap)
22 //
23
24 #include "SauvMedConvertor.hxx"
25
26 #include "CellModel.hxx"
27 #include "MEDFileMesh.hxx"
28 #include "MEDFileField.hxx"
29 #include "MEDFileData.hxx"
30 #include "MEDCouplingFieldDouble.hxx"
31
32 #include <iostream>
33 #include <cassert>
34 #include <cmath>
35 #include <queue>
36 #include <limits>
37
38 #include <cstdlib>
39 #include <cstring>
40 #include <fcntl.h>
41
42 #ifdef WIN32
43 #include <io.h>
44 #endif
45
46 #ifndef WIN32
47 #define HAS_XDR
48 #include <unistd.h>
49 #endif
50
51 #ifdef HAS_XDR
52 #include <rpc/xdr.h>
53 #endif
54
55 using namespace SauvUtilities;
56 using namespace ParaMEDMEM;
57 using namespace std;
58
59 namespace
60 {
61   // for ASCII file reader
62   const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
63   const int GIBI_BufferSize   = 16184; // for buffered reading
64
65   using namespace INTERP_KERNEL;
66
67   const size_t MaxMedCellType = NORM_ERROR;
68   const size_t NbGibiCellTypes = 47;
69   const TCellType GibiTypeToMed[NbGibiCellTypes] =
70     {
71       /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2   ,/*3 */ NORM_SEG3   ,/*4 */ NORM_TRI3   ,/*5 */ NORM_ERROR  ,
72       /*6 */ NORM_TRI6   ,/*7 */ NORM_ERROR  ,/*8 */ NORM_QUAD4  ,/*9 */ NORM_ERROR  ,/*10*/ NORM_QUAD8  ,
73       /*11*/ NORM_ERROR  ,/*12*/ NORM_ERROR  ,/*13*/ NORM_ERROR  ,/*14*/ NORM_HEXA8  ,/*15*/ NORM_HEXA20 ,
74       /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR  ,/*19*/ NORM_ERROR  ,/*20*/ NORM_ERROR  ,
75       /*21*/ NORM_ERROR  ,/*22*/ NORM_ERROR  ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5  ,
76       /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR  ,/*28*/ NORM_ERROR  ,/*29*/ NORM_ERROR  ,/*30*/ NORM_ERROR  ,
77       /*31*/ NORM_ERROR  ,/*32*/ NORM_ERROR  ,/*33*/ NORM_ERROR  ,/*34*/ NORM_ERROR  ,/*35*/ NORM_ERROR  ,
78       /*36*/ NORM_ERROR  ,/*37*/ NORM_ERROR  ,/*38*/ NORM_ERROR  ,/*39*/ NORM_ERROR  ,/*40*/ NORM_ERROR  ,
79       /*41*/ NORM_ERROR  ,/*42*/ NORM_ERROR  ,/*43*/ NORM_ERROR  ,/*44*/ NORM_ERROR  ,/*45*/ NORM_ERROR  ,
80       /*46*/ NORM_ERROR  ,/*47*/ NORM_ERROR
81     };
82
83   //================================================================================
84   /*!
85    * \brief Return dimension of a group
86    */
87   //================================================================================
88
89   unsigned getDim( const Group* grp )
90   {
91     return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
92   }
93
94   //================================================================================
95   /*!
96    * \brief Converts connectivity of quadratic elements
97    */
98   //================================================================================
99
100   inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
101                                 const Cell &                            aCell )
102   {
103     if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
104       {
105         Cell* ma = const_cast<Cell*>(&aCell);
106         vector< Node* > new_nodes( ma->_nodes.size() );
107         for ( size_t i = 0; i < new_nodes.size(); ++i )
108           new_nodes[ i ] = ma->_nodes[ conn[ i ]];
109         ma->_nodes.swap( new_nodes );
110       }
111   }
112
113   //================================================================================
114   /*!
115    * \brief Returns a vector of pairs of node indices to inverse a med volume element
116    */
117   //================================================================================
118
119   void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
120                          vector<pair<int,int> > &                swapVec )
121   {
122     swapVec.clear();
123
124     switch ( type )
125       {
126       case NORM_TETRA4:
127         swapVec.resize(1);
128         swapVec[0] = make_pair( 1, 2 );
129         break;
130       case NORM_PYRA5:
131         swapVec.resize(1);
132         swapVec[0] = make_pair( 1, 3 );
133         break;
134       case NORM_PENTA6:
135         swapVec.resize(2);
136         swapVec[0] = make_pair( 1, 2 );
137         swapVec[1] = make_pair( 4, 5 );
138         break;
139       case NORM_HEXA8:
140         swapVec.resize(2);
141         swapVec[0] = make_pair( 1, 3 );
142         swapVec[1] = make_pair( 5, 7 );
143         break;
144       case NORM_TETRA10:
145         swapVec.resize(3);
146         swapVec[0] = make_pair( 1, 2 );
147         swapVec[1] = make_pair( 4, 6 );
148         swapVec[2] = make_pair( 8, 9 );
149         break;
150       case NORM_PYRA13:
151         swapVec.resize(4);
152         swapVec[0] = make_pair( 1, 3 );
153         swapVec[1] = make_pair( 5, 8 );
154         swapVec[2] = make_pair( 6, 7 );
155         swapVec[3] = make_pair( 10, 12 );
156         break;
157       case NORM_PENTA15:
158         swapVec.resize(4);
159         swapVec[0] = make_pair( 1, 2 );
160         swapVec[1] = make_pair( 4, 5 );
161         swapVec[2] = make_pair( 6, 8 );
162         swapVec[3] = make_pair( 9, 11 );
163         break;
164       case NORM_HEXA20:
165         swapVec.resize(7);
166         swapVec[0] = make_pair( 1, 3 );
167         swapVec[1] = make_pair( 5, 7 );
168         swapVec[2] = make_pair( 8, 11 );
169         swapVec[3] = make_pair( 9, 10 );
170         swapVec[4] = make_pair( 12, 15 );
171         swapVec[5] = make_pair( 13, 14 );
172         swapVec[6] = make_pair( 17, 19 );
173         break;
174         //   case NORM_SEG3: no need to reverse edges
175         //     swapVec.resize(1);
176         //     swapVec[0] = make_pair( 1, 2 );
177         //     break;
178       case NORM_TRI6:
179         swapVec.resize(2);
180         swapVec[0] = make_pair( 1, 2 );
181         swapVec[1] = make_pair( 3, 5 );
182         break;
183       case NORM_QUAD8:
184         swapVec.resize(3);
185         swapVec[0] = make_pair( 1, 3 );
186         swapVec[1] = make_pair( 4, 7 );
187         swapVec[2] = make_pair( 5, 6 );
188         break;
189       default:;
190       }
191   }
192
193   //================================================================================
194   /*!
195    * \brief Inverses element orientation using vector of indices to swap
196    */
197   //================================================================================
198
199   inline void reverse(const Cell & aCell, const vector<pair<int,int> > & swapVec )
200   {
201     Cell* ma = const_cast<Cell*>(&aCell);
202     for ( unsigned i = 0; i < swapVec.size(); ++i )
203       std::swap( ma->_nodes[ swapVec[i].first ],
204                  ma->_nodes[ swapVec[i].second ]);
205     if ( swapVec.empty() )
206       ma->_reverse = true;
207     else
208       ma->_reverse = false;
209   }
210   //================================================================================
211   /*!
212    * \brief Comparator of cells by number used for ordering cells within a med group
213    */
214   struct TCellByIDCompare
215   {
216     bool operator () (const Cell* i1, const Cell* i2)
217     {
218       return i1->_number < i2->_number;
219     }
220   };
221   typedef map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
222
223   //================================================================================
224   /*!
225    * \brief Fill Group::_relocTable if necessary
226    */
227   //================================================================================
228
229   void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
230   {
231     if ( !grp->_isProfile ) return;
232
233     // check if relocation table is necessary
234     bool isRelocated = false;
235     unsigned newOrder = 0;
236     TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
237     for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
238       isRelocated = ( c2oIt->second != newOrder );
239
240     if ( isRelocated )
241     {
242       grp->_relocTable.resize( cell2order.size() );
243       for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
244         grp->_relocTable[ c2oIt->second ] = newOrder;
245     }
246   }
247 }
248
249 //================================================================================
250 /*!
251  * \brief Return dimension for the given cell type
252  */
253 //================================================================================
254
255 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
256 {
257   return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
258 }
259
260 //================================================================================
261 /*!
262  * \brief Returns interlace array to transform a quadratic GIBI element to a MED one.
263  *        i-th array item gives node index in GIBI connectivity for i-th MED node
264  */
265 //================================================================================
266
267 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
268 {
269   static vector<const int*> conn;
270   static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
271   static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,7,3};
272   static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
273   static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
274   static const int quad8  [] = {0,2,4,6, 1,3,5,7};
275   static const int tria6  [] = {0,2,4, 1,3,5};
276   static const int seg3   [] = {0,2,1};
277   if ( conn.empty() )
278   {
279     conn.resize( MaxMedCellType + 1, 0 );
280     conn[ NORM_HEXA20 ] = hexa20;
281     conn[ NORM_PENTA15] = penta15;
282     conn[ NORM_PYRA13 ] = pyra13;
283     conn[ NORM_TETRA10] = tetra10;
284     conn[ NORM_SEG3   ] = seg3;
285     conn[ NORM_TRI6   ] = tria6;
286     conn[ NORM_QUAD8  ] = quad8;
287   }
288   return conn[ type ];
289 }
290
291 //================================================================================
292 /*!
293  * \brief Avoid coping sortedNodeIDs
294  */
295 //================================================================================
296
297 Cell::Cell(const Cell& ma)
298   : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
299 {
300   if ( ma._sortedNodeIDs )
301     {
302       _sortedNodeIDs = new int[ _nodes.size() ];
303       std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
304     }
305 }
306
307 //================================================================================
308 /*!
309  * \brief Rerturn the i-th link of face
310  */
311 //================================================================================
312
313 SauvUtilities::Link Cell::link(int i) const
314 {
315   int i2 = ( i + 1 ) % _nodes.size();
316   if ( _reverse )
317     return make_pair( _nodes[i2]->_number, _nodes[i]->_number );
318   else
319     return make_pair( _nodes[i]->_number, _nodes[i2]->_number );
320 }
321
322 //================================================================================
323 /*!
324  * \brief creates if needed and return _sortedNodeIDs
325  */
326 //================================================================================
327
328 const TID* Cell::getSortedNodes() const
329 {
330   if ( !_sortedNodeIDs )
331   {
332     size_t l=_nodes.size();
333     _sortedNodeIDs = new int[ l ];
334
335     for (size_t i=0; i!=l; ++i)
336       _sortedNodeIDs[i]=_nodes[i]->_number;
337     std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
338   }
339   return _sortedNodeIDs;
340 }
341
342 //================================================================================
343 /*!
344  * \brief Compare sorted ids of cell nodes
345  */
346 //================================================================================
347
348 bool Cell::operator< (const Cell& ma) const
349 {
350   if ( _nodes.size() == 1 )
351     return _nodes[0] < ma._nodes[0];
352
353   const int* v1 = getSortedNodes();
354   const int* v2 = ma.getSortedNodes();
355   for ( const int* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
356     if(*v1 != *v2)
357       return *v1 < *v2;
358   return false;
359 }
360
361 //================================================================================
362 /*!
363  * \brief Dump a Cell
364  */
365 //================================================================================
366
367 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
368 {
369   os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0]->_number;
370   for( size_t i=1; i!=ma._nodes.size(); ++i)
371     os << ", " << ma._nodes[i]->_number;
372 #ifdef _DEBUG_
373   os << " > sortedNodes: ";
374   if ( ma._sortedNodeIDs ) {
375     os << "< ";
376     for( size_t i=0; i!=ma._nodes.size(); ++i)
377       os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
378     os << " >";
379   }
380   else {
381     os << "NULL";
382   }
383 #endif
384   return os;
385 }
386
387 //================================================================================
388 /*!
389  * \brief Return nb of elements in the group
390  */
391 //================================================================================
392
393 int Group::size() const
394 {
395   int sizze = 0;
396   if ( !_relocTable.empty() )
397     sizze =  _relocTable.size();
398   else if ( _medGroup )
399     sizze = _medGroup->getNumberOfTuples();
400   else if ( !_cells.empty() )
401     sizze = _cells.size();
402   else
403     for ( size_t i = 0; i < _groups.size(); ++i )
404       sizze += _groups[i]->size();
405   return sizze;
406 }
407
408 //================================================================================
409 /*!
410  * \brief Conver gibi element type to med one
411  */
412 //================================================================================
413
414 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
415 {
416   if ( gibiType < 1 || gibiType > NbGibiCellTypes )
417     return NORM_ERROR;
418
419   return GibiTypeToMed[ gibiType - 1 ];
420 }
421
422 //================================================================================
423 /*!
424  * \brief Conver med element type to gibi one
425  */
426 //================================================================================
427
428 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
429 {
430   for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
431     if ( GibiTypeToMed[ i ] == medGeomType )
432       return i + 1;
433
434   return -1;
435 }
436
437 //================================================================================
438 /*!
439  * \brief Remember the file name
440  */
441 //================================================================================
442
443 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
444 {
445 }
446
447 //================================================================================
448 /*!
449  * \brief Constructor of ASCII sauve file reader
450  */
451 //================================================================================
452
453 ASCIIReader::ASCIIReader(const char* fileName)
454   :FileReader(fileName),
455    _file(-1)
456 {
457 }
458
459 //================================================================================
460 /*!
461  * \brief Return true
462  */
463 //================================================================================
464
465 bool ASCIIReader::isASCII() const
466 {
467   return true;
468 }
469
470 //================================================================================
471 /*!
472  * \brief Try to open an ASCII file
473  */
474 //================================================================================
475
476 bool ASCIIReader::open()
477 {
478 #ifdef WIN32
479   _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
480 #else
481   _file = ::open (_fileName.c_str(), O_RDONLY);
482 #endif
483   if (_file >= 0)
484     {
485       _start  = new char [GIBI_BufferSize]; // working buffer beginning
486       //_tmpBuf = new char [GIBI_MaxOutputLen];
487       _ptr    = _start;
488       _eptr   = _start;
489       _lineNb = 0;
490     }
491   else
492     {
493       //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
494     }
495   return (_file >= 0);
496 }
497
498 //================================================================================
499 /*!
500  * \brief Close the file
501  */
502 //================================================================================
503
504 ASCIIReader::~ASCIIReader()
505 {
506   if (_file >= 0)
507     {
508       ::close (_file);
509       if (_start != 0L)
510         {
511           delete [] _start;
512           //delete [] _tmpBuf;
513           _start = 0;
514         }
515       _file = -1;
516     }
517 }
518
519 //================================================================================
520 /*!
521  * \brief Return a next line of the file
522  */
523 //================================================================================
524
525 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
526 {
527   if ( getLine( line )) return true;
528   if ( raiseOEF )
529     THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
530   return false;
531 }
532
533 //================================================================================
534 /*!
535  * \brief Read a next line of the file if necessary
536  */
537 //================================================================================
538
539 bool ASCIIReader::getLine(char* & line)
540 {
541   bool aResult = true;
542   // Check the state of the buffer;
543   // if there is too little left, read the next portion of data
544   int nBytesRest = _eptr - _ptr;
545   if (nBytesRest < GIBI_MaxOutputLen)
546     {
547       if (nBytesRest > 0)
548         {
549           // move the remaining portion to the buffer beginning
550           for ( int i = 0; i < nBytesRest; ++i )
551             _start[i] = _ptr[i];
552           //memcpy (_tmpBuf, _ptr, nBytesRest);
553           //memcpy (_start, _tmpBuf, nBytesRest);
554         }
555       else
556         {
557           nBytesRest = 0;
558         }
559       _ptr = _start;
560       const int nBytesRead = ::read (_file,
561                                      &_start [nBytesRest],
562                                      GIBI_BufferSize - nBytesRest);
563       nBytesRest += nBytesRead;
564       _eptr = &_start [nBytesRest];
565     }
566   // Check the buffer for the end-of-line
567   char * ptr = _ptr;
568   while (true)
569     {
570       // Check for end-of-the-buffer, the ultimate criterion for termination
571       if (ptr >= _eptr)
572         {
573           if (nBytesRest <= 0)
574             aResult = false;
575           else
576             _eptr[-1] = '\0';
577           break;
578         }
579       // seek the line-feed character
580       if (ptr[0] == '\n')
581         {
582           if (ptr[-1] == '\r')
583             ptr[-1] = '\0';
584           ptr[0] = '\0';
585           ++ptr;
586           break;
587         }
588       ++ptr;
589     }
590   // Output the result
591   line = _ptr;
592   _ptr = ptr;
593   _lineNb++;
594
595   return aResult;
596 }
597
598 //================================================================================
599 /*!
600  * \brief Prepare for iterating over given nb of values
601  *  \param nbToRead - nb of fields to read
602  *  \param nbPosInLine - nb of fields in one line
603  *  \param width - field width
604  *  \param shift - shift from line beginning to the field start
605  */
606 //================================================================================
607
608 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
609 {
610   _nbToRead    = nbToRead;
611   _nbPosInLine = nbPosInLine;
612   _width       = width;
613   _shift       = shift;
614   _iPos = _iRead = 0;
615   if ( _nbToRead )
616     {
617       getNextLine( _curPos );
618       _curPos = _curPos + _shift;
619     }
620   else
621     {
622       _curPos = 0;
623     }
624   _curLocale.clear();
625 }
626
627 //================================================================================
628 /*!
629  * \brief Prepare for iterating over given nb of string values
630  */
631 //================================================================================
632
633 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
634 {
635   init( nbValues, 72 / ( width + 1 ), width, 1 );
636 }
637
638 //================================================================================
639 /*!
640  * \brief Prepare for iterating over given nb of integer values
641  */
642 //================================================================================
643
644 void ASCIIReader::initIntReading(int nbValues)
645 {
646   init( nbValues, 10, 8 );
647 }
648
649 //================================================================================
650 /*!
651  * \brief Prepare for iterating over given nb of real values
652  */
653 //================================================================================
654
655 void ASCIIReader::initDoubleReading(int nbValues)
656 {
657   init( nbValues, 3, 22 );
658
659   // Correction 2 of getDouble(): set "C" numeric locale to read numbers
660   // with dot decimal point separator, as it is in SAUVE files
661   _curLocale = setlocale(LC_NUMERIC, "C");
662 }
663
664 //================================================================================
665 /*!
666  * \brief Return true if not all values have been read
667  */
668 //================================================================================
669
670 bool ASCIIReader::more() const
671 {
672   bool result = false;
673   if ( _iRead < _nbToRead)
674     {
675       if ( _curPos ) result = true;
676     }
677   return result;
678 }
679
680 //================================================================================
681 /*!
682  * \brief Go to the nex value
683  */
684 //================================================================================
685
686 void ASCIIReader::next()
687 {
688   if ( !more() )
689     THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
690   ++_iRead;
691   ++_iPos;
692   if ( _iRead < _nbToRead )
693     {
694       if ( _iPos >= _nbPosInLine )
695         {
696           getNextLine( _curPos );
697           _curPos = _curPos + _shift;
698           _iPos = 0;
699         }
700       else
701         {
702           _curPos = _curPos + _width + _shift;
703         }
704     }
705   else
706     {
707       _curPos = 0;
708       if ( !_curLocale.empty() )
709         {
710           setlocale(LC_NUMERIC, _curLocale.c_str());
711           _curLocale.clear();
712         }
713     }
714 }
715
716 //================================================================================
717 /*!
718  * \brief Return the current integer value
719  */
720 //================================================================================
721
722 int ASCIIReader::getInt() const
723 {
724   // fix for two glued ints (issue 0021009):
725   // Line nb    |   File contents
726   // ------------------------------------------------------------------------------------
727   // 53619905   |       1       2       6       8
728   // 53619906   |                                                                SCALAIRE
729   // 53619907   |    -63312600499       1       0       0       0      -2       0       2
730   //   where -63312600499 is actualy -633 and 12600499
731   char hold=_curPos[_width];
732   _curPos[_width] = '\0';
733   int result = atoi( _curPos );
734   _curPos[_width] = hold;
735   return result;
736   //return atoi(str());
737 }
738
739 //================================================================================
740 /*!
741  * \brief Return the current float value
742  */
743 //================================================================================
744
745 float ASCIIReader::getFloat() const
746 {
747   return getDouble();
748 }
749
750 //================================================================================
751 /*!
752  * \brief Return the current double value
753  */
754 //================================================================================
755
756 double ASCIIReader::getDouble() const
757 {
758   //std::string aStr (_curPos);
759
760   // Correction: add missing 'E' specifier
761   // int aPosStart = aStr.find_first_not_of(" \t");
762   // if (aPosStart < (int)aStr.length()) {
763   //   int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
764   //   if (aPosSign < (int)aStr.length()) {
765   //     if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
766   //       aStr.insert(aPosSign, "E", 1);
767   //   }
768   // }
769
770   // Different Correction (more optimal)
771   // Sample:
772   //  0.00000000000000E+00 -2.37822406690632E+01  6.03062748797469E+01
773   //  7.70000000000000-100  7.70000000000000+100  7.70000000000000+100
774   //0123456789012345678901234567890123456789012345678901234567890123456789
775   const size_t posE = 18;
776   std::string aStr (_curPos);
777   if ( aStr.find('E') < 0 && aStr.find('e') < 0 )
778     {
779       if ( aStr.size() < posE+1 )
780         THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
781       aStr.insert( posE, "E", 1 );
782       return atof(aStr.c_str());
783     }
784   return atof( _curPos );
785 }
786
787 //================================================================================
788 /*!
789  * \brief Return the current string value
790  */
791 //================================================================================
792
793 string ASCIIReader::getName() const
794 {
795   int len = _width;
796   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
797     len--;
798   return string( _curPos, len );
799 }
800
801 //================================================================================
802 /*!
803  * \brief Constructor of a binary sauve file reader
804  */
805 //================================================================================
806
807 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
808 {
809 }
810
811 //================================================================================
812 /*!
813  * \brief Close the XDR sauve file
814  */
815 //================================================================================
816
817 XDRReader::~XDRReader()
818 {
819 #ifdef HAS_XDR  
820   if ( _xdrs_file )
821     {
822       xdr_destroy((XDR*)_xdrs);
823       free((XDR*)_xdrs);
824       ::fclose(_xdrs_file);
825       _xdrs_file = NULL;
826     }
827 #endif
828 }
829
830 //================================================================================
831 /*!
832  * \brief Return false
833  */
834 //================================================================================
835
836 bool XDRReader::isASCII() const
837 {
838   return false;
839 }
840
841 //================================================================================
842 /*!
843  * \brief Try to open an XRD file
844  */
845 //================================================================================
846
847 bool XDRReader::open()
848 {
849   bool xdr_ok = false;
850 #ifdef HAS_XDR
851   if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
852     {
853       _xdrs = (XDR *)malloc(sizeof(XDR));
854       xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
855
856       const int maxsize = 10;
857       char icha[maxsize+1];
858       char* icha2 = icha;
859       if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
860         {
861           icha[maxsize] = '\0';
862           xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
863         }
864       if ( !xdr_ok )
865         {
866           xdr_destroy((XDR*)_xdrs);
867           free((XDR*)_xdrs);
868           fclose(_xdrs_file);
869           _xdrs_file = NULL;
870         }
871     }
872 #endif
873   return xdr_ok;
874 }
875
876 //================================================================================
877 /*!
878  * \brief A stub
879  */
880 //================================================================================
881
882 bool XDRReader::getNextLine (char* &, bool )
883 {
884   return true;
885 }
886
887 //================================================================================
888 /*!
889  * \brief Prepare for iterating over given nb of values
890  *  \param nbToRead - nb of fields to read
891  *  \param width - field width
892  */
893 //================================================================================
894
895 void XDRReader::init( int nbToRead, int width/*=0*/ )
896 {
897   if(_iRead < _nbToRead)
898     {
899       cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << endl;
900       cout << "Unfinished iteration before new one !" << endl;
901       THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
902     }
903   _iRead    = 0;
904   _nbToRead = nbToRead;
905   _width    = width;
906 }
907
908 //================================================================================
909 /*!
910  * \brief Prepare for iterating over given nb of string values
911  */
912 //================================================================================
913
914 void XDRReader::initNameReading(int nbValues, int width)
915 {
916   init( nbValues, width );
917   _xdr_kind = _xdr_kind_char;
918   if(nbValues*width)
919     {
920       unsigned int nels = nbValues*width;
921       _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
922 #ifdef HAS_XDR
923       xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
924 #endif
925       _xdr_cvals[nels] = '\0';
926     }
927 }
928
929 //================================================================================
930 /*!
931  * \brief Prepare for iterating over given nb of integer values
932  */
933 //================================================================================
934
935 void XDRReader::initIntReading(int nbValues)
936 {
937   init( nbValues );
938   _xdr_kind = _xdr_kind_int;
939   if(nbValues)
940     {
941 #ifdef HAS_XDR
942       unsigned int nels = nbValues;
943       unsigned int actual_nels;
944       _xdr_ivals = (int*)malloc(nels*sizeof(int));
945       xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
946 #endif
947     }
948 }
949
950 //================================================================================
951 /*!
952  * \brief Prepare for iterating over given nb of real values
953  */
954 //================================================================================
955
956 void XDRReader::initDoubleReading(int nbValues)
957 {
958   init( nbValues );
959   _xdr_kind = _xdr_kind_double;
960   if(nbValues)
961     {
962 #ifdef HAS_XDR
963       unsigned int nels = nbValues;
964       unsigned int actual_nels;
965       _xdr_dvals = (double*)malloc(nels*sizeof(double));
966       xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
967 #endif
968     }
969 }
970
971 //================================================================================
972 /*!
973  * \brief Return true if not all values have been read
974  */
975 //================================================================================
976
977 bool XDRReader::more() const
978 {
979   return _iRead < _nbToRead;
980 }
981
982 //================================================================================
983 /*!
984  * \brief Go to the nex value
985  */
986 //================================================================================
987
988 void XDRReader::next()
989 {
990   if ( !more() )
991     THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
992
993   ++_iRead;
994   if ( _iRead < _nbToRead )
995     {
996     }
997   else
998     {
999       if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
1000       if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
1001       if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
1002       _xdr_kind = _xdr_kind_null;
1003     }
1004 }
1005
1006 //================================================================================
1007 /*!
1008  * \brief Return the current integer value
1009  */
1010 //================================================================================
1011
1012 int XDRReader::getInt() const
1013 {
1014   if(_iRead < _nbToRead)
1015     {
1016       return _xdr_ivals[_iRead];
1017     }
1018   else
1019     {
1020       int result = 0;
1021 #ifdef HAS_XDR
1022       xdr_int((XDR*)_xdrs, &result);
1023 #endif
1024       return result;
1025     }
1026 }
1027
1028 //================================================================================
1029 /*!
1030  * \brief Return the current float value
1031  */
1032 //================================================================================
1033
1034 float  XDRReader::getFloat() const
1035 {
1036   float result = 0;
1037 #ifdef HAS_XDR
1038   xdr_float((XDR*)_xdrs, &result);
1039 #endif
1040   return result;
1041 }
1042
1043 //================================================================================
1044 /*!
1045  * \brief Return the current double value
1046  */
1047 //================================================================================
1048
1049 double XDRReader::getDouble() const
1050 {
1051   if(_iRead < _nbToRead)
1052     {
1053       return _xdr_dvals[_iRead];
1054     }
1055   else
1056     {
1057       double result = 0;
1058 #ifdef HAS_XDR
1059       xdr_double((XDR*)_xdrs, &result);
1060 #endif
1061       return result;
1062     }
1063 }
1064
1065 //================================================================================
1066 /*!
1067  * \brief Return the current string value
1068  */
1069 //================================================================================
1070
1071 std::string XDRReader::getName() const
1072 {
1073   int len = _width;
1074   char* s = _xdr_cvals + _iRead*_width;
1075   while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
1076     len--;
1077   return string( s, len );
1078 }
1079
1080 //================================================================================
1081 /*!
1082  * \brief Throw an exception if not all needed data is present
1083  */
1084 //================================================================================
1085
1086 void IntermediateMED::checkDataAvailability() const throw(INTERP_KERNEL::Exception)
1087 {
1088   if ( _spaceDim == 0 )
1089     THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
1090
1091   if ( _groups.empty() )
1092     THROW_IK_EXCEPTION("No elements have been read");
1093
1094   if ( _points.empty() || _nbNodes == 0 )
1095     THROW_IK_EXCEPTION("Nodes of elements are not filled");
1096
1097   if ( _coords.empty() )
1098     THROW_IK_EXCEPTION("Node coordinates are missing");
1099
1100   if ( _coords.size() < _nbNodes * _spaceDim )
1101     THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
1102 }
1103
1104 //================================================================================
1105 /*!
1106  * \brief Makes ParaMEDMEM::MEDFileData from self
1107  */
1108 //================================================================================
1109
1110 ParaMEDMEM::MEDFileData* IntermediateMED::convertInMEDFileDS()
1111 {
1112   MEDCouplingAutoRefCountObjectPtr< MEDFileUMesh >  mesh   = makeMEDFileMesh();
1113   MEDCouplingAutoRefCountObjectPtr< MEDFileFields > fields = makeMEDFileFields(mesh);
1114
1115   MEDCouplingAutoRefCountObjectPtr< MEDFileMeshes > meshes = MEDFileMeshes::New();
1116   MEDCouplingAutoRefCountObjectPtr< MEDFileData >  medData = MEDFileData::New();
1117   meshes->pushMesh( mesh );
1118   medData->setMeshes( meshes );
1119   if ( fields ) medData->setFields( fields );
1120
1121   medData->incrRef();
1122   return medData;
1123 }
1124
1125 //================================================================================
1126 /*!
1127  * \brief Creates ParaMEDMEM::MEDFileUMesh from its data
1128  */
1129 //================================================================================
1130
1131 ParaMEDMEM::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
1132 {
1133   // check if all needed piles are present
1134   checkDataAvailability();
1135
1136   // set long names
1137   setGroupLongNames();
1138
1139   // fix element orientation
1140   if ( _spaceDim == 2 || _spaceDim == 1 )
1141     orientElements2D();
1142   else if ( _spaceDim == 3 )
1143     orientElements3D();
1144
1145   // process groups
1146   decreaseHierarchicalDepthOfSubgroups();
1147   eraseUselessGroups();
1148   //detectMixDimGroups();
1149
1150   // assign IDs
1151   _points.numberNodes();
1152   numberElements();
1153
1154   // make the med mesh
1155
1156   MEDFileUMesh* mesh = MEDFileUMesh::New();
1157
1158   DataArrayDouble *coords = getCoords();
1159   setConnectivity( mesh, coords );
1160   setGroups( mesh );
1161
1162   coords->decrRef();
1163
1164   if ( !mesh->getName() || strlen( mesh->getName() ) == 0 )
1165     mesh->setName( "MESH" );
1166
1167   return mesh;
1168 }
1169
1170 //================================================================================
1171 /*!
1172  * \brief Set long names to groups
1173  */
1174 //================================================================================
1175
1176 void IntermediateMED::setGroupLongNames()
1177 {
1178   // IMP 0020434: mapping GIBI names to MED names
1179   // set med names to objects (mesh, fields, support, group or other)
1180
1181   set<int> treatedGroups;
1182
1183   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
1184   for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
1185     {
1186       if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
1187
1188       SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
1189
1190       // if there are several names for grp then the 1st name is the name
1191       // of grp and the rest ones are names of groups referring grp (issue 0021311)
1192       const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
1193       if ( !isRefName )
1194         {
1195           grp._name = _mapStrings[ itGIBItoMED->med_id ];
1196         }
1197       else if ( !grp._refNames.empty() && grp._refNames.back().empty() )
1198         {
1199           for ( unsigned i = 0; i < grp._refNames.size(); ++i )
1200             if ( grp._refNames[i].empty() )
1201               grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
1202         }
1203       else
1204         {
1205           grp._refNames.push_back( _mapStrings[ (*itGIBItoMED).med_id ]);
1206         }
1207     }
1208 }
1209
1210 //================================================================================
1211 /*!
1212  * \brief Set long names to fields
1213  */
1214 //================================================================================
1215
1216 void IntermediateMED::setFieldLongNames(set< string >& usedNames)
1217 {
1218   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
1219   for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
1220     {
1221       if (itGIBItoMED->gibi_pile == PILE_FIELD)
1222         {
1223           _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
1224         }
1225       else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
1226         {
1227           _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
1228         }
1229     } // iterate on _listGIBItoMED_cham
1230
1231   for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
1232     {
1233       string medName  = _mapStrings[itGIBItoMED->med_id];
1234       string gibiName = _mapStrings[itGIBItoMED->gibi_id];
1235
1236       bool name_found = false;
1237       for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
1238         {
1239           vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
1240           for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
1241             {
1242               if (medName.find( fields[ifi]->_name + "." ) == 0 )
1243                 {
1244                   vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
1245                   int nbSub = aSubDs.size();
1246                   for (int isu = 0; isu < nbSub; isu++)
1247                     for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
1248                       {
1249                         if (aSubDs[isu].compName(ico) == gibiName)
1250                           {
1251                             string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
1252                             fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
1253                           }
1254                       }
1255                 }
1256             }
1257         }
1258     } // iterate on _listGIBItoMED_comp
1259
1260   for ( size_t i = 0; i < _nodeFields.size() ; i++)
1261     usedNames.insert( _nodeFields[i]->_name );
1262   for ( size_t i = 0; i < _cellFields.size() ; i++)
1263     usedNames.insert( _cellFields[i]->_name );
1264 }
1265
1266 //================================================================================
1267 /*!
1268  * \brief Decrease hierarchical depth of subgroups
1269  */
1270 //================================================================================
1271
1272 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
1273 {
1274   for (size_t i=0; i!=_groups.size(); ++i)
1275   {
1276     Group& grp = _groups[i];
1277     for (size_t j = 0; j < grp._groups.size(); ++j )
1278     {
1279       Group & sub_grp = *grp._groups[j];
1280       if ( !sub_grp._groups.empty() )
1281       {
1282         // replace j with its 1st subgroup
1283         grp._groups[j] = sub_grp._groups[0];
1284         // push back the rest subs
1285         grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
1286       }
1287     }
1288     // remove empty sub-_groups
1289     vector< Group* > newSubGroups;
1290     newSubGroups.reserve( grp._groups.size() );
1291     for (size_t j = 0; j < grp._groups.size(); ++j )
1292       if ( !grp._groups[j]->empty() )
1293         newSubGroups.push_back( grp._groups[j] );
1294     if ( newSubGroups.size() < grp._groups.size() )
1295       grp._groups.swap( newSubGroups );
1296   }
1297 }
1298
1299 //================================================================================
1300 /*!
1301  * \brief Erase _groups that won't be converted
1302  */
1303 //================================================================================
1304
1305 void IntermediateMED::eraseUselessGroups()
1306 {
1307   // propagate _isProfile=true to sub-groups of composite groups
1308   // for (size_t int i=0; i!=_groups.size(); ++i)
1309   // {
1310   //   Group* grp = _groups[i];
1311   //   if ( grp->_isProfile && !grp->_groups.empty() )
1312   //     for (size_t j = 0; j < grp->_groups.size(); ++j )
1313   //       grp->_groups[j]->_isProfile=true;
1314   // }
1315   std::set<Group*> groups2convert;
1316   // keep not named sub-groups of field supports
1317   for (size_t i=0; i!=_groups.size(); ++i)
1318   {
1319     Group& grp = _groups[i];
1320     if ( grp._isProfile && !grp._groups.empty() )
1321       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
1322   }
1323
1324   // keep named groups and their subgroups
1325   for (size_t i=0; i!=_groups.size(); ++i)
1326   {
1327     Group& grp = _groups[i];
1328     if ( !grp._name.empty() && !grp.empty() )
1329     {
1330       groups2convert.insert( &grp );
1331       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
1332     }
1333   }
1334   // erase groups that are not in groups2convert and not _isProfile
1335   for (size_t i=0; i!=_groups.size(); ++i)
1336   {
1337     Group* grp = &_groups[i];
1338     if ( !grp->_isProfile && !groups2convert.count( grp ) )
1339     {
1340       grp->_cells.clear();
1341       grp->_groups.clear();
1342     }
1343   }
1344 }
1345
1346 //================================================================================
1347 /*!
1348  * \brief Detect _groups of mixed dimension
1349  */
1350 //================================================================================
1351
1352 void IntermediateMED::detectMixDimGroups()
1353 {
1354   //hasMixedCells = false;
1355   for ( size_t i=0; i < _groups.size(); ++i )
1356   {
1357     Group& grp = _groups[i];
1358     if ( grp._groups.size() < 2 )
1359       continue;
1360
1361     // check if sub-groups have different dimension
1362     unsigned dim1 = getDim( &grp );
1363     for ( size_t j = 1; j  < grp._groups.size(); ++j )
1364     {
1365       unsigned dim2 = getDim( grp._groups[j] );
1366       if ( dim1 != dim2 )
1367       {
1368         grp._cells.clear();
1369         grp._groups.clear();
1370         if ( !grp._name.empty() )
1371           cout << "Erase a group with elements of different dim |" << grp._name << "|"<< endl;
1372         break;
1373       }
1374     }
1375   }
1376 }
1377
1378 //================================================================================
1379 /*!
1380  * \brief Fix connectivity of elements in 2D space
1381  */
1382 //================================================================================
1383
1384 void IntermediateMED::orientElements2D()
1385 {
1386   set<Cell>::const_iterator elemIt, elemEnd;
1387   vector< pair<int,int> > swapVec;
1388
1389   // ------------------------------------
1390   // fix connectivity of quadratic edges
1391   // ------------------------------------
1392   set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
1393   if ( !quadEdges.empty() )
1394     {
1395       elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
1396       for ( ; elemIt != elemEnd; ++elemIt )
1397         ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
1398     }
1399
1400   CellsByDimIterator faceIt( *this, 2 );
1401   while ( const set<Cell > * faces = faceIt.nextType() )
1402     {
1403       TCellType cellType = faceIt.type();
1404       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
1405
1406       getReverseVector( cellType, swapVec );
1407
1408       // ------------------------------------
1409       // fix connectivity of quadratic faces
1410       // ------------------------------------
1411       if ( isQuadratic )
1412         for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
1413           ConvertQuadratic( cellType, *elemIt );
1414
1415       // --------------------------
1416       // orient faces clockwise
1417       // --------------------------
1418       int iQuad = isQuadratic ? 2 : 1;
1419       for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
1420         {
1421           // look for index of the most left node
1422           int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
1423           double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
1424           for ( iNode = 1; iNode < nbNodes; ++iNode )
1425             if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
1426               minX = x, iLeft = iNode;
1427
1428           // indeces of the nodes neighboring the most left one
1429           int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
1430           int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
1431           // find components of prev-left and left-next vectors
1432           double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
1433           double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
1434           double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
1435           double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
1436           double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
1437           double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
1438           double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
1439           double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
1440           // normalise y of the vectors
1441           double modPL = sqrt ( xPL * xPL + yPL * yPL );
1442           double modLN = sqrt ( xLN * xLN + yLN * yLN );
1443           if ( modLN > std::numeric_limits<double>::min() &&
1444                modPL > std::numeric_limits<double>::min() )
1445             {
1446               yPL /= modPL;
1447               yLN /= modLN;
1448               // summary direction of neighboring links must be positive
1449               bool clockwise = ( yPL + yLN > 0 );
1450               if ( !clockwise )
1451                 reverse( *elemIt, swapVec );
1452             }
1453         }
1454     }
1455 }
1456
1457 //================================================================================
1458 /*!
1459  * \brief Fix connectivity of elements in 3D space
1460  */
1461 //================================================================================
1462
1463 void IntermediateMED::orientElements3D()
1464 {
1465   // set _reverse flags of faces
1466   orientFaces3D();
1467
1468   // -----------------
1469   // fix connectivity
1470   // -----------------
1471
1472   set<Cell>::const_iterator elemIt, elemEnd;
1473   vector< pair<int,int> > swapVec;
1474
1475   for ( int dim = 1; dim <= 3; ++dim )
1476   {
1477     CellsByDimIterator cellsIt( *this, dim );
1478     while ( const set<Cell > * elems = cellsIt.nextType() )
1479     {
1480       TCellType cellType = cellsIt.type();
1481       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
1482       getReverseVector( cellType, swapVec );
1483
1484       elemIt = elems->begin(), elemEnd = elems->end();
1485       for ( ; elemIt != elemEnd; elemIt++ )
1486       {
1487         // GIBI connectivity -> MED one
1488         if( isQuadratic )
1489           ConvertQuadratic( cellType, *elemIt );
1490
1491         // reverse faces
1492         if ( elemIt->_reverse )
1493           reverse ( *elemIt, swapVec );
1494       }
1495     }
1496   }
1497
1498   orientVolumes();
1499 }
1500
1501 //================================================================================
1502 /*!
1503  * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
1504  */
1505 //================================================================================
1506
1507 void IntermediateMED::orientFaces3D()
1508 {
1509   // fill map of links and their faces
1510   set<const Cell*> faces;
1511   map<const Cell*, Group*> fgm;
1512   map<Link, list<const Cell*> > linkFacesMap;
1513   map<Link, list<const Cell*> >::iterator lfIt, lfIt2;
1514
1515   for (size_t i=0; i!=_groups.size(); ++i)
1516     {
1517       Group& grp = _groups[i];
1518       if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
1519         for ( size_t j = 0; j < grp._cells.size(); ++j )
1520           if ( faces.insert( grp._cells[j] ).second )
1521             {
1522               for ( size_t k = 0; k < grp._cells[j]->_nodes.size(); ++k )
1523                 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
1524               fgm.insert( make_pair( grp._cells[j], &grp ));
1525             }
1526     }
1527   // dump linkFacesMap
1528   //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
1529   //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
1530   //       list<const Cell*> & fList = lfIt->second;
1531   //       list<const Cell*>::iterator fIt = fList.begin();
1532   //       for ( ; fIt != fList.end(); fIt++ )
1533   //         cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
1534   //     }
1535
1536   // Each oriented link must appear in one face only, else a face is reversed.
1537
1538   queue<const Cell*> faceQueue; /* the queue contains well oriented faces
1539                                      whose neighbors orientation is to be checked */
1540   bool manifold = true;
1541   while ( !linkFacesMap.empty() )
1542     {
1543       if ( faceQueue.empty() )
1544         {
1545           assert( !linkFacesMap.begin()->second.empty() );
1546           faceQueue.push( linkFacesMap.begin()->second.front() );
1547         }
1548       while ( !faceQueue.empty() )
1549         {
1550           const Cell* face = faceQueue.front();
1551           faceQueue.pop();
1552
1553           // loop on links of <face>
1554           for ( int i = 0; i < (int)face->_nodes.size(); ++i )
1555             {
1556               Link link = face->link( i );
1557               // find the neighbor faces
1558               lfIt = linkFacesMap.find( link );
1559               int nbFaceByLink = 0;
1560               list< const Cell* > ml;
1561               if ( lfIt != linkFacesMap.end() )
1562                 {
1563                   list<const Cell*> & fList = lfIt->second;
1564                   list<const Cell*>::iterator fIt = fList.begin();
1565                   assert( fIt != fList.end() );
1566                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
1567                     {
1568                       ml.push_back( *fIt );
1569                       if ( *fIt != face ) // wrongly oriented neighbor face
1570                         {
1571                           const Cell* badFace = *fIt;
1572                           // reverse and remove badFace from linkFacesMap
1573                           for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
1574                             {
1575                               Link badlink = badFace->link( j );
1576                               if ( badlink == link ) continue;
1577                               lfIt2 = linkFacesMap.find( badlink );
1578                               if ( lfIt2 != linkFacesMap.end() )
1579                                 {
1580                                   list<const Cell*> & ff = lfIt2->second;
1581                                   ff.erase( find( ff.begin(), ff.end(), badFace ));
1582                                   if ( ff.empty() )
1583                                     linkFacesMap.erase( lfIt2 );
1584                                 }
1585                             }
1586                           badFace->_reverse = true; // reverse
1587                           //INFOS_MED( "REVERSE " << *badFace );
1588                           faceQueue.push( badFace );
1589                         }
1590                     }
1591                   linkFacesMap.erase( lfIt );
1592                 }
1593               // add good neighbors to the queue
1594               Link revLink( link.second, link.first );
1595               lfIt = linkFacesMap.find( revLink );
1596               if ( lfIt != linkFacesMap.end() )
1597                 {
1598                   list<const Cell*> & fList = lfIt->second;
1599                   list<const Cell*>::iterator fIt = fList.begin();
1600                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
1601                     {
1602                       ml.push_back( *fIt );
1603                       if ( *fIt != face )
1604                         faceQueue.push( *fIt );
1605                     }
1606                   linkFacesMap.erase( lfIt );
1607                 }
1608               if ( nbFaceByLink > 2 )
1609                 {
1610                   if ( manifold )
1611                     {
1612                       list<const Cell*>::iterator ii = ml.begin();
1613                       cout << nbFaceByLink << " faces by 1 link:" << endl;
1614                       for( ; ii!= ml.end(); ii++ )
1615                         cout << "in sub-mesh <" << fgm[ *ii ]->_name << "> " << **ii << endl;
1616                     }
1617                   manifold = false;
1618                 }
1619             } // loop on links of the being checked face
1620         } // loop on the face queue
1621     } // while ( !linkFacesMap.empty() )
1622
1623   if ( !manifold )
1624     cout << " -> Non manifold mesh, faces orientation may be incorrect" << endl;
1625 }
1626
1627 //================================================================================
1628 /*!
1629  * \brief Orient volumes according to MED conventions:
1630  * normal of a bottom (first) face should be outside
1631  */
1632 //================================================================================
1633
1634 void IntermediateMED::orientVolumes()
1635 {
1636   set<Cell>::const_iterator elemIt, elemEnd;
1637   vector< pair<int,int> > swapVec;
1638
1639   CellsByDimIterator cellsIt( *this, 3 );
1640   while ( const set<Cell > * elems = cellsIt.nextType() )
1641     {
1642       TCellType cellType = cellsIt.type();
1643       elemIt = elems->begin(), elemEnd = elems->end();
1644       int nbBottomNodes = 0;
1645       switch ( cellType )
1646         {
1647         case NORM_TETRA4:
1648         case NORM_TETRA10:
1649         case NORM_PENTA6:
1650         case NORM_PENTA15:
1651           nbBottomNodes = 3; break;
1652         case NORM_PYRA5:
1653         case NORM_PYRA13:
1654         case NORM_HEXA8:
1655         case NORM_HEXA20:
1656           nbBottomNodes = 4; break;
1657         default: continue;
1658         }
1659       getReverseVector( cellType, swapVec );
1660
1661       for ( ; elemIt != elemEnd; elemIt++ )
1662         {
1663           // find a normal to the bottom face
1664           const double* n[4];
1665           n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
1666           n[1] = nodeCoords( elemIt->_nodes[1]);
1667           n[2] = nodeCoords( elemIt->_nodes[2]);
1668           n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
1669           double vec01[3]; // vector n[0]-n[1]
1670           vec01[0] = n[1][0] - n[0][0];
1671           vec01[1] = n[1][1] - n[0][1];
1672           vec01[2] = n[1][2] - n[0][2];
1673           double vec02 [3]; // vector n[0]-n[2]
1674           vec02[0] = n[2][0] - n[0][0];
1675           vec02[1] = n[2][1] - n[0][1];
1676           vec02[2] = n[2][2] - n[0][2];
1677           double normal [3]; // vec01 ^ vec02
1678           normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
1679           normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
1680           normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
1681           // check if the 102 angle is convex
1682           if ( nbBottomNodes > 3 )
1683             {
1684               const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
1685               double vec03 [3];  // vector n[0]-n3
1686               vec03[0] = n3[0] - n[0][0];
1687               vec03[1] = n3[1] - n[0][1];
1688               vec03[2] = n3[2] - n[0][2];
1689               if ( fabs( normal[0]+normal[1]+normal[2] ) <= numeric_limits<double>::max() ) // vec01 || vec02
1690                 {
1691                   normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
1692                   normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
1693                   normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
1694                 }
1695               else
1696                 {
1697                   double vec [3]; // normal ^ vec01
1698                   vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
1699                   vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
1700                   vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
1701                   double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
1702                   if ( dot2 < 0 ) // concave -> reverse normal
1703                     {
1704                       normal[0] *= -1;
1705                       normal[1] *= -1;
1706                       normal[2] *= -1;
1707                     }
1708                 }
1709             }
1710           // direction from top to bottom
1711           double tbDir[3];
1712           tbDir[0] = n[0][0] - n[3][0];
1713           tbDir[1] = n[0][1] - n[3][1];
1714           tbDir[2] = n[0][2] - n[3][2];
1715
1716           // compare 2 directions: normal and top-bottom
1717           double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
1718           if ( dot < 0. ) // need reverse
1719             reverse( *elemIt, swapVec );
1720
1721         } // loop on volumes of one geometry
1722     } // loop on 3D geometry types
1723
1724 }
1725
1726 //================================================================================
1727 /*!
1728  * \brief Assign new IDs to nodes by skipping not used nodes and return their number
1729  */
1730 //================================================================================
1731
1732 int NodeContainer::numberNodes()
1733 {
1734   int id = 1;
1735   for ( size_t i = 0; i < _nodes.size(); ++i )
1736     for ( size_t j = 0; j < _nodes[i].size(); ++j )
1737       if ( _nodes[i][j].isUsed() )
1738         _nodes[i][j]._number = id++;
1739   return id-1;
1740 }
1741
1742
1743 //================================================================================
1744 /*!
1745  * \brief Assign new IDs to elements
1746  */
1747 //================================================================================
1748
1749 void IntermediateMED::numberElements()
1750 {
1751   set<Cell>::const_iterator elemIt, elemEnd;
1752
1753   // numbering _cells of type NORM_POINT1 by node number
1754   {
1755     const set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
1756     elemIt = points.begin(), elemEnd = points.end();
1757     for ( ; elemIt != elemEnd; ++elemIt )
1758       elemIt->_number = elemIt->_nodes[0]->_number;
1759   }
1760
1761   // numbering 1D-3D _cells
1762   for ( int dim = 1; dim <= 3; ++dim )
1763     {
1764       // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
1765       bool ok = true, renumEntity = false;
1766       CellsByDimIterator cellsIt( *this, dim );
1767       int prevNbElems = 0;
1768       while ( const set<Cell> * typeCells = cellsIt.nextType() )
1769         {
1770           TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
1771           for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1772             {
1773               if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
1774               if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
1775             }
1776           TID typeSize = typeCells->size();
1777           if ( typeSize != maxNumber - minNumber + 1 )
1778             ok = false;
1779           if ( prevNbElems != 0 ) {
1780             if ( minNumber == 1 )
1781               renumEntity = true;
1782             else if ( prevNbElems+1 != (int)minNumber )
1783               ok = false;
1784           }
1785           prevNbElems += typeSize;
1786         }
1787
1788       if ( ok && renumEntity ) // each geom type was numerated separately
1789         {
1790           cellsIt.init( dim );
1791           prevNbElems = cellsIt.nextType()->size(); // no need to renumber the first type
1792           while ( const set<Cell> * typeCells = cellsIt.nextType() )
1793             {
1794               for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1795                 elemIt->_number += prevNbElems;
1796               prevNbElems += typeCells->size();
1797             }
1798         }
1799       if ( !ok )
1800         {
1801           int cellID=1;
1802           cellsIt.init( dim );
1803           while ( const set<Cell> * typeCells = cellsIt.nextType() )
1804             for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
1805               elemIt->_number = cellID++;
1806         }
1807     }
1808 }
1809
1810 //================================================================================
1811 /*!
1812  * \brief Creates coord array
1813  */
1814 //================================================================================
1815
1816 ParaMEDMEM::DataArrayDouble * IntermediateMED::getCoords()
1817 {
1818   DataArrayDouble* coordArray = DataArrayDouble::New();
1819   coordArray->alloc( _nbNodes, _spaceDim );
1820   double * coordPrt = coordArray->getPointer();
1821   for ( int i = 0, nb = _points.size(); i < nb; ++i )
1822     {
1823       Node* n = getNode( i+1 );
1824       if ( n->isUsed() )
1825         {
1826           const double* nCoords = nodeCoords( n );
1827           std::copy( nCoords, nCoords+_spaceDim, coordPrt );
1828           coordPrt += _spaceDim;
1829         }
1830     }
1831   return coordArray;
1832 }
1833
1834 //================================================================================
1835 /*!
1836  * \brief Sets connectivity of elements to the mesh
1837  *  \param mesh - mesh to fill in
1838  *  \param coords - coordinates that must be shared by all meshes of different dim
1839  */
1840 //================================================================================
1841
1842 void IntermediateMED::setConnectivity( ParaMEDMEM::MEDFileUMesh*    mesh,
1843                                        ParaMEDMEM::DataArrayDouble* coords )
1844 {
1845   int meshDim = 0;
1846
1847   mesh->setCoords( coords );
1848
1849   set<Cell>::const_iterator elemIt, elemEnd;
1850   for ( int dim = 3; dim > 0; --dim )
1851   {
1852     CellsByDimIterator dimCells( *this, dim );
1853
1854     int nbOfCells = 0;
1855     while ( const std::set<Cell > * cells = dimCells.nextType() )
1856       nbOfCells += cells->size();
1857     if ( nbOfCells == 0 )
1858       continue;
1859
1860     if ( !meshDim ) meshDim = dim;
1861
1862     MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
1863     dimMesh->setCoords( coords );
1864     dimMesh->setMeshDimension( dim );
1865     dimMesh->allocateCells( nbOfCells );
1866
1867     int prevNbCells = 0;
1868     dimCells.init( dim );
1869     while ( const std::set<Cell > * cells = dimCells.nextType() )
1870       {
1871         // fill connectivity array to take into account order of elements in the sauv file
1872         const int nbCellNodes = cells->begin()->_nodes.size();
1873         vector< TID > connectivity( cells->size() * nbCellNodes );
1874         int * nodalConnOfCell;
1875         for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
1876           {
1877             const Cell& cell = *elemIt;
1878             const int index = cell._number - 1 - prevNbCells;
1879             nodalConnOfCell = &connectivity[ index * nbCellNodes ];
1880             if ( cell._reverse )
1881               for ( int i = nbCellNodes-1; i >= 0; --i )
1882                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
1883             else
1884               for ( int i = 0; i < nbCellNodes; ++i )
1885                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
1886           }
1887         prevNbCells += cells->size();
1888
1889         // fill dimMesh
1890         TCellType cellType = dimCells.type();
1891         nodalConnOfCell = &connectivity[0];
1892         for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
1893           dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
1894       }
1895     dimMesh->finishInsertingCells();
1896     mesh->setMeshAtLevel( dim - meshDim, dimMesh );
1897     dimMesh->decrRef();
1898   }
1899 }
1900
1901 //================================================================================
1902 /*!
1903  * \brief Fill in the mesh with groups
1904  *  \param mesh - mesh to fill in
1905  */
1906 //================================================================================
1907
1908 void IntermediateMED::setGroups( ParaMEDMEM::MEDFileUMesh* mesh )
1909 {
1910   bool isMeshNameSet = false;
1911   const int meshDim = mesh->getMeshDimension();
1912   for ( int dim = 0; dim <= meshDim; ++dim )
1913     {
1914       const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
1915
1916       vector<const DataArrayInt *> medGroups;
1917       vector<MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > refGroups;
1918       for ( size_t i = 0; i < _groups.size(); ++i )
1919         {
1920           Group& grp = _groups[i];
1921           if ( (int)getDim( &grp ) != dim &&
1922                grp._groups.empty() ) // to allow groups on diff dims
1923             continue;
1924           // convert only named groups or field supports
1925           if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
1926             continue;
1927           //if ( grp._medGroup ) continue; // already converted
1928
1929           // sort cells by ID and remember their initial order in the group
1930           TCellToOrderMap cell2order;
1931           unsigned orderInGroup = 0;
1932           vector< Group* > groupVec;
1933           if ( grp._groups.empty() ) groupVec.push_back( & grp );
1934           else                       groupVec = grp._groups;
1935           for ( size_t iG = 0; iG < groupVec.size(); ++iG )
1936             {
1937               Group* aG = groupVec[ iG ];
1938               if ( (int)getDim( aG ) != dim )
1939                 continue;
1940               for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
1941                 cell2order.insert( cell2order.end(), make_pair( aG->_cells[iC], orderInGroup++ ));
1942             }
1943           if ( cell2order.empty() )
1944             continue;
1945           bool isSelfIntersect = ( orderInGroup != cell2order.size() );
1946           if ( isSelfIntersect ) // self intersecting group
1947             {
1948               ostringstream msg;
1949               msg << "Self intersecting sub-mesh: id = " << i+1
1950                   << ", name = |" << grp._name << "|" << endl
1951                   << " nb unique elements = " << cell2order.size() << endl
1952                   << " total nb elements  = " << orderInGroup;
1953               if ( grp._isProfile )
1954                 {
1955                   THROW_IK_EXCEPTION( msg.str() );
1956                 }
1957               else
1958                 {
1959                   cout << msg.str() << endl;
1960                 }
1961             }
1962           // create a med group
1963           grp._medGroup = DataArrayInt::New();
1964           grp._medGroup->setName( grp._name.c_str() );
1965           grp._medGroup->alloc( cell2order.size(), /*nbOfCompo=*/1 );
1966           int * idsPrt = grp._medGroup->getPointer();
1967           TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
1968           for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
1969             *idsPrt++ = (*cell2orderIt).first->_number - 1;
1970
1971           // try to set the mesh name
1972           if ( !isMeshNameSet &&
1973                dim == meshDim &&
1974                !grp._name.empty() &&
1975                grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
1976             {
1977               mesh->setName( grp._name.c_str() );
1978               isMeshNameSet = true;
1979             }
1980           if ( !grp._name.empty() )
1981             {
1982               medGroups.push_back( grp._medGroup );
1983             }
1984           // set relocation table
1985           setRelocationTable( &grp, cell2order );
1986
1987           // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
1988           // and several names (pile 27) refer (pile 10) to this group.
1989           // We create a copy of this group per each named reference
1990           for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
1991             if ( !grp._refNames[ iRef ].empty() )
1992               {
1993                 refGroups.push_back( grp._medGroup->deepCpy() );
1994                 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
1995                 medGroups.push_back( refGroups.back() );
1996               }
1997         }
1998       mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
1999     }
2000 }
2001
2002 //================================================================================
2003 /*!
2004  * \brief Return true if the group is on all elements and return its relative dimension
2005  */
2006 //================================================================================
2007
2008 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
2009 {
2010   int dim = getDim( grp );
2011
2012   int nbElems = 0;
2013   CellsByDimIterator dimCells( *this, dim );
2014   while ( const set<Cell > * cells = dimCells.nextType() )
2015     nbElems += cells->size();
2016
2017   const bool onAll = ( nbElems == grp->size() );
2018
2019   if ( dim == 0 )
2020     dimRel = 0;
2021   else
2022     {
2023       int meshDim = 3;
2024       for ( ; meshDim > 0; --meshDim )
2025         {
2026           dimCells.init( meshDim );
2027           if ( dimCells.nextType() )
2028             break;
2029         }
2030       dimRel = dim - meshDim;
2031     }
2032   return onAll;
2033 }
2034
2035 //================================================================================
2036 /*!
2037  * \brief Makes fields from own data
2038  */
2039 //================================================================================
2040
2041 ParaMEDMEM::MEDFileFields * IntermediateMED::makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh)
2042 {
2043   if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
2044
2045   // set long names
2046   set< string > usedFieldNames;
2047   setFieldLongNames(usedFieldNames);
2048
2049   MEDFileFields* fields = MEDFileFields::New();
2050
2051   for ( size_t i = 0; i < _nodeFields.size(); ++i )
2052     setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
2053
2054   for ( size_t i = 0; i < _cellFields.size(); ++i )
2055     setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
2056
2057   return fields;
2058 }
2059
2060 //================================================================================
2061 /*!
2062  * \brief Make med fields from a SauvUtilities::DoubleField
2063  */
2064 //================================================================================
2065
2066 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
2067                                  ParaMEDMEM::MEDFileFields*  medFields,
2068                                  ParaMEDMEM::MEDFileUMesh*   mesh,
2069                                  const TID                   castemID,
2070                                  set< string >&              usedFieldNames)
2071 {
2072   if ( !fld || !fld->isMedCompatible() ) return;
2073
2074   // if ( !fld->hasCommonSupport() ):
2075   //     each sub makes MEDFileFieldMultiTS
2076   // else:
2077   //     unite several subs into a MEDCouplingFieldDouble
2078
2079   const bool uniteSubs = fld->hasCommonSupport();
2080   if ( !uniteSubs )
2081     cout << "Castem field #" << castemID << " " << fld->_name
2082          << " is incompatible with MED format, so we split it into several fields" << endl;
2083
2084   for ( size_t iSub = 0; iSub < fld->_sub.size(); )
2085     {
2086       // set field name
2087       if ( !uniteSubs || fld->_name.empty() )
2088         makeFieldNewName( usedFieldNames, fld );
2089
2090       // allocate values
2091       DataArrayDouble * values = DataArrayDouble::New();
2092       values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
2093
2094       // set values
2095       double * valPtr = values->getPointer();
2096       if ( uniteSubs )
2097         {
2098           int nbElems = fld->_group->size();
2099           for ( int elemShift = 0; elemShift < nbElems; )
2100             elemShift += fld->setValues( valPtr, iSub++, elemShift );
2101           setTS( fld, values, medFields, mesh );
2102         }
2103       else
2104         {
2105           fld->setValues( valPtr, iSub++ );
2106           setTS( fld, values, medFields, mesh, iSub );
2107         }
2108     }
2109 }
2110
2111 //================================================================================
2112 /*!
2113  * \brief Store value array of a field into med fields
2114  */
2115 //================================================================================
2116
2117 void IntermediateMED::setTS( SauvUtilities::DoubleField*  fld,
2118                              ParaMEDMEM::DataArrayDouble* values,
2119                              ParaMEDMEM::MEDFileFields*   medFields,
2120                              ParaMEDMEM::MEDFileUMesh*    mesh,
2121                              const int                    iSub)
2122 {
2123   // analyze a field support
2124   const Group* support = fld->getSupport();
2125   int dimRel;
2126   const bool onAll = isOnAll( support, dimRel );
2127   if ( !onAll && support->_name.empty() )
2128     {
2129       const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
2130       support->_medGroup->setName( support->_name.c_str() );
2131     }
2132
2133   // make a time-stamp
2134   MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType(),
2135                                                                     fld->getMedTimeDisc() );
2136   timeStamp->setName( fld->_name.c_str() );
2137   timeStamp->setDescription( fld->_description.c_str() );
2138   MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
2139   timeStamp->setMesh( dimMesh );
2140   for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
2141     values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
2142   timeStamp->setArray( values );
2143   values->decrRef();
2144
2145   // get a field to add the time-stamp
2146   bool isNewMedField = false;
2147   if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
2148     {
2149       fld->_curMedField = MEDFileFieldMultiTS::New();
2150       isNewMedField = true;
2151     }
2152
2153   // set an order
2154   timeStamp->setOrder( fld->_curMedField->getNumberOfTS() );
2155
2156   // add the time-stamp
2157   if ( onAll )
2158     fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
2159   else
2160     fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
2161   timeStamp->decrRef();
2162
2163   if ( isNewMedField ) // timeStamp must be added before this
2164     {
2165       medFields->pushField( fld->_curMedField );
2166     }
2167 }
2168
2169 //================================================================================
2170 /*!
2171  * \brief Make a new unique name for a field
2172  */
2173 //================================================================================
2174
2175 void IntermediateMED::makeFieldNewName(std::set< std::string >&    usedNames,
2176                                        SauvUtilities::DoubleField* fld )
2177 {
2178   string base = fld->_name;
2179   if ( base.empty() )
2180     {
2181       base = "F_";
2182     }
2183   else
2184     {
2185       string::size_type pos = base.rfind('_');
2186       if ( pos != string::npos )
2187         base = base.substr( 0, pos+1 );
2188       else
2189         base += '_';
2190     }
2191
2192   int i = 1;
2193   do
2194     {
2195       fld->_name = base + SauvUtilities::toString( i++ );
2196     }
2197   while( !usedNames.insert( fld->_name ).second );
2198 }
2199
2200 //================================================================================
2201 /*!
2202  * \brief Return a vector ready to fill in
2203  */
2204 //================================================================================
2205
2206 std::vector< double >& DoubleField::addComponent( int nb_values )
2207 {
2208   _comp_values.push_back( std::vector< double >() );
2209   std::vector< double >& res = _comp_values.back();
2210   res.resize( nb_values );
2211   return res;
2212 }
2213
2214 DoubleField::~DoubleField()
2215 {
2216   if(_curMedField)
2217     _curMedField->decrRef();
2218 }
2219
2220 //================================================================================
2221 /*!
2222  * \brief Returns a supporting group
2223  */
2224 //================================================================================
2225
2226 const Group* DoubleField::getSupport( const int iSub ) const
2227 {
2228   return _group ? _group : _sub[iSub]._support;
2229 }
2230
2231 //================================================================================
2232 /*!
2233  * \brief Return true if each sub-component is a time stamp
2234  */
2235 //================================================================================
2236
2237 bool DoubleField::isMultiTimeStamps() const
2238 {
2239   if ( _sub.size() < 2 )
2240     return false;
2241   bool sameSupports = true;
2242   Group* grp1 = _sub[0]._support;
2243   for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
2244     sameSupports = ( grp1 == _sub[i]._support );
2245
2246   return sameSupports;
2247 }
2248
2249 //================================================================================
2250 /*!
2251  * \brief True if the field can be converted into the med field
2252  */
2253 //================================================================================
2254
2255 bool DoubleField::isMedCompatible() const
2256 {
2257   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
2258     {
2259       if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
2260         THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
2261
2262       if ( !_sub[iSub].isValidNbGauss() )
2263         {
2264           cout << "Skip field <" << _name << "> : different nb of gauss points in components" <<endl;
2265           return false;
2266         }
2267     }
2268   // check if there are no gauss or nbGauss() == nbCellNodes,
2269   // else we lack info on gauss point localization
2270   // TODO?
2271   return true;
2272 }
2273
2274 //================================================================================
2275 /*!
2276  * \brief return true if all sub-components has same components and same nbGauss
2277  */
2278 //================================================================================
2279
2280 bool DoubleField::hasSameComponentsBySupport() const
2281 {
2282   vector< _Sub_data >::const_iterator sub_data = _sub.begin();
2283   const _Sub_data& first_sub_data = *sub_data;
2284   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
2285   {
2286     if ( first_sub_data._comp_names != sub_data->_comp_names )
2287       return false; // diff names of components
2288
2289     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
2290          first_sub_data._support->_cellType == sub_data->_support->_cellType)
2291       return false; // diff nb of gauss points on same cell type
2292   }
2293   return true;
2294 }
2295
2296 //================================================================================
2297 /*!
2298  * \brief Return type of MEDCouplingFieldDouble
2299  */
2300 //================================================================================
2301
2302 ParaMEDMEM::TypeOfField DoubleField::getMedType( const int iSub ) const
2303 {
2304   using namespace INTERP_KERNEL;
2305
2306   const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
2307   if ( _sub[iSub].nbGauss() > 1 )
2308     {
2309       const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
2310       return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
2311     }
2312   else
2313     {
2314       return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
2315     }
2316 }
2317
2318 //================================================================================
2319 /*!
2320  * \brief Return TypeOfTimeDiscretization
2321  */
2322 //================================================================================
2323
2324 ParaMEDMEM::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
2325 {
2326   return ONE_TIME;
2327   // NO_TIME = 4,
2328   // ONE_TIME = 5,
2329   // LINEAR_TIME = 6,
2330   // CONST_ON_TIME_INTERVAL = 7
2331 }
2332
2333 //================================================================================
2334 /*!
2335  * \brief Return nb tuples to be used to allocate DataArrayDouble
2336  */
2337 //================================================================================
2338
2339 int DoubleField::getNbTuples( const int iSub ) const
2340 {
2341   int nb = 0;
2342   if ( hasCommonSupport() && !_group->_groups.empty() )
2343     for ( size_t i = 0; i < _group->_groups.size(); ++i )
2344       nb += _sub[i].nbGauss() * _sub[i]._support->size();
2345   else
2346     nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
2347   return nb;
2348 }
2349
2350 //================================================================================
2351 /*!
2352  * \brief Store values of a sub-component and return nb of elements in the iSub
2353  */
2354 //================================================================================
2355
2356 int DoubleField::setValues( double * valPtr, const int iSub, const int elemShift ) const
2357 {
2358   // find values for iSub
2359   int iComp = 0;
2360   for ( int iS = 0; iS < iSub; ++iS )
2361     iComp += _sub[iS].nbComponents();
2362   const vector< double > * compValues = &_comp_values[ iComp ];
2363
2364   const vector< unsigned >& relocTable = getSupport( iSub )->_relocTable;
2365
2366   // Set values
2367
2368   const int nbElems      = _sub[iSub]._support->size();
2369   const int nbGauss      = _sub[iSub].nbGauss();
2370   const int nbComponents = _sub[iSub].nbComponents();
2371   const int nbValsByElem = nbComponents * nbGauss;
2372   // check nb values
2373   int nbVals = 0;
2374   for ( iComp = 0; iComp < nbComponents; ++iComp )
2375     nbVals += compValues[iComp].size();
2376   if ( nbVals != nbElems * nbValsByElem )
2377     THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
2378   // compute nb values in previous subs
2379   int valsShift = 0;
2380   for ( int iS = iSub-1, shift = elemShift; shift > 0; )
2381   {
2382     int nbE = _sub[iS]._support->size();
2383     shift -= nbE;
2384     valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
2385   }
2386   for ( int iE = 0; iE < nbElems; ++iE )
2387     {
2388       int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
2389       for ( iComp = 0; iComp < nbComponents; ++iComp )
2390         for ( int iG = 0; iG < nbGauss; ++iG )
2391           valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
2392     }
2393   return nbElems;
2394 }
2395
2396 //================================================================================
2397 /*!
2398  * \brief Destructor of IntermediateMED
2399  */
2400 //================================================================================
2401
2402 IntermediateMED::~IntermediateMED()
2403 {
2404   for ( size_t i = 0; i < _nodeFields.size(); ++i )
2405     if ( _nodeFields[i] )
2406       delete _nodeFields[i];
2407   _nodeFields.clear();
2408
2409   for ( size_t i = 0; i < _cellFields.size(); ++i )
2410     if ( _cellFields[i] )
2411       delete _cellFields[i];
2412   _cellFields.clear();
2413
2414   for ( size_t i = 0; i < _groups.size(); ++i )
2415     if ( _groups[i]._medGroup )
2416       _groups[i]._medGroup->decrRef();
2417 }
2418
2419 //================================================================================
2420 /*!
2421  * \brief CellsByDimIterator constructor
2422  */
2423 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
2424 {
2425   myImed = & medi;
2426   init( dimm );
2427 }
2428 /*!
2429  * \brief Initialize iteration on cells of given dimention
2430  */
2431 void CellsByDimIterator::init(const int  dimm)
2432 {
2433   myCurType = -1;
2434   myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
2435   myDim = dimm;
2436 }
2437 /*!
2438  * \brief return next set of Cell's of required dimension
2439  */
2440 const std::set< Cell > * CellsByDimIterator::nextType()
2441 {
2442   while ( ++myCurType < myTypeEnd )
2443     if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
2444       return & myImed->_cellsByType[myCurType];
2445   return 0;
2446 }
2447 /*!
2448  * \brief return dimension of cells returned by the last or further next()
2449  */
2450 int CellsByDimIterator::dim(const bool last) const
2451 {
2452   int typp = myCurType;
2453   if ( !last )
2454     while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
2455       ++typp;
2456   return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
2457 }
2458 // END CellsByDimIterator ========================================================
2459