]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/MEDLoader/SauvMedConvertor.cxx
Salome HOME
XDR detection. The aim is to enable it even on win7.
[tools/medcoupling.git] / src / MEDLoader / SauvMedConvertor.cxx
1 // Copyright (C) 2007-2013  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 #else
45 #include <unistd.h>
46 #endif
47
48 #ifdef HAS_XDR
49 #include <rpc/xdr.h>
50 #endif
51
52 using namespace SauvUtilities;
53 using namespace ParaMEDMEM;
54 using namespace std;
55
56 namespace
57 {
58   // for ASCII file reader
59   const int GIBI_MaxOutputLen = 150; // max length of a line in the sauve file
60   const int GIBI_BufferSize   = 16184; // for buffered reading
61
62   using namespace INTERP_KERNEL;
63
64   const size_t MaxMedCellType = NORM_ERROR;
65   const size_t NbGibiCellTypes = 47;
66   const TCellType GibiTypeToMed[NbGibiCellTypes] =
67     {
68       /*1 */ NORM_POINT1 ,/*2 */ NORM_SEG2   ,/*3 */ NORM_SEG3   ,/*4 */ NORM_TRI3   ,/*5 */ NORM_ERROR  ,
69       /*6 */ NORM_TRI6   ,/*7 */ NORM_ERROR  ,/*8 */ NORM_QUAD4  ,/*9 */ NORM_ERROR  ,/*10*/ NORM_QUAD8  ,
70       /*11*/ NORM_ERROR  ,/*12*/ NORM_ERROR  ,/*13*/ NORM_ERROR  ,/*14*/ NORM_HEXA8  ,/*15*/ NORM_HEXA20 ,
71       /*16*/ NORM_PENTA6 ,/*17*/ NORM_PENTA15,/*18*/ NORM_ERROR  ,/*19*/ NORM_ERROR  ,/*20*/ NORM_ERROR  ,
72       /*21*/ NORM_ERROR  ,/*22*/ NORM_ERROR  ,/*23*/ NORM_TETRA4 ,/*24*/ NORM_TETRA10,/*25*/ NORM_PYRA5  ,
73       /*26*/ NORM_PYRA13 ,/*27*/ NORM_ERROR  ,/*28*/ NORM_ERROR  ,/*29*/ NORM_ERROR  ,/*30*/ NORM_ERROR  ,
74       /*31*/ NORM_ERROR  ,/*32*/ NORM_ERROR  ,/*33*/ NORM_ERROR  ,/*34*/ NORM_ERROR  ,/*35*/ NORM_ERROR  ,
75       /*36*/ NORM_ERROR  ,/*37*/ NORM_ERROR  ,/*38*/ NORM_ERROR  ,/*39*/ NORM_ERROR  ,/*40*/ NORM_ERROR  ,
76       /*41*/ NORM_ERROR  ,/*42*/ NORM_ERROR  ,/*43*/ NORM_ERROR  ,/*44*/ NORM_ERROR  ,/*45*/ NORM_ERROR  ,
77       /*46*/ NORM_ERROR  ,/*47*/ NORM_ERROR
78     };
79
80   //================================================================================
81   /*!
82    * \brief Return dimension of a group
83    */
84   //================================================================================
85
86   unsigned getDim( const Group* grp )
87   {
88     return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
89   }
90
91   //================================================================================
92   /*!
93    * \brief Converts connectivity of quadratic elements
94    */
95   //================================================================================
96
97   inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
98                                 const Cell &                            aCell )
99   {
100     if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
101       {
102         Cell* ma = const_cast<Cell*>(&aCell);
103         vector< Node* > new_nodes( ma->_nodes.size() );
104         for ( size_t i = 0; i < new_nodes.size(); ++i )
105           new_nodes[ i ] = ma->_nodes[ conn[ i ]];
106         ma->_nodes.swap( new_nodes );
107       }
108   }
109
110   //================================================================================
111   /*!
112    * \brief Returns a vector of pairs of node indices to inverse a med volume element
113    */
114   //================================================================================
115
116   void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
117                          vector<pair<int,int> > &                swapVec )
118   {
119     swapVec.clear();
120
121     switch ( type )
122       {
123       case NORM_TETRA4:
124         swapVec.resize(1);
125         swapVec[0] = make_pair( 1, 2 );
126         break;
127       case NORM_PYRA5:
128         swapVec.resize(1);
129         swapVec[0] = make_pair( 1, 3 );
130         break;
131       case NORM_PENTA6:
132         swapVec.resize(2);
133         swapVec[0] = make_pair( 1, 2 );
134         swapVec[1] = make_pair( 4, 5 );
135         break;
136       case NORM_HEXA8:
137         swapVec.resize(2);
138         swapVec[0] = make_pair( 1, 3 );
139         swapVec[1] = make_pair( 5, 7 );
140         break;
141       case NORM_TETRA10:
142         swapVec.resize(3);
143         swapVec[0] = make_pair( 1, 2 );
144         swapVec[1] = make_pair( 4, 6 );
145         swapVec[2] = make_pair( 8, 9 );
146         break;
147       case NORM_PYRA13:
148         swapVec.resize(4);
149         swapVec[0] = make_pair( 1, 3 );
150         swapVec[1] = make_pair( 5, 8 );
151         swapVec[2] = make_pair( 6, 7 );
152         swapVec[3] = make_pair( 10, 12 );
153         break;
154       case NORM_PENTA15:
155         swapVec.resize(4);
156         swapVec[0] = make_pair( 1, 2 );
157         swapVec[1] = make_pair( 4, 5 );
158         swapVec[2] = make_pair( 6, 8 );
159         swapVec[3] = make_pair( 9, 11 );
160         break;
161       case NORM_HEXA20:
162         swapVec.resize(7);
163         swapVec[0] = make_pair( 1, 3 );
164         swapVec[1] = make_pair( 5, 7 );
165         swapVec[2] = make_pair( 8, 11 );
166         swapVec[3] = make_pair( 9, 10 );
167         swapVec[4] = make_pair( 12, 15 );
168         swapVec[5] = make_pair( 13, 14 );
169         swapVec[6] = make_pair( 17, 19 );
170         break;
171         //   case NORM_SEG3: no need to reverse edges
172         //     swapVec.resize(1);
173         //     swapVec[0] = make_pair( 1, 2 );
174         //     break;
175       case NORM_TRI6:
176         swapVec.resize(2);
177         swapVec[0] = make_pair( 1, 2 );
178         swapVec[1] = make_pair( 3, 5 );
179         break;
180       case NORM_QUAD8:
181         swapVec.resize(3);
182         swapVec[0] = make_pair( 1, 3 );
183         swapVec[1] = make_pair( 4, 7 );
184         swapVec[2] = make_pair( 5, 6 );
185         break;
186       default:;
187       }
188   }
189
190   //================================================================================
191   /*!
192    * \brief Inverses element orientation using vector of indices to swap
193    */
194   //================================================================================
195
196   inline void reverse(const Cell & aCell, const vector<pair<int,int> > & swapVec )
197   {
198     Cell* ma = const_cast<Cell*>(&aCell);
199     for ( unsigned i = 0; i < swapVec.size(); ++i )
200       std::swap( ma->_nodes[ swapVec[i].first ],
201                  ma->_nodes[ swapVec[i].second ]);
202     if ( swapVec.empty() )
203       ma->_reverse = true;
204     else
205       ma->_reverse = false;
206   }
207   //================================================================================
208   /*!
209    * \brief Comparator of cells by number used for ordering cells within a med group
210    */
211   struct TCellByIDCompare
212   {
213     bool operator () (const Cell* i1, const Cell* i2)
214     {
215       return i1->_number < i2->_number;
216     }
217   };
218   typedef map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
219
220   //================================================================================
221   /*!
222    * \brief Fill Group::_relocTable if necessary
223    */
224   //================================================================================
225
226   void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
227   {
228     if ( !grp->_isProfile ) return;
229
230     // check if relocation table is necessary
231     bool isRelocated = false;
232     unsigned newOrder = 0;
233     TCellToOrderMap::iterator c2oIt = cell2order.begin(), c2oEnd = cell2order.end();
234     for ( ; !isRelocated && c2oIt != c2oEnd; ++c2oIt, ++newOrder )
235       isRelocated = ( c2oIt->second != newOrder );
236
237     if ( isRelocated )
238     {
239       grp->_relocTable.resize( cell2order.size() );
240       for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
241         grp->_relocTable[ c2oIt->second ] = newOrder;
242     }
243   }
244 }
245
246 namespace // define default GAUSS points
247 {
248   typedef std::vector<double> TDoubleVector;
249   typedef double*             TCoordSlice;
250   typedef int                 TInt;
251   //---------------------------------------------------------------
252   //! Shape function definitions
253   //---------------------------------------------------------------
254   struct TShapeFun
255   {
256     TInt myDim;
257     TInt myNbRef;
258     TDoubleVector myRefCoord;
259
260     TShapeFun(TInt theDim = 0, TInt theNbRef = 0)
261       :myDim(theDim),myNbRef(theNbRef),myRefCoord(theNbRef*theDim) {}
262
263     TInt GetNbRef() const { return myNbRef; }
264
265     TCoordSlice GetCoord(TInt theRefId) { return &myRefCoord[0] + theRefId * myDim; }
266   };
267   //---------------------------------------------------------------
268   /*!
269    * \brief Description of family of integration points
270    */
271   //---------------------------------------------------------------
272   struct TGaussDef
273   {
274     int           myType;      //!< element geometry (EGeometrieElement or med_geometrie_element)
275     TDoubleVector myRefCoords; //!< description of reference points
276     TDoubleVector myCoords;    //!< coordinates of Gauss points
277     TDoubleVector myWeights;   //!< weights, len(weights)==<nb of gauss points>
278
279     /*!
280      * \brief Creates definition of gauss points family
281      *  \param geomType - element geometry (EGeometrieElement or med_geometrie_element)
282      *  \param nbPoints - nb gauss point
283      *  \param variant - [1-3] to choose the variant of definition
284      * 
285      * Throws in case of invalid parameters
286      * variant == 1 refers to "Fonctions de forme et points d'integration 
287      *              des elements finis" v7.4 by J. PELLET, X. DESROCHES, 15/09/05
288      * variant == 2 refers to the same doc v6.4 by J.P. LEFEBVRE, X. DESROCHES, 03/07/03
289      * variant == 3 refers to the same doc v6.4, second variant for 2D elements
290      */
291     TGaussDef(const int geomType, const int nbPoints, const int variant=1);
292
293     int dim() const { return SauvUtilities::getDimension( NormalizedCellType( myType )); }
294     int nbPoints() const { return myWeights.capacity(); }
295
296   private:
297     void add(const double x, const double weight);
298     void add(const double x, const double y, const double weight);
299     void add(const double x, const double y, const double z, const double weight);
300     void setRefCoords(const TShapeFun& aShapeFun) { myRefCoords = aShapeFun.myRefCoord; }
301   };
302   struct TSeg2a: TShapeFun {
303     TSeg2a();
304   };
305   struct TSeg3a: TShapeFun {
306     TSeg3a();
307   };
308   struct TTria3a: TShapeFun {
309     TTria3a();
310   };
311   struct TTria6a: TShapeFun {
312     TTria6a();
313   };
314   struct TTria3b: TShapeFun {
315     TTria3b();
316   };
317   struct TTria6b: TShapeFun {
318     TTria6b();
319   };
320   struct TQuad4a: TShapeFun {
321     TQuad4a();
322   };
323   struct TQuad8a: TShapeFun {
324     TQuad8a();
325   };
326   struct TQuad4b: TShapeFun {
327     TQuad4b();
328   };
329   struct TQuad8b: TShapeFun {
330     TQuad8b();
331   };
332   struct TTetra4a: TShapeFun {
333     TTetra4a();
334   };
335   struct TTetra10a: TShapeFun {
336     TTetra10a();
337   };
338   struct TTetra4b: TShapeFun {
339     TTetra4b();
340   };
341   struct TTetra10b: TShapeFun {
342     TTetra10b();
343   };
344   struct THexa8a: TShapeFun {
345     THexa8a();
346   };
347   struct THexa20a: TShapeFun {
348     THexa20a(TInt theDim = 3, TInt theNbRef = 20);
349   };
350   struct THexa27a: THexa20a {
351     THexa27a();
352   };
353   struct THexa8b: TShapeFun {
354     THexa8b();
355   };
356   struct THexa20b: TShapeFun {
357     THexa20b(TInt theDim = 3, TInt theNbRef = 20);
358   };
359   struct TPenta6a: TShapeFun {
360     TPenta6a();
361   };
362   struct TPenta6b: TShapeFun {
363     TPenta6b();
364   };
365   struct TPenta15a: TShapeFun {
366     TPenta15a();
367   };
368   struct TPenta15b: TShapeFun {
369     TPenta15b();
370   };
371   struct TPyra5a: TShapeFun {
372     TPyra5a();
373   };
374   struct TPyra5b: TShapeFun {
375     TPyra5b();
376   };
377   struct TPyra13a: TShapeFun {
378     TPyra13a();
379   };
380   struct TPyra13b: TShapeFun {
381     TPyra13b();
382   };
383
384   void TGaussDef::add(const double x, const double weight)
385   {
386     if ( dim() != 1 )
387       THROW_IK_EXCEPTION("TGaussDef: dim() != 1");
388     if ( myWeights.capacity() == myWeights.size() )
389       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
390     myCoords.push_back( x );
391     myWeights.push_back( weight );
392   }
393   void TGaussDef::add(const double x, const double y, const double weight)
394   {
395     if ( dim() != 2 )
396       THROW_IK_EXCEPTION("TGaussDef: dim() != 2");
397     if ( myWeights.capacity() == myWeights.size() )
398       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
399     myCoords.push_back( x );
400     myCoords.push_back( y );
401     myWeights.push_back( weight );
402   }
403   void TGaussDef::add(const double x, const double y, const double z, const double weight)
404   {
405     if ( dim() != 3 )
406       THROW_IK_EXCEPTION("TGaussDef: dim() != 3");
407     if ( myWeights.capacity() == myWeights.size() )
408       THROW_IK_EXCEPTION("TGaussDef: Extra gauss point");
409     myCoords.push_back( x );
410     myCoords.push_back( y );
411     myCoords.push_back( z );
412     myWeights.push_back( weight );
413   }
414
415   //---------------------------------------------------------------
416   TSeg2a::TSeg2a():TShapeFun(1,2)
417   {
418     TInt aNbRef = GetNbRef();
419     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
420       TCoordSlice aCoord = GetCoord(aRefId);
421       switch(aRefId){
422       case  0: aCoord[0] = -1.0; break;
423       case  1: aCoord[0] =  1.0; break;
424       }
425     }
426   }
427   //---------------------------------------------------------------
428   TSeg3a::TSeg3a():TShapeFun(1,3)
429   {
430     TInt aNbRef = GetNbRef();
431     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
432       TCoordSlice aCoord = GetCoord(aRefId);
433       switch(aRefId){
434       case  0: aCoord[0] = -1.0; break;
435       case  1: aCoord[0] =  1.0; break;
436       case  2: aCoord[0] =  0.0; break;
437       }
438     }
439   }
440   //---------------------------------------------------------------
441   TTria3a::TTria3a():
442     TShapeFun(2,3)
443   {
444     TInt aNbRef = GetNbRef();
445     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
446       TCoordSlice aCoord = GetCoord(aRefId);
447       switch(aRefId){
448       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
449       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
450       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
451       }
452     }
453   }
454   //---------------------------------------------------------------
455   TTria6a::TTria6a():TShapeFun(2,6)
456   {
457     TInt aNbRef = GetNbRef();
458     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
459       TCoordSlice aCoord = GetCoord(aRefId);
460       switch(aRefId){
461       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
462       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
463       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
464
465       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
466       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
467       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
468       }
469     }
470   }
471   //---------------------------------------------------------------
472   TTria3b::TTria3b():
473     TShapeFun(2,3)
474   {
475     TInt aNbRef = GetNbRef();
476     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
477       TCoordSlice aCoord = GetCoord(aRefId);
478       switch(aRefId){
479       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
480       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
481       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
482       }
483     }
484   }
485   //---------------------------------------------------------------
486   TTria6b::TTria6b():
487     TShapeFun(2,6)
488   {
489     TInt aNbRef = GetNbRef();
490     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
491       TCoordSlice aCoord = GetCoord(aRefId);
492       switch(aRefId){
493       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
494       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
495       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
496
497       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
498       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
499       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
500       }
501     }
502   }
503   //---------------------------------------------------------------
504   TQuad4a::TQuad4a():
505     TShapeFun(2,4)
506   {
507     TInt aNbRef = GetNbRef();
508     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
509       TCoordSlice aCoord = GetCoord(aRefId);
510       switch(aRefId){
511       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
512       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
513       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
514       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
515       }
516     }
517   }
518   //---------------------------------------------------------------
519   TQuad8a::TQuad8a():
520     TShapeFun(2,8)
521   {
522     TInt aNbRef = GetNbRef();
523     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
524       TCoordSlice aCoord = GetCoord(aRefId);
525       switch(aRefId){
526       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
527       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
528       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
529       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
530
531       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
532       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
533       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
534       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
535       }
536     }
537   }
538   //---------------------------------------------------------------
539   TQuad4b::TQuad4b():
540     TShapeFun(2,4)
541   {
542     TInt aNbRef = GetNbRef();
543     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
544       TCoordSlice aCoord = GetCoord(aRefId);
545       switch(aRefId){
546       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
547       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
548       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
549       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
550       }
551     }
552   }
553   //---------------------------------------------------------------
554   TQuad8b::TQuad8b():
555     TShapeFun(2,8)
556   {
557     TInt aNbRef = GetNbRef();
558     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
559       TCoordSlice aCoord = GetCoord(aRefId);
560       switch(aRefId){
561       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
562       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
563       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
564       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
565
566       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
567       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
568       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
569       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
570       }
571     }
572   }
573   //---------------------------------------------------------------
574   TTetra4a::TTetra4a():
575     TShapeFun(3,4)
576   {
577     TInt aNbRef = GetNbRef();
578     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
579       TCoordSlice aCoord = GetCoord(aRefId);
580       switch(aRefId){
581       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
582       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
583       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
584       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
585       }
586     }
587   }
588   //---------------------------------------------------------------
589   TTetra10a::TTetra10a():
590     TShapeFun(3,10)
591   {
592     TInt aNbRef = GetNbRef();
593     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
594       TCoordSlice aCoord = GetCoord(aRefId);
595       switch(aRefId){
596       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
597       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
598       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
599       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
600
601       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
602       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
603       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
604
605       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
606       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
607       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
608       }
609     }
610   }
611   //---------------------------------------------------------------
612   TTetra4b::TTetra4b():
613     TShapeFun(3,4)
614   {
615     TInt aNbRef = GetNbRef();
616     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
617       TCoordSlice aCoord = GetCoord(aRefId);
618       switch(aRefId){
619       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
620       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
621       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
622       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
623       }
624     }
625   }
626   //---------------------------------------------------------------
627   TTetra10b::TTetra10b():
628     TShapeFun(3,10)
629   {
630     TInt aNbRef = GetNbRef();
631     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
632       TCoordSlice aCoord = GetCoord(aRefId);
633       switch(aRefId){
634       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
635       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
636       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
637       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
638
639       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
640       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
641       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
642
643       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
644       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
645       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
646       }
647     }
648   }
649   //---------------------------------------------------------------
650   THexa8a::THexa8a():
651     TShapeFun(3,8)
652   {
653     TInt aNbRef = GetNbRef();
654     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
655       TCoordSlice aCoord = GetCoord(aRefId);
656       switch(aRefId){
657       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
658       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
659       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
660       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
661       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
662       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
663       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
664       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
665       }
666     }
667   }
668   //---------------------------------------------------------------
669   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
670     TShapeFun(theDim,theNbRef)
671   {
672     TInt aNbRef = myRefCoord.size();
673     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
674       TCoordSlice aCoord = GetCoord(aRefId);
675       switch(aRefId){
676       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
677       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
678       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
679       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
680       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
681       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
682       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
683       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
684
685       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
686       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
687       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
688       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
689       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
690       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
691       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
692       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
693       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
694       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
695       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
696       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
697       }
698     }
699   }
700   //---------------------------------------------------------------
701   THexa27a::THexa27a():
702     THexa20a(3,27)
703   {
704     TInt aNbRef = myRefCoord.size();
705     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
706       TCoordSlice aCoord = GetCoord(aRefId);
707       switch(aRefId){
708       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
709       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
710       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
711       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
712       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
713       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
714       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
715       }
716     }
717   }
718   //---------------------------------------------------------------
719   THexa8b::THexa8b():
720     TShapeFun(3,8)
721   {
722     TInt aNbRef = GetNbRef();
723     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
724       TCoordSlice aCoord = GetCoord(aRefId);
725       switch(aRefId){
726       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
727       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
728       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
729       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
730       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
731       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
732       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
733       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
734       }
735     }
736   }
737   //---------------------------------------------------------------
738   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
739     TShapeFun(theDim,theNbRef)
740   {
741     TInt aNbRef = myRefCoord.size();
742     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
743       TCoordSlice aCoord = GetCoord(aRefId);
744       switch(aRefId){
745       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
746       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
747       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
748       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
749       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
750       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
751       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
752       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
753
754       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
755       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
756       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
757       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
758       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
759       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
760       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
761       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
762       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
763       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
764       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
765       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
766       }
767     }
768   }
769   //---------------------------------------------------------------
770   TPenta6a::TPenta6a():
771     TShapeFun(3,6)
772   {
773     TInt aNbRef = myRefCoord.size();
774     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
775       TCoordSlice aCoord = GetCoord(aRefId);
776       switch(aRefId){
777       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
778       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
779       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
780       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
781       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
782       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
783       }
784     }
785   }
786   //---------------------------------------------------------------
787   TPenta6b::TPenta6b():
788     TShapeFun(3,6)
789   {
790     TInt aNbRef = myRefCoord.size();
791     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
792       TCoordSlice aCoord = GetCoord(aRefId);
793       switch(aRefId){
794       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
795       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
796       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
797       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
798       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
799       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
800       }
801     }
802   }
803   //---------------------------------------------------------------
804   TPenta15a::TPenta15a():
805     TShapeFun(3,15)
806   {
807     TInt aNbRef = myRefCoord.size();
808     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
809       TCoordSlice aCoord = GetCoord(aRefId);
810       switch(aRefId){
811       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
812       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
813       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
814       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
815       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
816       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
817
818       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
819       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
820       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
821       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
822       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
823       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
824       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
825       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
826       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
827       }
828     }
829   }
830   //---------------------------------------------------------------
831   TPenta15b::TPenta15b():
832     TShapeFun(3,15)
833   {
834     TInt aNbRef = myRefCoord.size();
835     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
836       TCoordSlice aCoord = GetCoord(aRefId);
837       switch(aRefId){
838       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
839       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
840       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
841       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
842       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
843       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
844
845       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
846       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
847       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
848       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
849       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
850       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
851       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
852       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
853       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
854       }
855     }
856   }
857   //---------------------------------------------------------------
858   TPyra5a::TPyra5a():
859     TShapeFun(3,5)
860   {
861     TInt aNbRef = myRefCoord.size();
862     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
863       TCoordSlice aCoord = GetCoord(aRefId);
864       switch(aRefId){
865       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
866       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
867       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
868       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
869       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
870       }
871     }
872   }
873   //---------------------------------------------------------------
874   TPyra5b::TPyra5b():
875     TShapeFun(3,5)
876   {
877     TInt aNbRef = myRefCoord.size();
878     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
879       TCoordSlice aCoord = GetCoord(aRefId);
880       switch(aRefId){        
881       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
882       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
883       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
884       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
885       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
886       }
887     }
888   }
889   //---------------------------------------------------------------
890   TPyra13a::TPyra13a():
891     TShapeFun(3,13)
892   {
893     TInt aNbRef = myRefCoord.size();
894     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
895       TCoordSlice aCoord = GetCoord(aRefId);
896       switch(aRefId){
897       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
898       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
899       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
900       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
901       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
902
903       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
904       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
905       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
906       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
907       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
908       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
909       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
910       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
911       }
912     }
913   }
914   //---------------------------------------------------------------
915   TPyra13b::TPyra13b():
916     TShapeFun(3,13)
917   {
918     TInt aNbRef = myRefCoord.size();
919     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
920       TCoordSlice aCoord = GetCoord(aRefId);
921       switch(aRefId){
922       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
923       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
924       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
925       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
926       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
927
928       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
929       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
930       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
931       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
932       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
933       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
934       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
935       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
936       }
937     }
938   }
939   /*!
940    * \brief Fill definition of gauss points family
941    */
942
943   TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
944   {
945     myType = geom;
946     myCoords .reserve( nbGauss * dim() );
947     myWeights.reserve( nbGauss );
948
949     switch ( geom ) {
950
951     case NORM_SEG2:
952     case NORM_SEG3:
953       if (geom == NORM_SEG2) setRefCoords( TSeg2a() );
954       else                   setRefCoords( TSeg3a() );
955       switch ( nbGauss ) {
956       case 1: {
957         add( 0.0, 2.0 ); break;
958       }
959       case 2: {
960         const double a = 0.577350269189626;
961         add(  a,  1.0 );
962         add( -a,  1.0 ); break;
963       }
964       case 3: {
965         const double a = 0.774596669241;
966         const double P1 = 1./1.8;
967         const double P2 = 1./1.125;
968         add( -a,  P1 );
969         add(  0,  P2 ); 
970         add(  a,  P1 ); break;
971       }
972       case 4: {
973         const double a  = 0.339981043584856, b  = 0.861136311594053;
974         const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
975         add(  a,  P1 );
976         add( -a,  P1 );
977         add(  b,  P2 ); 
978         add( -b,  P2 ); break;
979       }
980       default:
981         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for SEG"<<nbGauss);
982       }
983       break;
984
985     case NORM_TRI3:
986     case NORM_TRI6:
987       if ( variant == 1 ) {
988         if (geom == NORM_TRI3) setRefCoords( TTria3b() );
989         else                   setRefCoords( TTria6b() );
990         switch ( nbGauss ) {
991         case 1: { // FPG1
992           add( 1/3., 1/3., 1/2. ); break;
993         }
994         case 3: { // FPG3
995           // what about COT3 ???
996           add( 1/6., 1/6., 1/6. );
997           add( 2/3., 1/6., 1/6. );
998           add( 1/6., 2/3., 1/6. ); break;
999         }
1000         case 4: { // FPG4
1001           add( 1/5., 1/5.,  25/(24*4.) );
1002           add( 3/5., 1/5.,  25/(24*4.) );
1003           add( 1/5., 3/5.,  25/(24*4.) );
1004           add( 1/3., 1/3., -27/(24*4.) ); break;
1005         }
1006         case 6: { // FPG6
1007           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
1008           const double a  = 0.445948490915965, b = 0.091576213509771;
1009           add(     b,     b, P2 ); 
1010           add( 1-2*b,     b, P2 );
1011           add(     b, 1-2*b, P2 );
1012           add(     a, 1-2*a, P1 );
1013           add(     a,     a, P1 ); 
1014           add( 1-2*a,     a, P1 ); break;
1015         }
1016         case 7: { // FPG7
1017           const double A  = 0.470142064105115;
1018           const double B  = 0.101286507323456;
1019           const double P1 = 0.066197076394253;
1020           const double P2 = 0.062969590272413;
1021           add(  1/3.,  1/3., 9/80. ); 
1022           add(     A,     A, P1 ); 
1023           add( 1-2*A,     A, P1 );
1024           add(     A, 1-2*A, P1 );
1025           add(     B,     B, P2 ); 
1026           add( 1-2*B,     B, P2 );
1027           add(     B, 1-2*B, P2 ); break;
1028         }
1029         case 12: { // FPG12
1030           const double A  = 0.063089014491502;
1031           const double B  = 0.249286745170910;
1032           const double C  = 0.310352451033785;
1033           const double D  = 0.053145049844816;
1034           const double P1 = 0.025422453185103;
1035           const double P2 = 0.058393137863189;
1036           const double P3 = 0.041425537809187;
1037           add(     A,     A, P1 ); 
1038           add( 1-2*A,     A, P1 );
1039           add(     A, 1-2*A, P1 );
1040           add(     B,     B, P2 ); 
1041           add( 1-2*B,     B, P2 );
1042           add(     B, 1-2*B, P2 );
1043           add(     C,     D, P3 );
1044           add(     D,     C, P3 );
1045           add( 1-C-D,     C, P3 );
1046           add( 1-C-D,     D, P3 );
1047           add(     C, 1-C-D, P3 );
1048           add(     D, 1-C-D, P3 ); break;
1049         }
1050         default:
1051           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 1: "
1052                      <<nbGauss);
1053         }
1054       }
1055       else if ( variant == 2 ) {
1056         if (geom == NORM_TRI3) setRefCoords( TTria3a() );
1057         else                   setRefCoords( TTria6a() );
1058         switch ( nbGauss ) {
1059         case 1: {
1060           add( -1/3., -1/3., 2. ); break;
1061         }
1062         case 3: {
1063           add( -2/3.,  1/3., 2/3. );
1064           add( -2/3., -2/3., 2/3. );
1065           add(  1/3., -2/3., 2/3. ); break;
1066         }
1067         case 6: {
1068           const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
1069           const double A  = 0.445948490915965, B = 0.091576213509771;
1070           add( 2*B-1, 1-4*B, 4*P2 ); 
1071           add( 2*B-1, 2*B-1, 4*P2 );
1072           add( 1-4*B, 2*B-1, 4*P2 );
1073           add( 1-4*A, 2*A-1, 4*P1 );
1074           add( 2*A-1, 1-4*A, 4*P1 ); 
1075           add( 2*A-1, 2*A-1, 4*P1 ); break;
1076         }
1077         default:
1078           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 2: "
1079                      <<nbGauss);
1080         }
1081       }
1082       else if ( variant == 3 ) {
1083         if (geom == NORM_TRI3) setRefCoords( TTria3b() );
1084         else                   setRefCoords( TTria6b() );
1085         switch ( nbGauss ) {
1086         case 4: {
1087           add( 1/3., 1/3., -27/96 );
1088           add( 0.2 , 0.2 ,  25/96 );
1089           add( 0.6 , 0.2 ,  25/96 );
1090           add( 0.2 , 0.6 ,  25/96 ); break;
1091         }
1092         default:
1093           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 3: "
1094                      <<nbGauss);
1095         }
1096       }
1097       break;
1098
1099     case NORM_QUAD4:
1100     case NORM_QUAD8:
1101       if ( variant == 1 ) {
1102         if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1103         else                    setRefCoords( TQuad8b() );
1104         switch ( nbGauss ) {
1105         case 1: { // FPG1
1106           add(  0,  0,  4 ); break;
1107         }
1108         case 4: { // FPG4
1109           const double a = 1/sqrt(3.);
1110           add( -a, -a,  1 );
1111           add(  a, -a,  1 );
1112           add(  a,  a,  1 );
1113           add( -a,  a,  1 ); break;
1114         }
1115         case 5: { // out from the 3 specs
1116           const double a = 0.774596669241483;
1117           add( -a, -a,  0.5 );
1118           add(  a, -a,  0.5 );
1119           add(  a,  a,  0.5 );
1120           add( -a,  a,  0.5 );
1121           add(  0,  0,  2.0 ); break;
1122         }
1123         case 9: { // FPG9
1124           const double a = 0.774596669241483;
1125           add( -a, -a,  25/81. );
1126           add(  a, -a,  25/81. );
1127           add(  a,  a,  25/81. );
1128           add( -a,  a,  25/81. );
1129           add( 0., -a,  40/81. );
1130           add(  a, 0.,  40/81. );
1131           add( 0.,  a,  40/81. );
1132           add( -a, 0.,  40/81. );
1133           add( 0., 0.,  64/81. ); break;
1134         }
1135         default:
1136           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1137                      <<nbGauss);
1138         }
1139       }
1140       else if ( variant == 2 ) {
1141         if (geom == NORM_QUAD4) setRefCoords( TQuad4a() );
1142         else                    setRefCoords( TQuad8a() );
1143         switch ( nbGauss ) {
1144         case 4: {
1145           const double a = 1/sqrt(3.);
1146           add( -a,  a,  1 );
1147           add( -a, -a,  1 );
1148           add(  a, -a,  1 );
1149           add(  a,  a,  1 ); break;
1150         }
1151         case 9: {
1152           const double a = 0.774596669241483;
1153           add( -a,  a,  25/81. );
1154           add( -a, -a,  25/81. );
1155           add(  a, -a,  25/81. );
1156           add(  a,  a,  25/81. );
1157           add( -a, 0.,  40/81. );
1158           add( 0., -a,  40/81. );
1159           add(  a, 0.,  40/81. );
1160           add( 0.,  a,  40/81. );
1161           add( 0., 0.,  64/81. ); break;
1162         }
1163         default:
1164           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1165                      <<nbGauss);
1166         }
1167       }
1168       else if ( variant == 3 ) {
1169         if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1170         else                    setRefCoords( TQuad8b() );
1171         switch ( nbGauss ) {
1172         case 4: {
1173           const double a = 3/sqrt(3.);
1174           add( -a, -a,  1 );
1175           add( -a,  a,  1 );
1176           add(  a, -a,  1 );
1177           add(  a,  a,  1 ); break;
1178         }
1179         case 9: {
1180           const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
1181           const double c12 = c1*c2, c22 = c2*c2, c1c2 = c1*c2;
1182           add( -a, -a,  c12  );
1183           add( -a, 0.,  c1c2 );
1184           add( -a,  a,  c12  );
1185           add( 0., -a,  c1c2 );
1186           add( 0., 0.,  c22  );
1187           add( 0.,  a,  c1c2 );
1188           add(  a, -a,  c12  );
1189           add(  a, 0.,  c1c2 );
1190           add(  a,  a,  c12  ); break;
1191         }
1192         default:
1193           THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 3: "
1194                      <<nbGauss);
1195         }
1196       }
1197       break;
1198
1199     case NORM_TETRA4:
1200     case NORM_TETRA10:
1201       if (geom == NORM_TETRA4) setRefCoords( TTetra4a() );
1202       else                 setRefCoords( TTetra10a() );
1203       switch ( nbGauss ) {
1204       case 4: { // FPG4
1205         const double a = (5 - sqrt(5.))/20., b = (5 + 3*sqrt(5.))/20.;
1206         add(  a,  a,  a,  1/24. );
1207         add(  a,  a,  b,  1/24. );
1208         add(  a,  b,  a,  1/24. );
1209         add(  b,  a,  a,  1/24. ); break;
1210       }
1211       case 5: { // FPG5
1212         const double a = 0.25, b = 1/6., c = 0.5;
1213         add(  a,  a,  a, -2/15. );
1214         add(  b,  b,  b,  3/40. );
1215         add(  b,  b,  c,  3/40. );
1216         add(  b,  c,  b,  3/40. );
1217         add(  c,  b,  b,  3/40. ); break;
1218       }
1219       case 15: { // FPG15
1220         const double a = 0.25;
1221         const double b1 = (7 + sqrt(15.))/34., c1 = (13 + 3*sqrt(15.))/34., d = (5 - sqrt(15.))/20.;
1222         const double b2 = (7 - sqrt(15.))/34., c2 = (13 - 3*sqrt(15.))/34., e = (5 + sqrt(15.))/20.;
1223         const double P1 = (2665 - 14*sqrt(15.))/226800.;
1224         const double P2 = (2665 + 14*sqrt(15.))/226800.;
1225         add(  a,  a,  a,  8/405.);//_____
1226         add( b1, b1, b1,  P1    );
1227         add( b1, b1, c1,  P1    );
1228         add( b1, c1, b1,  P1    );
1229         add( c1, b1, b1,  P1    );//_____
1230         add( b2, b2, b2,  P2    );
1231         add( b2, b2, c2,  P2    );
1232         add( b2, c2, b2,  P2    );
1233         add( c2, b2, b2,  P2    );//_____
1234         add(  d,  d,  e,  5/567.);
1235         add(  d,  e,  d,  5/567.);
1236         add(  e,  d,  d,  5/567.);
1237         add(  d,  e,  e,  5/567.);
1238         add(  e,  d,  e,  5/567.);
1239         add(  e,  e,  d,  5/567.);
1240         break;
1241       }
1242       default:
1243         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TETRA: "<<nbGauss);
1244       }
1245       break;
1246
1247     case NORM_PYRA5:
1248     case NORM_PYRA13:
1249       if (geom == NORM_PYRA5) setRefCoords( TPyra5a() );
1250       else                setRefCoords( TPyra13a() );
1251       switch ( nbGauss ) {
1252       case 5: { // FPG5
1253         const double h1 = 0.1531754163448146;
1254         const double h2 = 0.6372983346207416;
1255         add(  .5,  0.,  h1,  2/15. );
1256         add(  0.,  .5,  h1,  2/15. );
1257         add( -.5,  0.,  h1,  2/15. );
1258         add(  0., -.5,  h1,  2/15. );
1259         add(  0.,  0.,  h2,  2/15. ); break;
1260       }
1261       case 6: { // FPG6
1262         const double p1 = 0.1024890634400000 ;
1263         const double p2 = 0.1100000000000000 ;
1264         const double p3 = 0.1467104129066667 ;
1265         const double a  = 0.5702963741068025 ;
1266         const double h1 = 0.1666666666666666 ;
1267         const double h2 = 0.08063183038464675;
1268         const double h3 = 0.6098484849057127 ;
1269         add(  a, 0.,  h1,  p1 );
1270         add( 0.,  a,  h1,  p1 );
1271         add( -a, 0.,  h1,  p1 );
1272         add( 0., -a,  h1,  p1 );
1273         add( 0., 0.,  h2,  p2 );
1274         add( 0., 0.,  h3,  p3 ); break;
1275       }
1276       case 27: { // FPG27
1277         const double a1  = 0.788073483; 
1278         const double b6  = 0.499369002; 
1279         const double b1  = 0.848418011; 
1280         const double c8  = 0.478508449; 
1281         const double c1  = 0.652816472; 
1282         const double d12 = 0.032303742; 
1283         const double d1  = 1.106412899;
1284         double z = 1/2., fz = b1/2*(1 - z);
1285         add(  0.,  0.,   z,  a1 ); // 1
1286         add(  fz,  fz,   z,  b6 ); // 2
1287         add( -fz,  fz,   z,  b6 ); // 3
1288         add( -fz, -fz,   z,  b6 ); // 4
1289         add(  fz, -fz,   z,  b6 ); // 5
1290         z = (1 - b1)/2.;
1291         add(  0.,  0.,   z,  b6 ); // 6
1292         z = (1 + b1)/2.;
1293         add(  0.,  0.,   z,  b6 ); // 7
1294         z = (1 - c1)/2.; fz = c1*(1 - z);
1295         add(  fz,  0.,   z,  c8 ); // 8
1296         add(  0.,  fz,   z,  c8 ); // 9
1297         add( -fz,  0.,   z,  c8 ); // 10
1298         add(  0., -fz,   z,  c8 ); // 11
1299         z = (1 + c1)/2.; fz = c1*(1 - z);
1300         add(  fz,  0.,   z,  c8 ); // 12
1301         add(  0.,  fz,   z,  c8 ); // 13
1302         add( -fz,  0.,   z,  c8 ); // 14
1303         add(  0., -fz,   z,  c8 ); // 15
1304         z = (1 - d1)/2., fz = d1/2*(1 - z);
1305         add(  fz,  fz,   z,  d12); // 16
1306         add( -fz,  fz,   z,  d12); // 17
1307         add( -fz, -fz,   z,  d12); // 18
1308         add(  fz, -fz,   z,  d12); // 19
1309         z = 1/2.; fz = d1*(1 - z);
1310         add(  fz,  0.,   z,  d12); // 20
1311         add(  0.,  fz,   z,  d12); // 21
1312         add( -fz,  0.,   z,  d12); // 22
1313         add(  0., -fz,   z,  d12); // 23
1314         z = (1 + d1)/2., fz = d1/2*(1 - z);
1315         add(  fz,  fz,   z,  d12); // 24
1316         add( -fz,  fz,   z,  d12); // 25
1317         add( -fz, -fz,   z,  d12); // 26
1318         add(  fz, -fz,   z,  d12); // 27
1319         break;
1320       }
1321       default:
1322         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PYRA: "<<nbGauss);
1323       }
1324       break;
1325     case NORM_PENTA6:
1326     case NORM_PENTA15:
1327       if (geom == NORM_PENTA6) setRefCoords( TPenta6a() );
1328       else                 setRefCoords( TPenta15a() );
1329       switch ( nbGauss ) {
1330       case 6: { // FPG6
1331         const double a = sqrt(3.)/3.;
1332         add( -a, .5, .5,  1/6. );
1333         add( -a, 0., .5,  1/6. );
1334         add( -a, .5, 0.,  1/6. );
1335         add(  a, .5, .5,  1/6. );
1336         add(  a, 0., .5,  1/6. );
1337         add(  a, .5, 0.,  1/6. ); break;
1338       }
1339       case 8: { // FPG8
1340         const double a = 0.577350269189626;
1341         add( -a, 1/3., 1/3., -27/96. );
1342         add( -a,  0.6,  0.2,  25/96. );
1343         add( -a,  0.2,  0.6,  25/96. );
1344         add( -a,  0.2,  0.2,  25/96. );
1345         add( +a, 1/3., 1/3., -27/96. );
1346         add( +a,  0.6,  0.2,  25/96. );
1347         add( +a,  0.2,  0.6,  25/96. );
1348         add( +a,  0.2,  0.2,  25/96. ); break;
1349       }
1350       case 21: { // FPG21
1351         const double d = sqrt(3/5.), c1 = 5/9., c2 = 8/9.; // d <=> alfa
1352         const double a = (6 + sqrt(15.))/21.;
1353         const double b = (6 - sqrt(15.))/21.;
1354         const double P1 = (155 + sqrt(15.))/2400.;
1355         const double P2 = (155 - sqrt(15.))/2400.;  //___
1356         add( -d,  1/3.,  1/3., c1*9/80. );//___
1357         add( -d,     a,     a, c1*P1    );
1358         add( -d, 1-2*a,     a, c1*P1    );
1359         add( -d,     a, 1-2*a, c1*P1    );//___
1360         add( -d,     b,     b, c1*P2    );
1361         add( -d, 1-2*b,     b, c1*P2    );
1362         add( -d,     b, 1-2*b, c1*P2    );//___
1363         add( 0.,  1/3.,  1/3., c2*9/80. );//___
1364         add( 0.,     a,     a, c2*P1    );
1365         add( 0., 1-2*a,     a, c2*P1    );
1366         add( 0.,     a, 1-2*a, c2*P1    );//___
1367         add( 0.,     b,     b, c2*P2    );
1368         add( 0., 1-2*b,     b, c2*P2    );
1369         add( 0.,     b, 1-2*b, c2*P2    );//___
1370         add(  d,  1/3.,  1/3., c1*9/80. );//___
1371         add(  d,     a,     a, c1*P1    );
1372         add(  d, 1-2*a,     a, c1*P1    );
1373         add(  d,     a, 1-2*a, c1*P1    );//___
1374         add(  d,     b,     b, c1*P2    );
1375         add(  d, 1-2*b,     b, c1*P2    );
1376         add(  d,     b, 1-2*b, c1*P2    );//___
1377         break;
1378       }
1379       default:
1380         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1381       }
1382       break;
1383
1384     case NORM_HEXA8:
1385     case NORM_HEXA20:
1386       if (geom == NORM_HEXA8) setRefCoords( THexa8a() );
1387       else                    setRefCoords( THexa20a() );
1388       switch ( nbGauss ) {
1389       case 8: { // FPG8
1390         const double a = sqrt(3.)/3.;
1391         add( -a, -a, -a,  1. );
1392         add( -a, -a,  a,  1. );
1393         add( -a,  a, -a,  1. );
1394         add( -a,  a,  a,  1. );
1395         add(  a, -a, -a,  1. );
1396         add(  a, -a,  a,  1. );
1397         add(  a,  a, -a,  1. );
1398         add(  a,  a,  a,  1. ); break;
1399       }
1400       case 27: { // FPG27
1401         const double a = sqrt(3/5.), c1 = 5/9., c2 = 8/9.;
1402         const double c12 = c1*c1, c13 = c1*c1*c1;
1403         const double c22 = c2*c2, c23 = c2*c2*c2;
1404         add( -a, -a, -a,   c13  ); // 1
1405         add( -a, -a, 0., c12*c2 ); // 2
1406         add( -a, -a,  a,   c13  ); // 3
1407         add( -a, 0., -a, c12*c2 ); // 4
1408         add( -a, 0., 0., c1*c22 ); // 5
1409         add( -a, 0.,  a, c12*c2 ); // 6
1410         add( -a,  a, -a,   c13  ); // 7
1411         add( -a,  a, 0., c12*c2 ); // 8
1412         add( -a,  a,  a,   c13  ); // 9
1413         add( 0., -a, -a, c12*c2 ); // 10
1414         add( 0., -a, 0., c1*c22 ); // 11
1415         add( 0., -a,  a, c12*c2 ); // 12
1416         add( 0., 0., -a, c1*c22 ); // 13
1417         add( 0., 0., 0.,   c23  ); // 14
1418         add( 0., 0.,  a, c1*c22 ); // 15
1419         add( 0.,  a, -a, c12*c2 ); // 16
1420         add( 0.,  a, 0., c1*c22 ); // 17
1421         add( 0.,  a,  a, c12*c2 ); // 18
1422         add(  a, -a, -a,   c13  ); // 19
1423         add(  a, -a, 0., c12*c2 ); // 20
1424         add(  a, -a,  a,   c13  ); // 21
1425         add(  a, 0., -a, c12*c2 ); // 22
1426         add(  a, 0., 0., c1*c22 ); // 23
1427         add(  a, 0.,  a, c12*c2 ); // 24
1428         add(  a,  a, -a,   c13  ); // 25
1429         add(  a,  a, 0., c12*c2 ); // 26
1430         add(  a,  a,  a,   c13  ); // 27
1431         break;
1432       }
1433       default:
1434         THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1435       }
1436       break;
1437
1438     default:
1439       THROW_IK_EXCEPTION("TGaussDef: unexpected EGeometrieElement: "<< geom);
1440     }
1441
1442     if ( myWeights.capacity() != myWeights.size() )
1443       THROW_IK_EXCEPTION("TGaussDef: Not all gauss points defined");
1444   }
1445 }
1446   
1447 //================================================================================
1448 /*!
1449  * \brief Return dimension for the given cell type
1450  */
1451 //================================================================================
1452
1453 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
1454 {
1455   return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
1456 }
1457
1458 //================================================================================
1459 /*!
1460  * \brief Returns interlace array to transform a quadratic GIBI element to a MED one.
1461  *        i-th array item gives node index in GIBI connectivity for i-th MED node
1462  */
1463 //================================================================================
1464
1465 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
1466 {
1467   static vector<const int*> conn;
1468   static const int hexa20 [] = {0,6,4,2, 12,18,16,14, 7,5,3,1, 19,17,15,13, 8,11,10,9};
1469   static const int penta15[] = {0,2,4, 9,11,13, 1,3,5, 10,12,14, 6,8,7};
1470   static const int pyra13 [] = {0,2,4,6, 12, 1,3,5,7, 8,9,10,11};
1471   static const int tetra10[] = {0,2,4, 9, 1,3,5, 6,7,8};
1472   static const int quad8  [] = {0,2,4,6, 1,3,5,7};
1473   static const int tria6  [] = {0,2,4, 1,3,5};
1474   static const int seg3   [] = {0,2,1};
1475   if ( conn.empty() )
1476   {
1477     conn.resize( MaxMedCellType + 1, 0 );
1478     conn[ NORM_HEXA20 ] = hexa20;
1479     conn[ NORM_PENTA15] = penta15;
1480     conn[ NORM_PYRA13 ] = pyra13;
1481     conn[ NORM_TETRA10] = tetra10;
1482     conn[ NORM_SEG3   ] = seg3;
1483     conn[ NORM_TRI6   ] = tria6;
1484     conn[ NORM_QUAD8  ] = quad8;
1485   }
1486   return conn[ type ];
1487 }
1488
1489 //================================================================================
1490 /*!
1491  * \brief Avoid coping sortedNodeIDs
1492  */
1493 //================================================================================
1494
1495 Cell::Cell(const Cell& ma)
1496   : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
1497 {
1498   if ( ma._sortedNodeIDs )
1499     {
1500       _sortedNodeIDs = new int[ _nodes.size() ];
1501       std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
1502     }
1503 }
1504
1505 //================================================================================
1506 /*!
1507  * \brief Rerturn the i-th link of face
1508  */
1509 //================================================================================
1510
1511 SauvUtilities::Link Cell::link(int i) const
1512 {
1513   int i2 = ( i + 1 ) % _nodes.size();
1514   if ( _reverse )
1515     return make_pair( _nodes[i2]->_number, _nodes[i]->_number );
1516   else
1517     return make_pair( _nodes[i]->_number, _nodes[i2]->_number );
1518 }
1519
1520 //================================================================================
1521 /*!
1522  * \brief creates if needed and return _sortedNodeIDs
1523  */
1524 //================================================================================
1525
1526 const TID* Cell::getSortedNodes() const
1527 {
1528   if ( !_sortedNodeIDs )
1529   {
1530     size_t l=_nodes.size();
1531     _sortedNodeIDs = new int[ l ];
1532
1533     for (size_t i=0; i!=l; ++i)
1534       _sortedNodeIDs[i]=_nodes[i]->_number;
1535     std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
1536   }
1537   return _sortedNodeIDs;
1538 }
1539
1540 //================================================================================
1541 /*!
1542  * \brief Compare sorted ids of cell nodes
1543  */
1544 //================================================================================
1545
1546 bool Cell::operator< (const Cell& ma) const
1547 {
1548   if ( _nodes.size() == 1 )
1549     return _nodes[0] < ma._nodes[0];
1550
1551   const int* v1 = getSortedNodes();
1552   const int* v2 = ma.getSortedNodes();
1553   for ( const int* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
1554     if(*v1 != *v2)
1555       return *v1 < *v2;
1556   return false;
1557 }
1558
1559 //================================================================================
1560 /*!
1561  * \brief Dump a Cell
1562  */
1563 //================================================================================
1564
1565 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
1566 {
1567   os << "cell " << ma._number << " (" << ma._nodes.size() << " nodes) : < " << ma._nodes[0]->_number;
1568   for( size_t i=1; i!=ma._nodes.size(); ++i)
1569     os << ", " << ma._nodes[i]->_number;
1570 #ifdef _DEBUG_
1571   os << " > sortedNodes: ";
1572   if ( ma._sortedNodeIDs ) {
1573     os << "< ";
1574     for( size_t i=0; i!=ma._nodes.size(); ++i)
1575       os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
1576     os << " >";
1577   }
1578   else {
1579     os << "NULL";
1580   }
1581 #endif
1582   return os;
1583 }
1584
1585 //================================================================================
1586 /*!
1587  * \brief Return nb of elements in the group
1588  */
1589 //================================================================================
1590
1591 int Group::size() const
1592 {
1593   int sizze = 0;
1594   if ( !_relocTable.empty() )
1595     sizze =  _relocTable.size();
1596   else if ( _medGroup )
1597     sizze = _medGroup->getNumberOfTuples();
1598   else if ( !_cells.empty() )
1599     sizze = _cells.size();
1600   else
1601     for ( size_t i = 0; i < _groups.size(); ++i )
1602       sizze += _groups[i]->size();
1603   return sizze;
1604 }
1605
1606 //================================================================================
1607 /*!
1608  * \brief Conver gibi element type to med one
1609  */
1610 //================================================================================
1611
1612 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
1613 {
1614   if ( gibiType < 1 || gibiType > NbGibiCellTypes )
1615     return NORM_ERROR;
1616
1617   return GibiTypeToMed[ gibiType - 1 ];
1618 }
1619
1620 //================================================================================
1621 /*!
1622  * \brief Conver med element type to gibi one
1623  */
1624 //================================================================================
1625
1626 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
1627 {
1628   for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
1629     if ( GibiTypeToMed[ i ] == medGeomType )
1630       return i + 1;
1631
1632   return -1;
1633 }
1634
1635 //================================================================================
1636 /*!
1637  * \brief Remember the file name
1638  */
1639 //================================================================================
1640
1641 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
1642 {
1643 }
1644
1645 //================================================================================
1646 /*!
1647  * \brief Constructor of ASCII sauve file reader
1648  */
1649 //================================================================================
1650
1651 ASCIIReader::ASCIIReader(const char* fileName)
1652   :FileReader(fileName),
1653    _file(-1)
1654 {
1655 }
1656
1657 //================================================================================
1658 /*!
1659  * \brief Return true
1660  */
1661 //================================================================================
1662
1663 bool ASCIIReader::isASCII() const
1664 {
1665   return true;
1666 }
1667
1668 //================================================================================
1669 /*!
1670  * \brief Try to open an ASCII file
1671  */
1672 //================================================================================
1673
1674 bool ASCIIReader::open()
1675 {
1676 #ifdef WIN32
1677   _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
1678 #else
1679   _file = ::open (_fileName.c_str(), O_RDONLY);
1680 #endif
1681   if (_file >= 0)
1682     {
1683       _start  = new char [GIBI_BufferSize]; // working buffer beginning
1684       //_tmpBuf = new char [GIBI_MaxOutputLen];
1685       _ptr    = _start;
1686       _eptr   = _start;
1687       _lineNb = 0;
1688     }
1689   else
1690     {
1691       //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
1692     }
1693   return (_file >= 0);
1694 }
1695
1696 //================================================================================
1697 /*!
1698  * \brief Close the file
1699  */
1700 //================================================================================
1701
1702 ASCIIReader::~ASCIIReader()
1703 {
1704   if (_file >= 0)
1705     {
1706       ::close (_file);
1707       if (_start != 0L)
1708         {
1709           delete [] _start;
1710           //delete [] _tmpBuf;
1711           _start = 0;
1712         }
1713       _file = -1;
1714     }
1715 }
1716
1717 //================================================================================
1718 /*!
1719  * \brief Return a next line of the file
1720  */
1721 //================================================================================
1722
1723 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
1724 {
1725   if ( getLine( line )) return true;
1726   if ( raiseOEF )
1727     THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
1728   return false;
1729 }
1730
1731 //================================================================================
1732 /*!
1733  * \brief Read a next line of the file if necessary
1734  */
1735 //================================================================================
1736
1737 bool ASCIIReader::getLine(char* & line)
1738 {
1739   bool aResult = true;
1740   // Check the state of the buffer;
1741   // if there is too little left, read the next portion of data
1742   int nBytesRest = _eptr - _ptr;
1743   if (nBytesRest < GIBI_MaxOutputLen)
1744     {
1745       if (nBytesRest > 0)
1746         {
1747           // move the remaining portion to the buffer beginning
1748           for ( int i = 0; i < nBytesRest; ++i )
1749             _start[i] = _ptr[i];
1750           //memcpy (_tmpBuf, _ptr, nBytesRest);
1751           //memcpy (_start, _tmpBuf, nBytesRest);
1752         }
1753       else
1754         {
1755           nBytesRest = 0;
1756         }
1757       _ptr = _start;
1758       const int nBytesRead = ::read (_file,
1759                                      &_start [nBytesRest],
1760                                      GIBI_BufferSize - nBytesRest);
1761       nBytesRest += nBytesRead;
1762       _eptr = &_start [nBytesRest];
1763     }
1764   // Check the buffer for the end-of-line
1765   char * ptr = _ptr;
1766   while (true)
1767     {
1768       // Check for end-of-the-buffer, the ultimate criterion for termination
1769       if (ptr >= _eptr)
1770         {
1771           if (nBytesRest <= 0)
1772             aResult = false;
1773           else
1774             _eptr[-1] = '\0';
1775           break;
1776         }
1777       // seek the line-feed character
1778       if (ptr[0] == '\n')
1779         {
1780           if (ptr[-1] == '\r')
1781             ptr[-1] = '\0';
1782           ptr[0] = '\0';
1783           ++ptr;
1784           break;
1785         }
1786       ++ptr;
1787     }
1788   // Output the result
1789   line = _ptr;
1790   _ptr = ptr;
1791   _lineNb++;
1792
1793   return aResult;
1794 }
1795
1796 //================================================================================
1797 /*!
1798  * \brief Prepare for iterating over given nb of values
1799  *  \param nbToRead - nb of fields to read
1800  *  \param nbPosInLine - nb of fields in one line
1801  *  \param width - field width
1802  *  \param shift - shift from line beginning to the field start
1803  */
1804 //================================================================================
1805
1806 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
1807 {
1808   _nbToRead    = nbToRead;
1809   _nbPosInLine = nbPosInLine;
1810   _width       = width;
1811   _shift       = shift;
1812   _iPos = _iRead = 0;
1813   if ( _nbToRead )
1814     {
1815       getNextLine( _curPos );
1816       _curPos = _curPos + _shift;
1817     }
1818   else
1819     {
1820       _curPos = 0;
1821     }
1822   _curLocale.clear();
1823 }
1824
1825 //================================================================================
1826 /*!
1827  * \brief Prepare for iterating over given nb of string values
1828  */
1829 //================================================================================
1830
1831 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
1832 {
1833   init( nbValues, 72 / ( width + 1 ), width, 1 );
1834 }
1835
1836 //================================================================================
1837 /*!
1838  * \brief Prepare for iterating over given nb of integer values
1839  */
1840 //================================================================================
1841
1842 void ASCIIReader::initIntReading(int nbValues)
1843 {
1844   init( nbValues, 10, 8 );
1845 }
1846
1847 //================================================================================
1848 /*!
1849  * \brief Prepare for iterating over given nb of real values
1850  */
1851 //================================================================================
1852
1853 void ASCIIReader::initDoubleReading(int nbValues)
1854 {
1855   init( nbValues, 3, 22 );
1856
1857   // Correction 2 of getDouble(): set "C" numeric locale to read numbers
1858   // with dot decimal point separator, as it is in SAUVE files
1859   _curLocale = setlocale(LC_NUMERIC, "C");
1860 }
1861
1862 //================================================================================
1863 /*!
1864  * \brief Return true if not all values have been read
1865  */
1866 //================================================================================
1867
1868 bool ASCIIReader::more() const
1869 {
1870   bool result = false;
1871   if ( _iRead < _nbToRead)
1872     {
1873       if ( _curPos ) result = true;
1874     }
1875   return result;
1876 }
1877
1878 //================================================================================
1879 /*!
1880  * \brief Go to the nex value
1881  */
1882 //================================================================================
1883
1884 void ASCIIReader::next()
1885 {
1886   if ( !more() )
1887     THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
1888   ++_iRead;
1889   ++_iPos;
1890   if ( _iRead < _nbToRead )
1891     {
1892       if ( _iPos >= _nbPosInLine )
1893         {
1894           getNextLine( _curPos );
1895           _curPos = _curPos + _shift;
1896           _iPos = 0;
1897         }
1898       else
1899         {
1900           _curPos = _curPos + _width + _shift;
1901         }
1902     }
1903   else
1904     {
1905       _curPos = 0;
1906       if ( !_curLocale.empty() )
1907         {
1908           setlocale(LC_NUMERIC, _curLocale.c_str());
1909           _curLocale.clear();
1910         }
1911     }
1912 }
1913
1914 //================================================================================
1915 /*!
1916  * \brief Return the current integer value
1917  */
1918 //================================================================================
1919
1920 int ASCIIReader::getInt() const
1921 {
1922   // fix for two glued ints (issue 0021009):
1923   // Line nb    |   File contents
1924   // ------------------------------------------------------------------------------------
1925   // 53619905   |       1       2       6       8
1926   // 53619906   |                                                                SCALAIRE
1927   // 53619907   |    -63312600499       1       0       0       0      -2       0       2
1928   //   where -63312600499 is actualy -633 and 12600499
1929   char hold=_curPos[_width];
1930   _curPos[_width] = '\0';
1931   int result = atoi( _curPos );
1932   _curPos[_width] = hold;
1933   return result;
1934   //return atoi(str());
1935 }
1936
1937 //================================================================================
1938 /*!
1939  * \brief Return the current float value
1940  */
1941 //================================================================================
1942
1943 float ASCIIReader::getFloat() const
1944 {
1945   return getDouble();
1946 }
1947
1948 //================================================================================
1949 /*!
1950  * \brief Return the current double value
1951  */
1952 //================================================================================
1953
1954 double ASCIIReader::getDouble() const
1955 {
1956   //std::string aStr (_curPos);
1957
1958   // Correction: add missing 'E' specifier
1959   // int aPosStart = aStr.find_first_not_of(" \t");
1960   // if (aPosStart < (int)aStr.length()) {
1961   //   int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
1962   //   if (aPosSign < (int)aStr.length()) {
1963   //     if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
1964   //       aStr.insert(aPosSign, "E", 1);
1965   //   }
1966   // }
1967
1968   // Different Correction (more optimal)
1969   // Sample:
1970   //  0.00000000000000E+00 -2.37822406690632E+01  6.03062748797469E+01
1971   //  7.70000000000000-100  7.70000000000000+100  7.70000000000000+100
1972   //0123456789012345678901234567890123456789012345678901234567890123456789
1973   const size_t posE = 18;
1974   std::string aStr (_curPos);
1975   if ( aStr.find('E') < 0 && aStr.find('e') < 0 )
1976     {
1977       if ( aStr.size() < posE+1 )
1978         THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
1979       aStr.insert( posE, "E", 1 );
1980       return atof(aStr.c_str());
1981     }
1982   return atof( _curPos );
1983 }
1984
1985 //================================================================================
1986 /*!
1987  * \brief Return the current string value
1988  */
1989 //================================================================================
1990
1991 string ASCIIReader::getName() const
1992 {
1993   int len = _width;
1994   while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
1995     len--;
1996   return string( _curPos, len );
1997 }
1998
1999 //================================================================================
2000 /*!
2001  * \brief Constructor of a binary sauve file reader
2002  */
2003 //================================================================================
2004
2005 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
2006 {
2007 }
2008
2009 //================================================================================
2010 /*!
2011  * \brief Close the XDR sauve file
2012  */
2013 //================================================================================
2014
2015 XDRReader::~XDRReader()
2016 {
2017 #ifdef HAS_XDR  
2018   if ( _xdrs_file )
2019     {
2020       xdr_destroy((XDR*)_xdrs);
2021       free((XDR*)_xdrs);
2022       ::fclose(_xdrs_file);
2023       _xdrs_file = NULL;
2024     }
2025 #endif
2026 }
2027
2028 //================================================================================
2029 /*!
2030  * \brief Return false
2031  */
2032 //================================================================================
2033
2034 bool XDRReader::isASCII() const
2035 {
2036   return false;
2037 }
2038
2039 //================================================================================
2040 /*!
2041  * \brief Try to open an XRD file
2042  */
2043 //================================================================================
2044
2045 bool XDRReader::open()
2046 {
2047   bool xdr_ok = false;
2048 #ifdef HAS_XDR
2049   if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
2050     {
2051       _xdrs = (XDR *)malloc(sizeof(XDR));
2052       xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
2053
2054       const int maxsize = 10;
2055       char icha[maxsize+1];
2056       char* icha2 = icha;
2057       if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
2058         {
2059           icha[maxsize] = '\0';
2060           xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
2061         }
2062       if ( !xdr_ok )
2063         {
2064           xdr_destroy((XDR*)_xdrs);
2065           free((XDR*)_xdrs);
2066           fclose(_xdrs_file);
2067           _xdrs_file = NULL;
2068         }
2069     }
2070 #endif
2071   return xdr_ok;
2072 }
2073
2074 //================================================================================
2075 /*!
2076  * \brief A stub
2077  */
2078 //================================================================================
2079
2080 bool XDRReader::getNextLine (char* &, bool )
2081 {
2082   return true;
2083 }
2084
2085 //================================================================================
2086 /*!
2087  * \brief Prepare for iterating over given nb of values
2088  *  \param nbToRead - nb of fields to read
2089  *  \param width - field width
2090  */
2091 //================================================================================
2092
2093 void XDRReader::init( int nbToRead, int width/*=0*/ )
2094 {
2095   if(_iRead < _nbToRead)
2096     {
2097       cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << endl;
2098       cout << "Unfinished iteration before new one !" << endl;
2099       THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
2100     }
2101   _iRead    = 0;
2102   _nbToRead = nbToRead;
2103   _width    = width;
2104 }
2105
2106 //================================================================================
2107 /*!
2108  * \brief Prepare for iterating over given nb of string values
2109  */
2110 //================================================================================
2111
2112 void XDRReader::initNameReading(int nbValues, int width)
2113 {
2114   init( nbValues, width );
2115   _xdr_kind = _xdr_kind_char;
2116   if(nbValues*width)
2117     {
2118       unsigned int nels = nbValues*width;
2119       _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
2120 #ifdef HAS_XDR
2121       xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
2122 #endif
2123       _xdr_cvals[nels] = '\0';
2124     }
2125 }
2126
2127 //================================================================================
2128 /*!
2129  * \brief Prepare for iterating over given nb of integer values
2130  */
2131 //================================================================================
2132
2133 void XDRReader::initIntReading(int nbValues)
2134 {
2135   init( nbValues );
2136   _xdr_kind = _xdr_kind_int;
2137   if(nbValues)
2138     {
2139 #ifdef HAS_XDR
2140       unsigned int nels = nbValues;
2141       unsigned int actual_nels;
2142       _xdr_ivals = (int*)malloc(nels*sizeof(int));
2143       xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
2144 #endif
2145     }
2146 }
2147
2148 //================================================================================
2149 /*!
2150  * \brief Prepare for iterating over given nb of real values
2151  */
2152 //================================================================================
2153
2154 void XDRReader::initDoubleReading(int nbValues)
2155 {
2156   init( nbValues );
2157   _xdr_kind = _xdr_kind_double;
2158   if(nbValues)
2159     {
2160 #ifdef HAS_XDR
2161       unsigned int nels = nbValues;
2162       unsigned int actual_nels;
2163       _xdr_dvals = (double*)malloc(nels*sizeof(double));
2164       xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
2165 #endif
2166     }
2167 }
2168
2169 //================================================================================
2170 /*!
2171  * \brief Return true if not all values have been read
2172  */
2173 //================================================================================
2174
2175 bool XDRReader::more() const
2176 {
2177   return _iRead < _nbToRead;
2178 }
2179
2180 //================================================================================
2181 /*!
2182  * \brief Go to the nex value
2183  */
2184 //================================================================================
2185
2186 void XDRReader::next()
2187 {
2188   if ( !more() )
2189     THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
2190
2191   ++_iRead;
2192   if ( _iRead < _nbToRead )
2193     {
2194     }
2195   else
2196     {
2197       if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
2198       if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
2199       if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
2200       _xdr_kind = _xdr_kind_null;
2201     }
2202 }
2203
2204 //================================================================================
2205 /*!
2206  * \brief Return the current integer value
2207  */
2208 //================================================================================
2209
2210 int XDRReader::getInt() const
2211 {
2212   if(_iRead < _nbToRead)
2213     {
2214       return _xdr_ivals[_iRead];
2215     }
2216   else
2217     {
2218       int result = 0;
2219 #ifdef HAS_XDR
2220       xdr_int((XDR*)_xdrs, &result);
2221 #endif
2222       return result;
2223     }
2224 }
2225
2226 //================================================================================
2227 /*!
2228  * \brief Return the current float value
2229  */
2230 //================================================================================
2231
2232 float  XDRReader::getFloat() const
2233 {
2234   float result = 0;
2235 #ifdef HAS_XDR
2236   xdr_float((XDR*)_xdrs, &result);
2237 #endif
2238   return result;
2239 }
2240
2241 //================================================================================
2242 /*!
2243  * \brief Return the current double value
2244  */
2245 //================================================================================
2246
2247 double XDRReader::getDouble() const
2248 {
2249   if(_iRead < _nbToRead)
2250     {
2251       return _xdr_dvals[_iRead];
2252     }
2253   else
2254     {
2255       double result = 0;
2256 #ifdef HAS_XDR
2257       xdr_double((XDR*)_xdrs, &result);
2258 #endif
2259       return result;
2260     }
2261 }
2262
2263 //================================================================================
2264 /*!
2265  * \brief Return the current string value
2266  */
2267 //================================================================================
2268
2269 std::string XDRReader::getName() const
2270 {
2271   int len = _width;
2272   char* s = _xdr_cvals + _iRead*_width;
2273   while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
2274     len--;
2275   return string( s, len );
2276 }
2277
2278 //================================================================================
2279 /*!
2280  * \brief Throw an exception if not all needed data is present
2281  */
2282 //================================================================================
2283
2284 void IntermediateMED::checkDataAvailability() const
2285 {
2286   if ( _spaceDim == 0 )
2287     THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
2288
2289   if ( _groups.empty() )
2290     THROW_IK_EXCEPTION("No elements have been read");
2291
2292   if ( _points.empty() || _nbNodes == 0 )
2293     THROW_IK_EXCEPTION("Nodes of elements are not filled");
2294
2295   if ( _coords.empty() )
2296     THROW_IK_EXCEPTION("Node coordinates are missing");
2297
2298   if ( _coords.size() < _nbNodes * _spaceDim )
2299     THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
2300 }
2301
2302 //================================================================================
2303 /*!
2304  * \brief Makes ParaMEDMEM::MEDFileData from self
2305  */
2306 //================================================================================
2307
2308 ParaMEDMEM::MEDFileData* IntermediateMED::convertInMEDFileDS()
2309 {
2310   MEDCouplingAutoRefCountObjectPtr< MEDFileUMesh >  mesh   = makeMEDFileMesh();
2311   MEDCouplingAutoRefCountObjectPtr< MEDFileFields > fields = makeMEDFileFields(mesh);
2312
2313   MEDCouplingAutoRefCountObjectPtr< MEDFileMeshes > meshes = MEDFileMeshes::New();
2314   MEDCouplingAutoRefCountObjectPtr< MEDFileData >  medData = MEDFileData::New();
2315   meshes->pushMesh( mesh );
2316   medData->setMeshes( meshes );
2317   if ( fields ) medData->setFields( fields );
2318
2319   return medData.retn();
2320 }
2321
2322 //================================================================================
2323 /*!
2324  * \brief Creates ParaMEDMEM::MEDFileUMesh from its data
2325  */
2326 //================================================================================
2327
2328 ParaMEDMEM::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
2329 {
2330   // check if all needed piles are present
2331   checkDataAvailability();
2332
2333   // set long names
2334   setGroupLongNames();
2335
2336   // fix element orientation
2337   if ( _spaceDim == 2 || _spaceDim == 1 )
2338     orientElements2D();
2339   else if ( _spaceDim == 3 )
2340     orientElements3D();
2341
2342   // process groups
2343   decreaseHierarchicalDepthOfSubgroups();
2344   eraseUselessGroups();
2345   //detectMixDimGroups();
2346
2347   // assign IDs
2348   _points.numberNodes();
2349   numberElements();
2350
2351   // make the med mesh
2352
2353   MEDFileUMesh* mesh = MEDFileUMesh::New();
2354
2355   DataArrayDouble *coords = getCoords();
2356   setConnectivity( mesh, coords );
2357   setGroups( mesh );
2358
2359   coords->decrRef();
2360
2361   if ( !mesh->getName().c_str() || strlen( mesh->getName().c_str() ) == 0 )
2362     mesh->setName( "MESH" );
2363
2364   return mesh;
2365 }
2366
2367 //================================================================================
2368 /*!
2369  * \brief Set long names to groups
2370  */
2371 //================================================================================
2372
2373 void IntermediateMED::setGroupLongNames()
2374 {
2375   // IMP 0020434: mapping GIBI names to MED names
2376   // set med names to objects (mesh, fields, support, group or other)
2377
2378   set<int> treatedGroups;
2379
2380   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
2381   for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
2382     {
2383       if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
2384
2385       SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
2386
2387       // if there are several names for grp then the 1st name is the name
2388       // of grp and the rest ones are names of groups referring grp (issue 0021311)
2389       const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
2390       if ( !isRefName )
2391         {
2392           grp._name = _mapStrings[ itGIBItoMED->med_id ];
2393         }
2394       else if ( !grp._refNames.empty() && grp._refNames.back().empty() )
2395         {
2396           for ( unsigned i = 0; i < grp._refNames.size(); ++i )
2397             if ( grp._refNames[i].empty() )
2398               grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
2399         }
2400       else
2401         {
2402           grp._refNames.push_back( _mapStrings[ (*itGIBItoMED).med_id ]);
2403         }
2404     }
2405 }
2406
2407 //================================================================================
2408 /*!
2409  * \brief Set long names to fields
2410  */
2411 //================================================================================
2412
2413 void IntermediateMED::setFieldLongNames(set< string >& usedNames)
2414 {
2415   list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
2416   for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
2417     {
2418       if (itGIBItoMED->gibi_pile == PILE_FIELD)
2419         {
2420           _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2421         }
2422       else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
2423         {
2424           _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2425         }
2426     } // iterate on _listGIBItoMED_cham
2427
2428   for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
2429     {
2430       string medName  = _mapStrings[itGIBItoMED->med_id];
2431       string gibiName = _mapStrings[itGIBItoMED->gibi_id];
2432
2433       bool name_found = false;
2434       for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
2435         {
2436           vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
2437           for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
2438             {
2439               if (medName.find( fields[ifi]->_name + "." ) == 0 )
2440                 {
2441                   vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
2442                   int nbSub = aSubDs.size();
2443                   for (int isu = 0; isu < nbSub; isu++)
2444                     for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
2445                       {
2446                         if (aSubDs[isu].compName(ico) == gibiName)
2447                           {
2448                             string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
2449                             fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
2450                           }
2451                       }
2452                 }
2453             }
2454         }
2455     } // iterate on _listGIBItoMED_comp
2456
2457   for ( size_t i = 0; i < _nodeFields.size() ; i++)
2458     usedNames.insert( _nodeFields[i]->_name );
2459   for ( size_t i = 0; i < _cellFields.size() ; i++)
2460     usedNames.insert( _cellFields[i]->_name );
2461 }
2462
2463 //================================================================================
2464 /*!
2465  * \brief Decrease hierarchical depth of subgroups
2466  */
2467 //================================================================================
2468
2469 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
2470 {
2471   for (size_t i=0; i!=_groups.size(); ++i)
2472   {
2473     Group& grp = _groups[i];
2474     for (size_t j = 0; j < grp._groups.size(); ++j )
2475     {
2476       Group & sub_grp = *grp._groups[j];
2477       if ( !sub_grp._groups.empty() )
2478       {
2479         // replace j with its 1st subgroup
2480         grp._groups[j] = sub_grp._groups[0];
2481         // push back the rest subs
2482         grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
2483       }
2484     }
2485     // remove empty sub-_groups
2486     vector< Group* > newSubGroups;
2487     newSubGroups.reserve( grp._groups.size() );
2488     for (size_t j = 0; j < grp._groups.size(); ++j )
2489       if ( !grp._groups[j]->empty() )
2490         newSubGroups.push_back( grp._groups[j] );
2491     if ( newSubGroups.size() < grp._groups.size() )
2492       grp._groups.swap( newSubGroups );
2493   }
2494 }
2495
2496 //================================================================================
2497 /*!
2498  * \brief Erase _groups that won't be converted
2499  */
2500 //================================================================================
2501
2502 void IntermediateMED::eraseUselessGroups()
2503 {
2504   // propagate _isProfile=true to sub-groups of composite groups
2505   // for (size_t int i=0; i!=_groups.size(); ++i)
2506   // {
2507   //   Group* grp = _groups[i];
2508   //   if ( grp->_isProfile && !grp->_groups.empty() )
2509   //     for (size_t j = 0; j < grp->_groups.size(); ++j )
2510   //       grp->_groups[j]->_isProfile=true;
2511   // }
2512   std::set<Group*> groups2convert;
2513   // keep not named sub-groups of field supports
2514   for (size_t i=0; i!=_groups.size(); ++i)
2515   {
2516     Group& grp = _groups[i];
2517     if ( grp._isProfile && !grp._groups.empty() )
2518       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2519   }
2520
2521   // keep named groups and their subgroups
2522   for (size_t i=0; i!=_groups.size(); ++i)
2523   {
2524     Group& grp = _groups[i];
2525     if ( !grp._name.empty() && !grp.empty() )
2526     {
2527       groups2convert.insert( &grp );
2528       groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2529     }
2530   }
2531   // erase groups that are not in groups2convert and not _isProfile
2532   for (size_t i=0; i!=_groups.size(); ++i)
2533   {
2534     Group* grp = &_groups[i];
2535     if ( !grp->_isProfile && !groups2convert.count( grp ) )
2536     {
2537       grp->_cells.clear();
2538       grp->_groups.clear();
2539     }
2540   }
2541 }
2542
2543 //================================================================================
2544 /*!
2545  * \brief Detect _groups of mixed dimension
2546  */
2547 //================================================================================
2548
2549 void IntermediateMED::detectMixDimGroups()
2550 {
2551   //hasMixedCells = false;
2552   for ( size_t i=0; i < _groups.size(); ++i )
2553   {
2554     Group& grp = _groups[i];
2555     if ( grp._groups.size() < 2 )
2556       continue;
2557
2558     // check if sub-groups have different dimension
2559     unsigned dim1 = getDim( &grp );
2560     for ( size_t j = 1; j  < grp._groups.size(); ++j )
2561     {
2562       unsigned dim2 = getDim( grp._groups[j] );
2563       if ( dim1 != dim2 )
2564       {
2565         grp._cells.clear();
2566         grp._groups.clear();
2567         if ( !grp._name.empty() )
2568           cout << "Erase a group with elements of different dim |" << grp._name << "|"<< endl;
2569         break;
2570       }
2571     }
2572   }
2573 }
2574
2575 //================================================================================
2576 /*!
2577  * \brief Fix connectivity of elements in 2D space
2578  */
2579 //================================================================================
2580
2581 void IntermediateMED::orientElements2D()
2582 {
2583   set<Cell>::const_iterator elemIt, elemEnd;
2584   vector< pair<int,int> > swapVec;
2585
2586   // ------------------------------------
2587   // fix connectivity of quadratic edges
2588   // ------------------------------------
2589   set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
2590   if ( !quadEdges.empty() )
2591     {
2592       elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
2593       for ( ; elemIt != elemEnd; ++elemIt )
2594         ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
2595     }
2596
2597   CellsByDimIterator faceIt( *this, 2 );
2598   while ( const set<Cell > * faces = faceIt.nextType() )
2599     {
2600       TCellType cellType = faceIt.type();
2601       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2602
2603       getReverseVector( cellType, swapVec );
2604
2605       // ------------------------------------
2606       // fix connectivity of quadratic faces
2607       // ------------------------------------
2608       if ( isQuadratic )
2609         for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2610           ConvertQuadratic( cellType, *elemIt );
2611
2612       // --------------------------
2613       // orient faces clockwise
2614       // --------------------------
2615       int iQuad = isQuadratic ? 2 : 1;
2616       for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2617         {
2618           // look for index of the most left node
2619           int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
2620           double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
2621           for ( iNode = 1; iNode < nbNodes; ++iNode )
2622             if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
2623               minX = x, iLeft = iNode;
2624
2625           // indeces of the nodes neighboring the most left one
2626           int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
2627           int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
2628           // find components of prev-left and left-next vectors
2629           double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
2630           double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
2631           double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
2632           double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
2633           double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
2634           double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
2635           double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
2636           double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
2637           // normalise y of the vectors
2638           double modPL = sqrt ( xPL * xPL + yPL * yPL );
2639           double modLN = sqrt ( xLN * xLN + yLN * yLN );
2640           if ( modLN > std::numeric_limits<double>::min() &&
2641                modPL > std::numeric_limits<double>::min() )
2642             {
2643               yPL /= modPL;
2644               yLN /= modLN;
2645               // summary direction of neighboring links must be positive
2646               bool clockwise = ( yPL + yLN > 0 );
2647               if ( !clockwise )
2648                 reverse( *elemIt, swapVec );
2649             }
2650         }
2651     }
2652 }
2653
2654 //================================================================================
2655 /*!
2656  * \brief Fix connectivity of elements in 3D space
2657  */
2658 //================================================================================
2659
2660 void IntermediateMED::orientElements3D()
2661 {
2662   // set _reverse flags of faces
2663   orientFaces3D();
2664
2665   // -----------------
2666   // fix connectivity
2667   // -----------------
2668
2669   set<Cell>::const_iterator elemIt, elemEnd;
2670   vector< pair<int,int> > swapVec;
2671
2672   for ( int dim = 1; dim <= 3; ++dim )
2673   {
2674     CellsByDimIterator cellsIt( *this, dim );
2675     while ( const set<Cell > * elems = cellsIt.nextType() )
2676     {
2677       TCellType cellType = cellsIt.type();
2678       bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2679       getReverseVector( cellType, swapVec );
2680
2681       elemIt = elems->begin(), elemEnd = elems->end();
2682       for ( ; elemIt != elemEnd; elemIt++ )
2683       {
2684         // GIBI connectivity -> MED one
2685         if( isQuadratic )
2686           ConvertQuadratic( cellType, *elemIt );
2687
2688         // reverse faces
2689         if ( elemIt->_reverse )
2690           reverse ( *elemIt, swapVec );
2691       }
2692     }
2693   }
2694
2695   orientVolumes();
2696 }
2697
2698 //================================================================================
2699 /*!
2700  * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
2701  */
2702 //================================================================================
2703
2704 void IntermediateMED::orientFaces3D()
2705 {
2706   // fill map of links and their faces
2707   set<const Cell*> faces;
2708   map<const Cell*, Group*> fgm;
2709   map<Link, list<const Cell*> > linkFacesMap;
2710   map<Link, list<const Cell*> >::iterator lfIt, lfIt2;
2711
2712   for (size_t i=0; i!=_groups.size(); ++i)
2713     {
2714       Group& grp = _groups[i];
2715       if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
2716         for ( size_t j = 0; j < grp._cells.size(); ++j )
2717           if ( faces.insert( grp._cells[j] ).second )
2718             {
2719               for ( size_t k = 0; k < grp._cells[j]->_nodes.size(); ++k )
2720                 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
2721               fgm.insert( make_pair( grp._cells[j], &grp ));
2722             }
2723     }
2724   // dump linkFacesMap
2725   //     for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
2726   //       cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << endl;
2727   //       list<const Cell*> & fList = lfIt->second;
2728   //       list<const Cell*>::iterator fIt = fList.begin();
2729   //       for ( ; fIt != fList.end(); fIt++ )
2730   //         cout << "\t" << **fIt << fgm[*fIt]->nom << endl;
2731   //     }
2732
2733   // Each oriented link must appear in one face only, else a face is reversed.
2734
2735   queue<const Cell*> faceQueue; /* the queue contains well oriented faces
2736                                      whose neighbors orientation is to be checked */
2737   bool manifold = true;
2738   while ( !linkFacesMap.empty() )
2739     {
2740       if ( faceQueue.empty() )
2741         {
2742           assert( !linkFacesMap.begin()->second.empty() );
2743           faceQueue.push( linkFacesMap.begin()->second.front() );
2744         }
2745       while ( !faceQueue.empty() )
2746         {
2747           const Cell* face = faceQueue.front();
2748           faceQueue.pop();
2749
2750           // loop on links of <face>
2751           for ( int i = 0; i < (int)face->_nodes.size(); ++i )
2752             {
2753               Link link = face->link( i );
2754               // find the neighbor faces
2755               lfIt = linkFacesMap.find( link );
2756               int nbFaceByLink = 0;
2757               list< const Cell* > ml;
2758               if ( lfIt != linkFacesMap.end() )
2759                 {
2760                   list<const Cell*> & fList = lfIt->second;
2761                   list<const Cell*>::iterator fIt = fList.begin();
2762                   assert( fIt != fList.end() );
2763                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2764                     {
2765                       ml.push_back( *fIt );
2766                       if ( *fIt != face ) // wrongly oriented neighbor face
2767                         {
2768                           const Cell* badFace = *fIt;
2769                           // reverse and remove badFace from linkFacesMap
2770                           for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
2771                             {
2772                               Link badlink = badFace->link( j );
2773                               if ( badlink == link ) continue;
2774                               lfIt2 = linkFacesMap.find( badlink );
2775                               if ( lfIt2 != linkFacesMap.end() )
2776                                 {
2777                                   list<const Cell*> & ff = lfIt2->second;
2778                                   list<const Cell*>::iterator lfIt3 = find( ff.begin(), ff.end(), badFace );
2779                                   // check if badFace has been found,
2780                                   // else we can't erase it
2781                                   // case of degenerated face in edge
2782                                   if (lfIt3 != ff.end())
2783                                     {
2784                                       ff.erase( lfIt3 );
2785                                       if ( ff.empty() )
2786                                         linkFacesMap.erase( lfIt2 );
2787                                     }
2788                                 }
2789                             }
2790                           badFace->_reverse = true; // reverse
2791                           //INFOS_MED( "REVERSE " << *badFace );
2792                           faceQueue.push( badFace );
2793                         }
2794                     }
2795                   linkFacesMap.erase( lfIt );
2796                 }
2797               // add good neighbors to the queue
2798               Link revLink( link.second, link.first );
2799               lfIt = linkFacesMap.find( revLink );
2800               if ( lfIt != linkFacesMap.end() )
2801                 {
2802                   list<const Cell*> & fList = lfIt->second;
2803                   list<const Cell*>::iterator fIt = fList.begin();
2804                   for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2805                     {
2806                       ml.push_back( *fIt );
2807                       if ( *fIt != face )
2808                         faceQueue.push( *fIt );
2809                     }
2810                   linkFacesMap.erase( lfIt );
2811                 }
2812               if ( nbFaceByLink > 2 )
2813                 {
2814                   if ( manifold )
2815                     {
2816                       list<const Cell*>::iterator ii = ml.begin();
2817                       cout << nbFaceByLink << " faces by 1 link:" << endl;
2818                       for( ; ii!= ml.end(); ii++ )
2819                         cout << "in sub-mesh <" << fgm[ *ii ]->_name << "> " << **ii << endl;
2820                     }
2821                   manifold = false;
2822                 }
2823             } // loop on links of the being checked face
2824         } // loop on the face queue
2825     } // while ( !linkFacesMap.empty() )
2826
2827   if ( !manifold )
2828     cout << " -> Non manifold mesh, faces orientation may be incorrect" << endl;
2829 }
2830
2831 //================================================================================
2832 /*!
2833  * \brief Orient volumes according to MED conventions:
2834  * normal of a bottom (first) face should be outside
2835  */
2836 //================================================================================
2837
2838 void IntermediateMED::orientVolumes()
2839 {
2840   set<Cell>::const_iterator elemIt, elemEnd;
2841   vector< pair<int,int> > swapVec;
2842
2843   CellsByDimIterator cellsIt( *this, 3 );
2844   while ( const set<Cell > * elems = cellsIt.nextType() )
2845     {
2846       TCellType cellType = cellsIt.type();
2847       elemIt = elems->begin(), elemEnd = elems->end();
2848       int nbBottomNodes = 0;
2849       switch ( cellType )
2850         {
2851         case NORM_TETRA4:
2852         case NORM_TETRA10:
2853         case NORM_PENTA6:
2854         case NORM_PENTA15:
2855           nbBottomNodes = 3; break;
2856         case NORM_PYRA5:
2857         case NORM_PYRA13:
2858         case NORM_HEXA8:
2859         case NORM_HEXA20:
2860           nbBottomNodes = 4; break;
2861         default: continue;
2862         }
2863       getReverseVector( cellType, swapVec );
2864
2865       for ( ; elemIt != elemEnd; elemIt++ )
2866         {
2867           // find a normal to the bottom face
2868           const double* n[4];
2869           n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
2870           n[1] = nodeCoords( elemIt->_nodes[1]);
2871           n[2] = nodeCoords( elemIt->_nodes[2]);
2872           n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
2873           double vec01[3]; // vector n[0]-n[1]
2874           vec01[0] = n[1][0] - n[0][0];
2875           vec01[1] = n[1][1] - n[0][1];
2876           vec01[2] = n[1][2] - n[0][2];
2877           double vec02 [3]; // vector n[0]-n[2]
2878           vec02[0] = n[2][0] - n[0][0];
2879           vec02[1] = n[2][1] - n[0][1];
2880           vec02[2] = n[2][2] - n[0][2];
2881           double normal [3]; // vec01 ^ vec02
2882           normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
2883           normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
2884           normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
2885           // check if the 102 angle is convex
2886           if ( nbBottomNodes > 3 )
2887             {
2888               const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
2889               double vec03 [3];  // vector n[0]-n3
2890               vec03[0] = n3[0] - n[0][0];
2891               vec03[1] = n3[1] - n[0][1];
2892               vec03[2] = n3[2] - n[0][2];
2893               if ( fabs( normal[0]+normal[1]+normal[2] ) <= numeric_limits<double>::max() ) // vec01 || vec02
2894                 {
2895                   normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
2896                   normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
2897                   normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
2898                 }
2899               else
2900                 {
2901                   double vec [3]; // normal ^ vec01
2902                   vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
2903                   vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
2904                   vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
2905                   double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
2906                   if ( dot2 < 0 ) // concave -> reverse normal
2907                     {
2908                       normal[0] *= -1;
2909                       normal[1] *= -1;
2910                       normal[2] *= -1;
2911                     }
2912                 }
2913             }
2914           // direction from top to bottom
2915           double tbDir[3];
2916           tbDir[0] = n[0][0] - n[3][0];
2917           tbDir[1] = n[0][1] - n[3][1];
2918           tbDir[2] = n[0][2] - n[3][2];
2919
2920           // compare 2 directions: normal and top-bottom
2921           double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
2922           if ( dot < 0. ) // need reverse
2923             reverse( *elemIt, swapVec );
2924
2925         } // loop on volumes of one geometry
2926     } // loop on 3D geometry types
2927
2928 }
2929
2930 //================================================================================
2931 /*!
2932  * \brief Assign new IDs to nodes by skipping not used nodes and return their number
2933  */
2934 //================================================================================
2935
2936 int NodeContainer::numberNodes()
2937 {
2938   int id = 1;
2939   for ( size_t i = 0; i < _nodes.size(); ++i )
2940     for ( size_t j = 0; j < _nodes[i].size(); ++j )
2941       if ( _nodes[i][j].isUsed() )
2942         _nodes[i][j]._number = id++;
2943   return id-1;
2944 }
2945
2946
2947 //================================================================================
2948 /*!
2949  * \brief Assign new IDs to elements
2950  */
2951 //================================================================================
2952
2953 void IntermediateMED::numberElements()
2954 {
2955   set<Cell>::const_iterator elemIt, elemEnd;
2956
2957   // numbering _cells of type NORM_POINT1 by node number
2958   {
2959     const set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
2960     elemIt = points.begin(), elemEnd = points.end();
2961     for ( ; elemIt != elemEnd; ++elemIt )
2962       elemIt->_number = elemIt->_nodes[0]->_number;
2963   }
2964
2965   // numbering 1D-3D _cells
2966   for ( int dim = 1; dim <= 3; ++dim )
2967     {
2968       // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
2969       bool ok = true, renumEntity = false;
2970       CellsByDimIterator cellsIt( *this, dim );
2971       int prevNbElems = 0;
2972       while ( const set<Cell> * typeCells = cellsIt.nextType() )
2973         {
2974           TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
2975           for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
2976             {
2977               if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
2978               if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
2979             }
2980           TID typeSize = typeCells->size();
2981           if ( typeSize != maxNumber - minNumber + 1 )
2982             ok = false;
2983           if ( prevNbElems != 0 ) {
2984             if ( minNumber == 1 )
2985               renumEntity = true;
2986             else if ( prevNbElems+1 != (int)minNumber )
2987               ok = false;
2988           }
2989           prevNbElems += typeSize;
2990         }
2991
2992       if ( ok && renumEntity ) // each geom type was numerated separately
2993         {
2994           cellsIt.init( dim );
2995           prevNbElems = cellsIt.nextType()->size(); // no need to renumber the first type
2996           while ( const set<Cell> * typeCells = cellsIt.nextType() )
2997             {
2998               for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
2999                 elemIt->_number += prevNbElems;
3000               prevNbElems += typeCells->size();
3001             }
3002         }
3003       if ( !ok )
3004         {
3005           int cellID=1;
3006           cellsIt.init( dim );
3007           while ( const set<Cell> * typeCells = cellsIt.nextType() )
3008             for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3009               elemIt->_number = cellID++;
3010         }
3011     }
3012 }
3013
3014 //================================================================================
3015 /*!
3016  * \brief Creates coord array
3017  */
3018 //================================================================================
3019
3020 ParaMEDMEM::DataArrayDouble * IntermediateMED::getCoords()
3021 {
3022   DataArrayDouble* coordArray = DataArrayDouble::New();
3023   coordArray->alloc( _nbNodes, _spaceDim );
3024   double * coordPrt = coordArray->getPointer();
3025   for ( int i = 0, nb = _points.size(); i < nb; ++i )
3026     {
3027       Node* n = getNode( i+1 );
3028       if ( n->isUsed() )
3029         {
3030           const double* nCoords = nodeCoords( n );
3031           std::copy( nCoords, nCoords+_spaceDim, coordPrt );
3032           coordPrt += _spaceDim;
3033         }
3034     }
3035   return coordArray;
3036 }
3037
3038 //================================================================================
3039 /*!
3040  * \brief Sets connectivity of elements to the mesh
3041  *  \param mesh - mesh to fill in
3042  *  \param coords - coordinates that must be shared by all meshes of different dim
3043  */
3044 //================================================================================
3045
3046 void IntermediateMED::setConnectivity( ParaMEDMEM::MEDFileUMesh*    mesh,
3047                                        ParaMEDMEM::DataArrayDouble* coords )
3048 {
3049   int meshDim = 0;
3050
3051   mesh->setCoords( coords );
3052
3053   set<Cell>::const_iterator elemIt, elemEnd;
3054   for ( int dim = 3; dim > 0; --dim )
3055   {
3056     CellsByDimIterator dimCells( *this, dim );
3057
3058     int nbOfCells = 0;
3059     while ( const std::set<Cell > * cells = dimCells.nextType() )
3060       nbOfCells += cells->size();
3061     if ( nbOfCells == 0 )
3062       continue;
3063
3064     if ( !meshDim ) meshDim = dim;
3065
3066     MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
3067     dimMesh->setCoords( coords );
3068     dimMesh->setMeshDimension( dim );
3069     dimMesh->allocateCells( nbOfCells );
3070
3071     int prevNbCells = 0;
3072     dimCells.init( dim );
3073     while ( const std::set<Cell > * cells = dimCells.nextType() )
3074       {
3075         // fill connectivity array to take into account order of elements in the sauv file
3076         const int nbCellNodes = cells->begin()->_nodes.size();
3077         vector< TID > connectivity( cells->size() * nbCellNodes );
3078         int * nodalConnOfCell;
3079         for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
3080           {
3081             const Cell& cell = *elemIt;
3082             const int index = cell._number - 1 - prevNbCells;
3083             nodalConnOfCell = &connectivity[ index * nbCellNodes ];
3084             if ( cell._reverse )
3085               for ( int i = nbCellNodes-1; i >= 0; --i )
3086                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3087             else
3088               for ( int i = 0; i < nbCellNodes; ++i )
3089                 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3090           }
3091         prevNbCells += cells->size();
3092
3093         // fill dimMesh
3094         TCellType cellType = dimCells.type();
3095         nodalConnOfCell = &connectivity[0];
3096         for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
3097           dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
3098       }
3099     dimMesh->finishInsertingCells();
3100     mesh->setMeshAtLevel( dim - meshDim, dimMesh );
3101     dimMesh->decrRef();
3102   }
3103 }
3104
3105 //================================================================================
3106 /*!
3107  * \brief Fill in the mesh with groups
3108  *  \param mesh - mesh to fill in
3109  */
3110 //================================================================================
3111
3112 void IntermediateMED::setGroups( ParaMEDMEM::MEDFileUMesh* mesh )
3113 {
3114   bool isMeshNameSet = false;
3115   const int meshDim = mesh->getMeshDimension();
3116   for ( int dim = 0; dim <= meshDim; ++dim )
3117     {
3118       const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
3119
3120       vector<const DataArrayInt *> medGroups;
3121       vector<MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > refGroups;
3122       for ( size_t i = 0; i < _groups.size(); ++i )
3123         {
3124           Group& grp = _groups[i];
3125           if ( (int)getDim( &grp ) != dim &&
3126                grp._groups.empty() ) // to allow groups on diff dims
3127             continue;
3128           // convert only named groups or field supports
3129           if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
3130             continue;
3131           //if ( grp._medGroup ) continue; // already converted
3132
3133           // sort cells by ID and remember their initial order in the group
3134           TCellToOrderMap cell2order;
3135           unsigned orderInGroup = 0;
3136           vector< Group* > groupVec;
3137           if ( grp._groups.empty() ) groupVec.push_back( & grp );
3138           else                       groupVec = grp._groups;
3139           for ( size_t iG = 0; iG < groupVec.size(); ++iG )
3140             {
3141               Group* aG = groupVec[ iG ];
3142               if ( (int)getDim( aG ) != dim )
3143                 continue;
3144               for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
3145                 cell2order.insert( cell2order.end(), make_pair( aG->_cells[iC], orderInGroup++ ));
3146             }
3147           if ( cell2order.empty() )
3148             continue;
3149           bool isSelfIntersect = ( orderInGroup != cell2order.size() );
3150           if ( isSelfIntersect ) // self intersecting group
3151             {
3152               ostringstream msg;
3153               msg << "Self intersecting sub-mesh: id = " << i+1
3154                   << ", name = |" << grp._name << "|" << endl
3155                   << " nb unique elements = " << cell2order.size() << endl
3156                   << " total nb elements  = " << orderInGroup;
3157               if ( grp._isProfile )
3158                 {
3159                   THROW_IK_EXCEPTION( msg.str() );
3160                 }
3161               else
3162                 {
3163                   cout << msg.str() << endl;
3164                 }
3165             }
3166           // create a med group
3167           grp._medGroup = DataArrayInt::New();
3168           grp._medGroup->setName( grp._name.c_str() );
3169           grp._medGroup->alloc( cell2order.size(), /*nbOfCompo=*/1 );
3170           int * idsPrt = grp._medGroup->getPointer();
3171           TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
3172           for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
3173             *idsPrt++ = (*cell2orderIt).first->_number - 1;
3174
3175           // try to set the mesh name
3176           if ( !isMeshNameSet &&
3177                dim == meshDim &&
3178                !grp._name.empty() &&
3179                grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
3180             {
3181               mesh->setName( grp._name.c_str() );
3182               isMeshNameSet = true;
3183             }
3184           if ( !grp._name.empty() )
3185             {
3186               medGroups.push_back( grp._medGroup );
3187             }
3188           // set relocation table
3189           setRelocationTable( &grp, cell2order );
3190
3191           // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
3192           // and several names (pile 27) refer (pile 10) to this group.
3193           // We create a copy of this group per each named reference
3194           for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
3195             if ( !grp._refNames[ iRef ].empty() )
3196               {
3197                 refGroups.push_back( grp._medGroup->deepCpy() );
3198                 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
3199                 medGroups.push_back( refGroups.back() );
3200               }
3201         }
3202       mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
3203     }
3204 }
3205
3206 //================================================================================
3207 /*!
3208  * \brief Return true if the group is on all elements and return its relative dimension
3209  */
3210 //================================================================================
3211
3212 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
3213 {
3214   int dim = getDim( grp );
3215
3216   int nbElems = 0;
3217   CellsByDimIterator dimCells( *this, dim );
3218   while ( const set<Cell > * cells = dimCells.nextType() )
3219     nbElems += cells->size();
3220
3221   const bool onAll = ( nbElems == grp->size() );
3222
3223   if ( dim == 0 )
3224     dimRel = 0;
3225   else
3226     {
3227       int meshDim = 3;
3228       for ( ; meshDim > 0; --meshDim )
3229         {
3230           dimCells.init( meshDim );
3231           if ( dimCells.nextType() )
3232             break;
3233         }
3234       dimRel = dim - meshDim;
3235     }
3236   return onAll;
3237 }
3238
3239 //================================================================================
3240 /*!
3241  * \brief Makes fields from own data
3242  */
3243 //================================================================================
3244
3245 ParaMEDMEM::MEDFileFields * IntermediateMED::makeMEDFileFields(ParaMEDMEM::MEDFileUMesh* mesh)
3246 {
3247   if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
3248
3249   // set long names
3250   set< string > usedFieldNames;
3251   setFieldLongNames(usedFieldNames);
3252
3253   MEDFileFields* fields = MEDFileFields::New();
3254
3255   for ( size_t i = 0; i < _nodeFields.size(); ++i )
3256     setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
3257
3258   for ( size_t i = 0; i < _cellFields.size(); ++i )
3259     setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
3260
3261   return fields;
3262 }
3263
3264 //================================================================================
3265 /*!
3266  * \brief Make med fields from a SauvUtilities::DoubleField
3267  */
3268 //================================================================================
3269
3270 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
3271                                  ParaMEDMEM::MEDFileFields*  medFields,
3272                                  ParaMEDMEM::MEDFileUMesh*   mesh,
3273                                  const TID                   castemID,
3274                                  set< string >&              usedFieldNames)
3275 {
3276   bool sameNbGauss = true;
3277   if ( !fld || !fld->isMedCompatible( sameNbGauss )) return;
3278
3279   if ( !sameNbGauss )
3280     fld->splitSubWithDiffNbGauss();
3281
3282   // if ( !fld->hasCommonSupport() ):
3283   //     each sub makes MEDFileFieldMultiTS
3284   // else:
3285   //     unite several subs into a MEDCouplingFieldDouble
3286
3287   const bool uniteSubs = fld->hasCommonSupport() && sameNbGauss;
3288   if ( !uniteSubs )
3289     cout << "Castem field #" << castemID << " <" << fld->_name
3290          << "> is incompatible with MED format, so we split it into several fields:" << endl;
3291
3292   for ( size_t iSub = 0; iSub < fld->_sub.size(); )
3293     {
3294       // set field name
3295       if ( !uniteSubs || fld->_name.empty() )
3296         makeFieldNewName( usedFieldNames, fld );
3297
3298       // allocate values
3299       DataArrayDouble * values = DataArrayDouble::New();
3300       values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
3301
3302       // set values
3303       double * valPtr = values->getPointer();
3304       if ( uniteSubs )
3305         {
3306           int nbElems = fld->_group->size();
3307           for ( int elemShift = 0; elemShift < nbElems && iSub < fld->_sub.size(); )
3308             elemShift += fld->setValues( valPtr, iSub++, elemShift );
3309           setTS( fld, values, medFields, mesh );
3310         }
3311       else
3312         {
3313           fld->setValues( valPtr, iSub );
3314           setTS( fld, values, medFields, mesh, iSub++ );
3315
3316           cout << fld->_name << " with compoments";
3317           for ( size_t i = 0; i < (size_t)fld->_sub[iSub-1].nbComponents(); ++i )
3318             cout << " " << fld->_sub[iSub-1]._comp_names[ i ];
3319           cout << endl;
3320         }
3321     }
3322 }
3323
3324 //================================================================================
3325 /*!
3326  * \brief Store value array of a field into med fields
3327  */
3328 //================================================================================
3329
3330 void IntermediateMED::setTS( SauvUtilities::DoubleField*  fld,
3331                              ParaMEDMEM::DataArrayDouble* values,
3332                              ParaMEDMEM::MEDFileFields*   medFields,
3333                              ParaMEDMEM::MEDFileUMesh*    mesh,
3334                              const int                    iSub)
3335 {
3336   // treat a field support
3337   const Group* support = fld->getSupport( iSub );
3338   int dimRel;
3339   const bool onAll = isOnAll( support, dimRel );
3340   if ( !onAll && support->_name.empty() )
3341     {
3342       const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
3343       support->_medGroup->setName( support->_name.c_str() );
3344     }
3345
3346   // make and fill a time-stamp
3347
3348   MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType( iSub ),
3349                                                                     fld->getMedTimeDisc() );
3350   timeStamp->setName( fld->_name.c_str() );
3351   timeStamp->setDescription( fld->_description.c_str() );
3352   // set the mesh
3353   if ( onAll )
3354     {
3355       MEDCouplingAutoRefCountObjectPtr
3356         < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3357       timeStamp->setMesh( dimMesh );
3358     }
3359   else if ( timeStamp->getTypeOfField() == ParaMEDMEM::ON_NODES )
3360     {
3361       DataArrayDouble * coo = mesh->getCoords();
3362       MEDCouplingAutoRefCountObjectPtr
3363         <DataArrayDouble> subCoo = coo->selectByTupleId(support->_medGroup->begin(),
3364                                                         support->_medGroup->end()); 
3365       MEDCouplingAutoRefCountObjectPtr< MEDCouplingUMesh > nodeSubMesh =
3366         MEDCouplingUMesh::Build0DMeshFromCoords( subCoo );
3367       timeStamp->setMesh( nodeSubMesh );
3368     }
3369   else
3370     {
3371       MEDCouplingAutoRefCountObjectPtr
3372         < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3373       MEDCouplingAutoRefCountObjectPtr
3374         <MEDCouplingMesh> subMesh = dimMesh->buildPart(support->_medGroup->begin(),
3375                                                        support->_medGroup->end());
3376       timeStamp->setMesh( subMesh);
3377     }
3378   // set values
3379   for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
3380     values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
3381   timeStamp->setArray( values );
3382   values->decrRef();
3383   // set gauss points
3384   if ( timeStamp->getTypeOfField() == ParaMEDMEM::ON_GAUSS_PT )
3385     {
3386       TGaussDef gaussDef( support->_cellType, fld->_sub[iSub].nbGauss() );
3387       timeStamp->setGaussLocalizationOnType( support->_cellType,
3388                                              gaussDef.myRefCoords,
3389                                              gaussDef.myCoords,
3390                                              gaussDef.myWeights );
3391     }
3392   // get a field to add the time-stamp
3393   bool isNewMedField = false;
3394   if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
3395     {
3396       fld->_curMedField = MEDFileFieldMultiTS::New();
3397       isNewMedField = true;
3398     }
3399
3400   // set an order
3401   const int nbTS = fld->_curMedField->getNumberOfTS();
3402   if ( nbTS > 0 )
3403     timeStamp->setOrder( nbTS );
3404
3405   // add the time-stamp
3406   timeStamp->checkCoherency();
3407   if ( onAll )
3408     fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
3409   else
3410     fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
3411   timeStamp->decrRef();
3412
3413   if ( isNewMedField ) // timeStamp must be added before this
3414     {
3415       medFields->pushField( fld->_curMedField );
3416     }
3417 }
3418
3419 //================================================================================
3420 /*!
3421  * \brief Make a new unique name for a field
3422  */
3423 //================================================================================
3424
3425 void IntermediateMED::makeFieldNewName(std::set< std::string >&    usedNames,
3426                                        SauvUtilities::DoubleField* fld )
3427 {
3428   string base = fld->_name;
3429   if ( base.empty() )
3430     {
3431       base = "F_";
3432     }
3433   else
3434     {
3435       string::size_type pos = base.rfind('_');
3436       if ( pos != string::npos )
3437         base = base.substr( 0, pos+1 );
3438       else
3439         base += '_';
3440     }
3441
3442   int i = 1;
3443   do
3444     {
3445       fld->_name = base + SauvUtilities::toString( i++ );
3446     }
3447   while( !usedNames.insert( fld->_name ).second );
3448 }
3449
3450 //================================================================================
3451 /*!
3452  * \brief Split sub-components with different nb of gauss points into several sub-components
3453  *  \param [in,out] fld - a field to split if necessary
3454  */
3455 //================================================================================
3456
3457 void DoubleField::splitSubWithDiffNbGauss()
3458 {
3459   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
3460     {
3461       if ( _sub[iSub].isSameNbGauss() ) continue;
3462
3463       _sub.insert( _sub.begin() + iSub + 1, 1, _Sub_data() );
3464       _Sub_data & subToSplit = _sub[iSub];
3465       _Sub_data & subNew     = _sub[iSub+1];
3466       size_t iDiff = 1;
3467       while ( subToSplit._nb_gauss[ 0 ] == subToSplit._nb_gauss[ iDiff ] )
3468         ++iDiff;
3469       subNew._support = subToSplit._support;
3470       subNew._comp_names.assign( subToSplit._comp_names.begin() + iDiff,
3471                                  subToSplit._comp_names.end() );
3472       subNew._nb_gauss.assign  ( subToSplit._nb_gauss.begin() + iDiff,
3473                                  subToSplit._nb_gauss.end() );
3474       subToSplit._comp_names.resize( iDiff );
3475       subToSplit._nb_gauss.resize  ( iDiff );
3476     }
3477 }
3478
3479 //================================================================================
3480 /*!
3481  * \brief Return a vector ready to fill in
3482  */
3483 //================================================================================
3484
3485 std::vector< double >& DoubleField::addComponent( int nb_values )
3486 {
3487   _comp_values.push_back( std::vector< double >() );
3488   std::vector< double >& res = _comp_values.back();
3489   res.resize( nb_values );
3490   return res;
3491 }
3492
3493 DoubleField::~DoubleField()
3494 {
3495   if(_curMedField)
3496     _curMedField->decrRef();
3497 }
3498
3499 //================================================================================
3500 /*!
3501  * \brief Returns a supporting group
3502  */
3503 //================================================================================
3504
3505 const Group* DoubleField::getSupport( const int iSub ) const
3506 {
3507   return _group ? _group : _sub[iSub]._support;
3508 }
3509
3510 //================================================================================
3511 /*!
3512  * \brief Return true if each sub-component is a time stamp
3513  */
3514 //================================================================================
3515
3516 bool DoubleField::isMultiTimeStamps() const
3517 {
3518   if ( _sub.size() < 2 )
3519     return false;
3520   bool sameSupports = true;
3521   Group* grp1 = _sub[0]._support;
3522   for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
3523     sameSupports = ( grp1 == _sub[i]._support );
3524
3525   return sameSupports;
3526 }
3527
3528 //================================================================================
3529 /*!
3530  * \brief True if the field can be converted into the med field
3531  */
3532 //================================================================================
3533
3534 bool DoubleField::isMedCompatible(bool& sameNbGauss) const
3535 {
3536   for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
3537     {
3538       if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
3539         THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
3540
3541       sameNbGauss = true;
3542       if ( !_sub[iSub].isSameNbGauss() )
3543         {
3544           cout << "Field <" << _name << "> : different nb of gauss points in components" <<endl;
3545           sameNbGauss = false;
3546           //return false;
3547         }
3548     }
3549   return true;
3550 }
3551
3552 //================================================================================
3553 /*!
3554  * \brief return true if all sub-components has same components and same nbGauss
3555  */
3556 //================================================================================
3557
3558 bool DoubleField::hasSameComponentsBySupport() const
3559 {
3560   vector< _Sub_data >::const_iterator sub_data = _sub.begin();
3561   const _Sub_data& first_sub_data = *sub_data;
3562   for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
3563   {
3564     if ( first_sub_data._comp_names != sub_data->_comp_names )
3565       return false; // diff names of components
3566
3567     if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
3568          first_sub_data._support->_cellType == sub_data->_support->_cellType)
3569       return false; // diff nb of gauss points on same cell type
3570   }
3571   return true;
3572 }
3573
3574 //================================================================================
3575 /*!
3576  * \brief Return type of MEDCouplingFieldDouble
3577  */
3578 //================================================================================
3579
3580 ParaMEDMEM::TypeOfField DoubleField::getMedType( const int iSub ) const
3581 {
3582   using namespace INTERP_KERNEL;
3583
3584   const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
3585   if ( _sub[iSub].nbGauss() > 1 )
3586     {
3587       const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
3588       return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
3589     }
3590   else
3591     {
3592       return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
3593     }
3594 }
3595
3596 //================================================================================
3597 /*!
3598  * \brief Return TypeOfTimeDiscretization
3599  */
3600 //================================================================================
3601
3602 ParaMEDMEM::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
3603 {
3604   return ONE_TIME;
3605   // NO_TIME = 4,
3606   // ONE_TIME = 5,
3607   // LINEAR_TIME = 6,
3608   // CONST_ON_TIME_INTERVAL = 7
3609 }
3610
3611 //================================================================================
3612 /*!
3613  * \brief Return nb tuples to be used to allocate DataArrayDouble
3614  */
3615 //================================================================================
3616
3617 int DoubleField::getNbTuples( const int iSub ) const
3618 {
3619   int nb = 0;
3620   if ( hasCommonSupport() && !_group->_groups.empty() )
3621     for ( size_t i = 0; i < _group->_groups.size(); ++i )
3622       nb += _sub[i].nbGauss() * _sub[i]._support->size();
3623   else
3624     nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
3625   return nb;
3626 }
3627
3628 //================================================================================
3629 /*!
3630  * \brief Store values of a sub-component and return nb of elements in the iSub
3631  */
3632 //================================================================================
3633
3634 int DoubleField::setValues( double * valPtr, const int iSub, const int elemShift ) const
3635 {
3636   // find values for iSub
3637   int iComp = 0;
3638   for ( int iS = 0; iS < iSub; ++iS )
3639     iComp += _sub[iS].nbComponents();
3640   const vector< double > * compValues = &_comp_values[ iComp ];
3641
3642   // Set values
3643
3644   const vector< unsigned >& relocTable = getSupport( iSub )->_relocTable;
3645
3646   const int nbElems      = _sub[iSub]._support->size();
3647   const int nbGauss      = _sub[iSub].nbGauss();
3648   const int nbComponents = _sub[iSub].nbComponents();
3649   const int nbValsByElem = nbComponents * nbGauss;
3650
3651   // check nb values
3652   int nbVals = 0;
3653   for ( iComp = 0; iComp < nbComponents; ++iComp )
3654     nbVals += compValues[iComp].size();
3655   const bool isConstField = ( nbVals == nbComponents ); // one value per component (issue 22321)
3656   if ( !isConstField && nbVals != nbElems * nbValsByElem )
3657     THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
3658
3659   // compute nb values in previous subs
3660   int valsShift = 0;
3661   for ( int iS = iSub-1, shift = elemShift; shift > 0; --iS)
3662     {
3663       int nbE = _sub[iS]._support->size();
3664       shift -= nbE;
3665       valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
3666     }
3667
3668   if ( isConstField )
3669     for ( int iE = 0; iE < nbElems; ++iE )
3670       {
3671         int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3672         for ( iComp = 0; iComp < nbComponents; ++iComp )
3673           valPtr[ iMed + iComp ] = compValues[iComp][ 0 ];
3674       }
3675   else
3676     for ( int iE = 0; iE < nbElems; ++iE )
3677       {
3678         int iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3679         for ( iComp = 0; iComp < nbComponents; ++iComp )
3680           for ( int iG = 0; iG < nbGauss; ++iG )
3681             valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
3682       }
3683   return nbElems;
3684 }
3685
3686 //================================================================================
3687 /*!
3688  * \brief Destructor of IntermediateMED
3689  */
3690 //================================================================================
3691
3692 IntermediateMED::~IntermediateMED()
3693 {
3694   for ( size_t i = 0; i < _nodeFields.size(); ++i )
3695     if ( _nodeFields[i] )
3696       delete _nodeFields[i];
3697   _nodeFields.clear();
3698
3699   for ( size_t i = 0; i < _cellFields.size(); ++i )
3700     if ( _cellFields[i] )
3701       delete _cellFields[i];
3702   _cellFields.clear();
3703
3704   for ( size_t i = 0; i < _groups.size(); ++i )
3705     if ( _groups[i]._medGroup )
3706       _groups[i]._medGroup->decrRef();
3707 }
3708
3709 //================================================================================
3710 /*!
3711  * \brief CellsByDimIterator constructor
3712  */
3713 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
3714 {
3715   myImed = & medi;
3716   init( dimm );
3717 }
3718 /*!
3719  * \brief Initialize iteration on cells of given dimention
3720  */
3721 void CellsByDimIterator::init(const int  dimm)
3722 {
3723   myCurType = -1;
3724   myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
3725   myDim = dimm;
3726 }
3727 /*!
3728  * \brief return next set of Cell's of required dimension
3729  */
3730 const std::set< Cell > * CellsByDimIterator::nextType()
3731 {
3732   while ( ++myCurType < myTypeEnd )
3733     if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
3734       return & myImed->_cellsByType[myCurType];
3735   return 0;
3736 }
3737 /*!
3738  * \brief return dimension of cells returned by the last or further next()
3739  */
3740 int CellsByDimIterator::dim(const bool last) const
3741 {
3742   int typp = myCurType;
3743   if ( !last )
3744     while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
3745       ++typp;
3746   return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
3747 }
3748 // END CellsByDimIterator ========================================================
3749