1 // Copyright (C) 2007-2023 CEA, EDF
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, or (at your option) any later version.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : SauvMedConvertor.cxx
20 // Created : Tue Aug 16 14:43:20 2011
21 // Author : Edward AGAPOV (eap)
24 #include "SauvMedConvertor.hxx"
26 #include "CellModel.hxx"
27 #include "MEDFileMesh.hxx"
28 #include "MEDFileField.hxx"
29 #include "MEDFileData.hxx"
30 #include "MEDCouplingFieldDouble.hxx"
49 #include <rpc/types.h>
53 using namespace SauvUtilities;
54 using namespace MEDCoupling;
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
62 using namespace INTERP_KERNEL;
64 const size_t MaxMedCellType = NORM_ERROR;
65 const size_t NbGibiCellTypes = 47;
66 const TCellType GibiTypeToMed[NbGibiCellTypes] =
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
80 //================================================================================
82 * \brief Return dimension of a group
84 //================================================================================
86 unsigned getDim( const Group* grp )
88 return SauvUtilities::getDimension( grp->_groups.empty() ? grp->_cellType : grp->_groups[0]->_cellType );
91 //================================================================================
93 * \brief Converts connectivity of quadratic elements
95 //================================================================================
97 inline void ConvertQuadratic( const INTERP_KERNEL::NormalizedCellType type,
100 if ( const int * conn = getGibi2MedQuadraticInterlace( type ))
102 Cell* ma = const_cast<Cell*>(&aCell);
103 std::vector< Node* > new_nodes( ma->_nodes.size() );
104 for (std:: size_t i = 0; i < new_nodes.size(); ++i )
105 new_nodes[ i ] = ma->_nodes[ conn[ i ]];
106 ma->_nodes.swap( new_nodes );
110 //================================================================================
112 * \brief Returns a vector of pairs of node indices to inverse a med volume element
114 //================================================================================
116 void getReverseVector (const INTERP_KERNEL::NormalizedCellType type,
117 std::vector<std::pair<int,int> > & swapVec )
125 swapVec[0] = std::make_pair( 1, 2 );
129 swapVec[0] = std::make_pair( 1, 3 );
133 swapVec[0] = std::make_pair( 1, 2 );
134 swapVec[1] = std::make_pair( 4, 5 );
138 swapVec[0] = std::make_pair( 1, 3 );
139 swapVec[1] = std::make_pair( 5, 7 );
143 swapVec[0] = std::make_pair( 1, 2 );
144 swapVec[1] = std::make_pair( 4, 6 );
145 swapVec[2] = std::make_pair( 8, 9 );
149 swapVec[0] = std::make_pair( 1, 3 );
150 swapVec[1] = std::make_pair( 5, 8 );
151 swapVec[2] = std::make_pair( 6, 7 );
152 swapVec[3] = std::make_pair( 10, 12 );
156 swapVec[0] = std::make_pair( 1, 2 );
157 swapVec[1] = std::make_pair( 4, 5 );
158 swapVec[2] = std::make_pair( 6, 8 );
159 swapVec[3] = std::make_pair( 9, 11 );
163 swapVec[0] = std::make_pair( 1, 3 );
164 swapVec[1] = std::make_pair( 5, 7 );
165 swapVec[2] = std::make_pair( 8, 11 );
166 swapVec[3] = std::make_pair( 9, 10 );
167 swapVec[4] = std::make_pair( 12, 15 );
168 swapVec[5] = std::make_pair( 13, 14 );
169 swapVec[6] = std::make_pair( 17, 19 );
171 // case NORM_SEG3: no need to reverse edges
172 // swapVec.resize(1);
173 // swapVec[0] = std::make_pair( 1, 2 );
177 swapVec[0] = std::make_pair( 1, 2 );
178 swapVec[1] = std::make_pair( 3, 5 );
182 swapVec[0] = std::make_pair( 1, 3 );
183 swapVec[1] = std::make_pair( 4, 7 );
184 swapVec[2] = std::make_pair( 5, 6 );
190 //================================================================================
192 * \brief Inverses element orientation using vector of indices to swap
194 //================================================================================
196 inline void reverse(const Cell & aCell, const std::vector<std::pair<int,int> > & swapVec )
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() )
205 ma->_reverse = false;
207 //================================================================================
209 * \brief Comparator of cells by number used for ordering cells within a med group
211 struct TCellByIDCompare
213 bool operator () (const Cell* i1, const Cell* i2) const
215 return i1->_number < i2->_number;
218 typedef std::map< const Cell*, unsigned, TCellByIDCompare > TCellToOrderMap;
220 //================================================================================
222 * \brief Fill Group::_relocTable if necessary
224 //================================================================================
226 void setRelocationTable( Group* grp, TCellToOrderMap& cell2order )
228 if ( !grp->_isProfile ) return;
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 );
239 grp->_relocTable.resize( cell2order.size() );
240 for ( newOrder = 0, c2oIt = cell2order.begin(); c2oIt != c2oEnd; ++c2oIt, ++newOrder )
241 grp->_relocTable[ c2oIt->second ] = newOrder;
246 namespace // define default GAUSS points
248 typedef std::vector<double> TDoubleVector;
249 typedef double* TCoordSlice;
251 //---------------------------------------------------------------
252 //! Shape function definitions
253 //---------------------------------------------------------------
258 TDoubleVector myRefCoord;
260 TShapeFun(TInt theDim = 0, TInt theNbRef = 0)
261 :myDim(theDim),myNbRef(theNbRef),myRefCoord(theNbRef*theDim) {}
263 TInt GetNbRef() const { return myNbRef; }
265 TCoordSlice GetCoord(std::size_t theRefId) { return &myRefCoord[0] + theRefId * myDim; }
267 //---------------------------------------------------------------
269 * \brief Description of family of integration points
271 //---------------------------------------------------------------
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>
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
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
291 TGaussDef(const int geomType, const int nbPoints, const int variant=1);
293 int dim() const { return SauvUtilities::getDimension( NormalizedCellType( myType )); }
294 std::size_t nbPoints() const { return myWeights.capacity(); }
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; }
302 struct TSeg2a: TShapeFun {
305 struct TSeg3a: TShapeFun {
308 struct TTria3a: TShapeFun {
311 struct TTria6a: TShapeFun {
314 struct TTria3b: TShapeFun {
317 struct TTria6b: TShapeFun {
320 struct TQuad4a: TShapeFun {
323 struct TQuad8a: TShapeFun {
326 struct TQuad4b: TShapeFun {
329 struct TQuad8b: TShapeFun {
332 struct TTetra4a: TShapeFun {
335 struct TTetra10a: TShapeFun {
338 struct TTetra4b: TShapeFun {
341 struct TTetra10b: TShapeFun {
344 struct THexa8a: TShapeFun {
347 struct THexa20a: TShapeFun {
348 THexa20a(TInt theDim = 3, TInt theNbRef = 20);
350 struct THexa27a: THexa20a {
353 struct THexa8b: TShapeFun {
356 struct THexa20b: TShapeFun {
357 THexa20b(TInt theDim = 3, TInt theNbRef = 20);
359 struct TPenta6a: TShapeFun {
362 struct TPenta6b: TShapeFun {
365 struct TPenta15a: TShapeFun {
368 struct TPenta15b: TShapeFun {
371 struct TPyra5a: TShapeFun {
374 struct TPyra5b: TShapeFun {
377 struct TPyra13a: TShapeFun {
380 struct TPyra13b: TShapeFun {
384 void TGaussDef::add(const double x, const double weight)
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 );
393 void TGaussDef::add(const double x, const double y, const double weight)
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 );
403 void TGaussDef::add(const double x, const double y, const double z, const double weight)
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 );
415 //---------------------------------------------------------------
416 TSeg2a::TSeg2a():TShapeFun(1,2)
418 TInt aNbRef = GetNbRef();
419 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
420 TCoordSlice aCoord = GetCoord(aRefId);
422 case 0: aCoord[0] = -1.0; break;
423 case 1: aCoord[0] = 1.0; break;
427 //---------------------------------------------------------------
428 TSeg3a::TSeg3a():TShapeFun(1,3)
430 TInt aNbRef = GetNbRef();
431 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
432 TCoordSlice aCoord = GetCoord(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;
440 //---------------------------------------------------------------
444 TInt aNbRef = GetNbRef();
445 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
446 TCoordSlice aCoord = GetCoord(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;
454 //---------------------------------------------------------------
455 TTria6a::TTria6a():TShapeFun(2,6)
457 TInt aNbRef = GetNbRef();
458 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
459 TCoordSlice aCoord = GetCoord(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;
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;
471 //---------------------------------------------------------------
475 TInt aNbRef = GetNbRef();
476 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
477 TCoordSlice aCoord = GetCoord(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;
485 //---------------------------------------------------------------
489 TInt aNbRef = GetNbRef();
490 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
491 TCoordSlice aCoord = GetCoord(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;
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;
503 //---------------------------------------------------------------
507 TInt aNbRef = GetNbRef();
508 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
509 TCoordSlice aCoord = GetCoord(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;
518 //---------------------------------------------------------------
522 TInt aNbRef = GetNbRef();
523 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
524 TCoordSlice aCoord = GetCoord(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;
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;
538 //---------------------------------------------------------------
542 TInt aNbRef = GetNbRef();
543 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
544 TCoordSlice aCoord = GetCoord(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;
553 //---------------------------------------------------------------
557 TInt aNbRef = GetNbRef();
558 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
559 TCoordSlice aCoord = GetCoord(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;
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;
573 //---------------------------------------------------------------
574 TTetra4a::TTetra4a():
577 TInt aNbRef = GetNbRef();
578 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
579 TCoordSlice aCoord = GetCoord(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;
588 //---------------------------------------------------------------
589 TTetra10a::TTetra10a():
592 TInt aNbRef = GetNbRef();
593 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
594 TCoordSlice aCoord = GetCoord(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;
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;
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;
611 //---------------------------------------------------------------
612 TTetra4b::TTetra4b():
615 TInt aNbRef = GetNbRef();
616 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
617 TCoordSlice aCoord = GetCoord(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;
626 //---------------------------------------------------------------
627 TTetra10b::TTetra10b():
630 TInt aNbRef = GetNbRef();
631 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
632 TCoordSlice aCoord = GetCoord(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;
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;
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;
649 //---------------------------------------------------------------
653 TInt aNbRef = GetNbRef();
654 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
655 TCoordSlice aCoord = GetCoord(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;
668 //---------------------------------------------------------------
669 THexa20a::THexa20a(TInt theDim, TInt theNbRef):
670 TShapeFun(theDim,theNbRef)
672 std::size_t aNbRef = myRefCoord.size();
673 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
674 TCoordSlice aCoord = GetCoord(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;
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;
700 //---------------------------------------------------------------
701 THexa27a::THexa27a():
704 std::size_t aNbRef = myRefCoord.size();
705 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
706 TCoordSlice aCoord = GetCoord(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;
718 //---------------------------------------------------------------
722 TInt aNbRef = GetNbRef();
723 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
724 TCoordSlice aCoord = GetCoord(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;
737 //---------------------------------------------------------------
738 THexa20b::THexa20b(TInt theDim, TInt theNbRef):
739 TShapeFun(theDim,theNbRef)
741 std::size_t aNbRef = myRefCoord.size();
742 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
743 TCoordSlice aCoord = GetCoord(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;
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;
769 //---------------------------------------------------------------
770 TPenta6a::TPenta6a():
773 std::size_t aNbRef = myRefCoord.size();
774 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
775 TCoordSlice aCoord = GetCoord(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;
786 //---------------------------------------------------------------
787 TPenta6b::TPenta6b():
790 std::size_t aNbRef = myRefCoord.size();
791 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
792 TCoordSlice aCoord = GetCoord(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;
803 //---------------------------------------------------------------
804 TPenta15a::TPenta15a():
807 std::size_t aNbRef = myRefCoord.size();
808 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
809 TCoordSlice aCoord = GetCoord(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;
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;
830 //---------------------------------------------------------------
831 TPenta15b::TPenta15b():
834 std::size_t aNbRef = myRefCoord.size();
835 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
836 TCoordSlice aCoord = GetCoord(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;
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;
857 //---------------------------------------------------------------
861 std::size_t aNbRef = myRefCoord.size();
862 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
863 TCoordSlice aCoord = GetCoord(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;
873 //---------------------------------------------------------------
877 std::size_t aNbRef = myRefCoord.size();
878 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
879 TCoordSlice aCoord = GetCoord(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;
889 //---------------------------------------------------------------
890 TPyra13a::TPyra13a():
893 std::size_t aNbRef = myRefCoord.size();
894 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
895 TCoordSlice aCoord = GetCoord(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;
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;
914 //---------------------------------------------------------------
915 TPyra13b::TPyra13b():
918 std::size_t aNbRef = myRefCoord.size();
919 for(std::size_t aRefId = 0; aRefId < aNbRef; aRefId++){
920 TCoordSlice aCoord = GetCoord(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;
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;
940 * \brief Fill definition of gauss points family
943 TGaussDef::TGaussDef(const int geom, const int nbGauss, const int variant)
946 myCoords .reserve( nbGauss * dim() );
947 myWeights.reserve( nbGauss );
953 if (geom == NORM_SEG2) setRefCoords( TSeg2a() );
954 else setRefCoords( TSeg3a() );
957 add( 0.0, 2.0 ); break;
960 const double a = 0.577350269189626;
962 add( -a, 1.0 ); break;
965 const double a = 0.774596669241;
966 const double P1 = 1./1.8;
967 const double P2 = 1./1.125;
973 const double a = 0.339981043584856, b = 0.861136311594053;
974 const double P1 = 0.652145154862546, P2 = 0.347854845137454 ;
978 add( -b, P2 ); break;
981 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for SEG"<<nbGauss);
987 if ( variant == 1 ) {
988 if (geom == NORM_TRI3) setRefCoords( TTria3b() );
989 else setRefCoords( TTria6b() );
992 add( 1/3., 1/3., 1/2. ); break;
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;
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;
1007 const double P1 = 0.11169079483905, P2 = 0.0549758718227661;
1008 const double a = 0.445948490915965, b = 0.091576213509771;
1010 add( 1-2*b, b, P2 );
1011 add( b, 1-2*b, P2 );
1012 add( a, 1-2*a, P1 );
1014 add( 1-2*a, a, P1 ); break;
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. );
1023 add( 1-2*A, A, P1 );
1024 add( A, 1-2*A, P1 );
1026 add( 1-2*B, B, P2 );
1027 add( B, 1-2*B, P2 ); break;
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;
1038 add( 1-2*A, A, P1 );
1039 add( A, 1-2*A, P1 );
1041 add( 1-2*B, B, P2 );
1042 add( B, 1-2*B, P2 );
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;
1051 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 1: "
1055 else if ( variant == 2 ) {
1056 if (geom == NORM_TRI3) setRefCoords( TTria3a() );
1057 else setRefCoords( TTria6a() );
1058 switch ( nbGauss ) {
1060 add( -1/3., -1/3., 2. ); break;
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;
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;
1078 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 2: "
1082 else if ( variant == 3 ) {
1083 if (geom == NORM_TRI3) setRefCoords( TTria3b() );
1084 else setRefCoords( TTria6b() );
1085 switch ( nbGauss ) {
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;
1093 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TRIA, variant 3: "
1101 if ( variant == 1 ) {
1102 if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1103 else setRefCoords( TQuad8b() );
1104 switch ( nbGauss ) {
1106 add( 0, 0, 4 ); break;
1109 const double a = 1/sqrt(3.);
1113 add( -a, a, 1 ); break;
1115 case 5: { // out from the 3 specs
1116 const double a = 0.774596669241483;
1121 add( 0, 0, 2.0 ); break;
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;
1136 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1140 else if ( variant == 2 ) {
1141 if (geom == NORM_QUAD4) setRefCoords( TQuad4a() );
1142 else setRefCoords( TQuad8a() );
1143 switch ( nbGauss ) {
1145 const double a = 1/sqrt(3.);
1149 add( a, a, 1 ); break;
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;
1164 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 1: "
1168 else if ( variant == 3 ) {
1169 if (geom == NORM_QUAD4) setRefCoords( TQuad4b() );
1170 else setRefCoords( TQuad8b() );
1171 switch ( nbGauss ) {
1173 const double a = 3/sqrt(3.);
1177 add( a, a, 1 ); break;
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;
1183 add( -a, 0., c1c2 );
1185 add( 0., -a, c1c2 );
1190 add( a, a, c12 ); break;
1193 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for QUAD, variant 3: "
1201 if (geom == NORM_TETRA4) setRefCoords( TTetra4a() );
1202 else setRefCoords( TTetra10a() );
1203 switch ( nbGauss ) {
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;
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;
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.);
1243 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for TETRA: "<<nbGauss);
1249 if (geom == NORM_PYRA5) setRefCoords( TPyra5a() );
1250 else setRefCoords( TPyra13a() );
1251 switch ( nbGauss ) {
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;
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;
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
1291 add( 0., 0., z, b6 ); // 6
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
1322 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PYRA: "<<nbGauss);
1327 if (geom == NORM_PENTA6) setRefCoords( TPenta6a() );
1328 else setRefCoords( TPenta15a() );
1329 switch ( nbGauss ) {
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;
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;
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 );//___
1380 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1386 if (geom == NORM_HEXA8) setRefCoords( THexa8a() );
1387 else setRefCoords( THexa20a() );
1388 switch ( nbGauss ) {
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;
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
1434 THROW_IK_EXCEPTION("TGaussDef: Invalid nb of gauss points for PENTA: " <<nbGauss);
1439 THROW_IK_EXCEPTION("TGaussDef: unexpected EGeometrieElement: "<< geom);
1442 if ( myWeights.capacity() != myWeights.size() )
1443 THROW_IK_EXCEPTION("TGaussDef: Not all gauss points defined");
1447 //================================================================================
1449 * \brief Return dimension for the given cell type
1451 //================================================================================
1453 unsigned SauvUtilities::getDimension( INTERP_KERNEL::NormalizedCellType type )
1455 return type == NORM_ERROR ? -1 : INTERP_KERNEL::CellModel::GetCellModel( type ).getDimension();
1458 //================================================================================
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
1463 //================================================================================
1465 const int * SauvUtilities::getGibi2MedQuadraticInterlace( INTERP_KERNEL::NormalizedCellType type )
1467 static std::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};
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;
1486 return conn[ type ];
1489 //================================================================================
1491 * \brief Avoid coping sortedNodeIDs
1493 //================================================================================
1495 Cell::Cell(const Cell& ma)
1496 : _nodes(ma._nodes), _reverse(ma._reverse), _sortedNodeIDs(0), _number(ma._number)
1498 if ( ma._sortedNodeIDs )
1500 _sortedNodeIDs = new TID[ _nodes.size() ];
1501 std::copy( ma._sortedNodeIDs, ma._sortedNodeIDs + _nodes.size(), _sortedNodeIDs );
1505 //================================================================================
1507 * \brief Rerturn the i-th link of face
1509 //================================================================================
1511 SauvUtilities::Link Cell::link(int i) const
1513 std::size_t i2 = ( i + 1 ) % _nodes.size();
1515 return std::make_pair( _nodes[i2]->_number, _nodes[i]->_number );
1517 return std::make_pair( _nodes[i]->_number, _nodes[i2]->_number );
1520 //================================================================================
1522 * \brief creates if needed and return _sortedNodeIDs
1524 //================================================================================
1526 const TID* Cell::getSortedNodes() const
1528 if ( !_sortedNodeIDs )
1530 size_t l=_nodes.size();
1531 _sortedNodeIDs = new TID[ l ];
1533 for (size_t i=0; i!=l; ++i)
1534 _sortedNodeIDs[i]=_nodes[i]->_number;
1535 std::sort( _sortedNodeIDs, _sortedNodeIDs + l );
1537 return _sortedNodeIDs;
1540 //================================================================================
1542 * \brief Compare sorted ids of cell nodes
1544 //================================================================================
1546 bool Cell::operator< (const Cell& ma) const
1548 if ( _nodes.size() == 1 )
1549 return _nodes[0] < ma._nodes[0];
1551 const TID* v1 = getSortedNodes();
1552 const TID* v2 = ma.getSortedNodes();
1553 for ( const TID* vEnd = v1 + _nodes.size(); v1 < vEnd; ++v1, ++v2 )
1559 //================================================================================
1561 * \brief Dump a Cell
1563 //================================================================================
1565 std::ostream& SauvUtilities::operator<< (std::ostream& os, const SauvUtilities::Cell& ma)
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;
1571 os << " > sortedNodes: ";
1572 if ( ma._sortedNodeIDs ) {
1574 for( size_t i=0; i!=ma._nodes.size(); ++i)
1575 os << ( i ? ", " : "" ) << ma._sortedNodeIDs[i];
1585 //================================================================================
1587 * \brief Return nb of elements in the group
1589 //================================================================================
1591 mcIdType Group::size() const
1593 std::size_t 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();
1601 for ( size_t i = 0; i < _groups.size(); ++i )
1602 sizze += _groups[i]->size();
1603 return ToIdType( sizze );
1606 //================================================================================
1608 * \brief Conver gibi element type to med one
1610 //================================================================================
1612 INTERP_KERNEL::NormalizedCellType SauvUtilities::gibi2medGeom( size_t gibiType )
1614 if ( gibiType < 1 || gibiType > NbGibiCellTypes )
1617 return GibiTypeToMed[ gibiType - 1 ];
1620 //================================================================================
1622 * \brief Conver med element type to gibi one
1624 //================================================================================
1626 int SauvUtilities::med2gibiGeom( INTERP_KERNEL::NormalizedCellType medGeomType )
1628 for ( unsigned int i = 0; i < NbGibiCellTypes; i++ )
1629 if ( GibiTypeToMed[ i ] == medGeomType )
1635 //================================================================================
1637 * \brief Remember the file name
1639 //================================================================================
1641 FileReader::FileReader(const char* fileName):_fileName(fileName),_iRead(0),_nbToRead(0)
1645 //================================================================================
1647 * \brief Constructor of ASCII sauve file reader
1649 //================================================================================
1651 ASCIIReader::ASCIIReader(const char* fileName)
1652 :FileReader(fileName),
1657 //================================================================================
1659 * \brief Return true
1661 //================================================================================
1663 bool ASCIIReader::isASCII() const
1668 //================================================================================
1670 * \brief Try to open an ASCII file
1672 //================================================================================
1674 bool ASCIIReader::open()
1677 _file = ::_open (_fileName.c_str(), _O_RDONLY|_O_BINARY);
1679 _file = ::open (_fileName.c_str(), O_RDONLY);
1683 _start = new char [GIBI_BufferSize]; // working buffer beginning
1684 //_tmpBuf = new char [GIBI_MaxOutputLen];
1691 //THROW_IK_EXCEPTION("Can't open file "<<_fileName << " fd: " << _file);
1693 return (_file >= 0);
1696 //================================================================================
1698 * \brief Close the file
1700 //================================================================================
1702 ASCIIReader::~ASCIIReader()
1710 //delete [] _tmpBuf;
1717 //================================================================================
1719 * \brief Return a next line of the file
1721 //================================================================================
1723 bool ASCIIReader::getNextLine (char* & line, bool raiseOEF /*= true*/ )
1725 if ( getLine( line )) return true;
1727 THROW_IK_EXCEPTION("Unexpected EOF on ln "<<_lineNb);
1731 //================================================================================
1733 * \brief Read a next line of the file if necessary
1735 //================================================================================
1737 bool ASCIIReader::getLine(char* & line)
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 std::size_t nBytesRest = _eptr - _ptr;
1743 if (nBytesRest < GIBI_MaxOutputLen)
1747 // move the remaining portion to the buffer beginning
1748 for ( std::size_t i = 0; i < nBytesRest; ++i )
1749 _start[i] = _ptr[i];
1750 //memcpy (_tmpBuf, _ptr, nBytesRest);
1751 //memcpy (_start, _tmpBuf, nBytesRest);
1758 const std::size_t nBytesRead = ::read (_file,
1759 &_start [nBytesRest],
1760 GIBI_BufferSize - nBytesRest);
1761 nBytesRest += nBytesRead;
1762 _eptr = &_start [nBytesRest];
1764 // Check the buffer for the end-of-line
1768 // Check for end-of-the-buffer, the ultimate criterion for termination
1771 if (nBytesRest <= 0)
1777 // seek the line-feed character
1780 if (ptr[-1] == '\r')
1788 // Output the result
1796 //================================================================================
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
1804 //================================================================================
1806 void ASCIIReader::init( int nbToRead, int nbPosInLine, int width, int shift /*= 0*/ )
1808 _nbToRead = nbToRead;
1809 _nbPosInLine = nbPosInLine;
1815 getNextLine( _curPos );
1816 _curPos = _curPos + _shift;
1824 //================================================================================
1826 * \brief Prepare for iterating over given nb of string values
1828 //================================================================================
1830 void ASCIIReader::initNameReading(int nbValues, int width /*= 8*/)
1832 init( nbValues, 72 / ( width + 1 ), width, 1 );
1835 //================================================================================
1837 * \brief Prepare for iterating over given nb of integer values
1839 //================================================================================
1841 void ASCIIReader::initIntReading(int nbValues)
1843 init( nbValues, 10, 8 );
1846 //================================================================================
1848 * \brief Prepare for iterating over given nb of real values
1850 //================================================================================
1852 void ASCIIReader::initDoubleReading(int nbValues)
1854 init( nbValues, 3, 22 );
1857 //================================================================================
1859 * \brief Return true if not all values have been read
1861 //================================================================================
1863 bool ASCIIReader::more() const
1865 bool result = false;
1866 if ( _iRead < _nbToRead)
1868 if ( _curPos ) result = true;
1873 //================================================================================
1875 * \brief Go to the nex value
1877 //================================================================================
1879 void ASCIIReader::next()
1882 THROW_IK_EXCEPTION("SauvUtilities::ASCIIReader::next(): no more() values to read");
1885 if ( _iRead < _nbToRead )
1887 if ( _iPos >= _nbPosInLine )
1889 getNextLine( _curPos );
1890 _curPos = _curPos + _shift;
1895 _curPos = _curPos + _width + _shift;
1904 //================================================================================
1906 * \brief Return the current integer value
1908 //================================================================================
1910 int ASCIIReader::getInt() const
1912 // fix for two glued ints (issue 0021009):
1913 // Line nb | File contents
1914 // ------------------------------------------------------------------------------------
1915 // 53619905 | 1 2 6 8
1916 // 53619906 | SCALAIRE
1917 // 53619907 | -63312600499 1 0 0 0 -2 0 2
1918 // where -63312600499 is actually -633 and 12600499
1919 char hold=_curPos[_width];
1920 _curPos[_width] = '\0';
1921 int result = atoi( _curPos );
1922 _curPos[_width] = hold;
1924 //return atoi(str());
1927 //================================================================================
1929 * \brief Return the current float value
1931 //================================================================================
1933 float ASCIIReader::getFloat() const
1935 return (float)getDouble();
1938 //================================================================================
1940 * \brief Return the current double value
1942 //================================================================================
1944 double ASCIIReader::getDouble() const
1946 //std::string aStr (_curPos);
1948 // Correction: add missing 'E' specifier
1949 // int aPosStart = aStr.find_first_not_of(" \t");
1950 // if (aPosStart < (int)aStr.length()) {
1951 // int aPosSign = aStr.find_first_of("+-", aPosStart + 1); // pass the leading symbol, as it can be a sign
1952 // if (aPosSign < (int)aStr.length()) {
1953 // if (aStr[aPosSign - 1] != 'e' && aStr[aPosSign - 1] != 'E')
1954 // aStr.insert(aPosSign, "E", 1);
1958 // Different Correction (more optimal)
1960 // 0.00000000000000E+00 -2.37822406690632E+01 6.03062748797469E+01
1961 // 7.70000000000000-100 7.70000000000000+100 7.70000000000000+100
1962 //0123456789012345678901234567890123456789012345678901234567890123456789
1963 const size_t posE = 18;
1964 std::string aStr (_curPos);
1965 if ( aStr.find('E') == std::string::npos && aStr.find('e') == std::string::npos )
1967 if ( aStr.size() < posE+1 )
1968 THROW_IK_EXCEPTION("No more doubles (line #" << lineNb() << ")");
1969 aStr.insert( posE, "E", 1 );
1970 return atof(aStr.c_str());
1972 return atof( _curPos );
1975 //================================================================================
1977 * \brief Return the current string value
1979 //================================================================================
1981 std::string ASCIIReader::getName() const
1984 while (( _curPos[len-1] == ' ' || _curPos[len-1] == 0) && len > 0 )
1986 return std::string( _curPos, len );
1989 //================================================================================
1991 * \brief Constructor of a binary sauve file reader
1993 //================================================================================
1995 XDRReader::XDRReader(const char* fileName) :FileReader(fileName), _xdrs_file(NULL)
1999 //================================================================================
2001 * \brief Close the XDR sauve file
2003 //================================================================================
2005 XDRReader::~XDRReader()
2010 xdr_destroy((XDR*)_xdrs);
2012 ::fclose(_xdrs_file);
2018 //================================================================================
2020 * \brief Return false
2022 //================================================================================
2024 bool XDRReader::isASCII() const
2029 //================================================================================
2031 * \brief Try to open an XRD file
2033 //================================================================================
2035 bool XDRReader::open()
2037 bool xdr_ok = false;
2040 if ((_xdrs_file = ::fopen(_fileName.c_str(), "rb")))
2042 if ((_xdrs_file = ::fopen(_fileName.c_str(), "r")))
2045 _xdrs = (XDR *)malloc(sizeof(XDR));
2046 xdrstdio_create((XDR*)_xdrs, _xdrs_file, XDR_DECODE);
2048 const int maxsize = 10;
2049 char icha[maxsize+1];
2051 if (( xdr_ok = xdr_string((XDR*)_xdrs, &icha2, maxsize)))
2053 icha[maxsize] = '\0';
2054 xdr_ok = (strcmp(icha, "CASTEM XDR") == 0);
2058 xdr_destroy((XDR*)_xdrs);
2068 //================================================================================
2072 //================================================================================
2074 bool XDRReader::getNextLine (char* &, bool )
2079 //================================================================================
2081 * \brief Prepare for iterating over given nb of values
2082 * \param nbToRead - nb of fields to read
2083 * \param width - field width
2085 //================================================================================
2087 void XDRReader::init( int nbToRead, int width/*=0*/ )
2089 if(_iRead < _nbToRead)
2091 std::cout << "_iRead, _nbToRead : " << _iRead << " " << _nbToRead << std::endl;
2092 std::cout << "Unfinished iteration before new one !" << std::endl;
2093 THROW_IK_EXCEPTION("SauvUtilities::XDRReader::init(): Unfinished iteration before new one !");
2096 _nbToRead = nbToRead;
2100 //================================================================================
2102 * \brief Prepare for iterating over given nb of string values
2104 //================================================================================
2106 void XDRReader::initNameReading(int nbValues, int width)
2108 init( nbValues, width );
2109 _xdr_kind = _xdr_kind_char;
2112 unsigned int nels = nbValues*width;
2113 _xdr_cvals = (char*)malloc((nels+1)*sizeof(char));
2115 xdr_string((XDR*)_xdrs, &_xdr_cvals, nels);
2117 _xdr_cvals[nels] = '\0';
2121 //================================================================================
2123 * \brief Prepare for iterating over given nb of integer values
2125 //================================================================================
2127 void XDRReader::initIntReading(int nbValues)
2130 _xdr_kind = _xdr_kind_int;
2134 unsigned int nels = nbValues;
2135 unsigned int actual_nels;
2136 _xdr_ivals = (int*)malloc(nels*sizeof(int));
2137 xdr_array((XDR*)_xdrs, (char **)&_xdr_ivals, &actual_nels, nels, sizeof(int), (xdrproc_t)xdr_int);
2142 //================================================================================
2144 * \brief Prepare for iterating over given nb of real values
2146 //================================================================================
2148 void XDRReader::initDoubleReading(int nbValues)
2151 _xdr_kind = _xdr_kind_double;
2155 unsigned int nels = nbValues;
2156 unsigned int actual_nels;
2157 _xdr_dvals = (double*)malloc(nels*sizeof(double));
2158 xdr_array((XDR*)_xdrs, (char **)&_xdr_dvals, &actual_nels, nels, sizeof(double), (xdrproc_t)xdr_double);
2163 //================================================================================
2165 * \brief Return true if not all values have been read
2167 //================================================================================
2169 bool XDRReader::more() const
2171 return _iRead < _nbToRead;
2174 //================================================================================
2176 * \brief Go to the nex value
2178 //================================================================================
2180 void XDRReader::next()
2183 THROW_IK_EXCEPTION("SauvUtilities::XDRReader::next(): no more() values to read");
2186 if ( _iRead < _nbToRead )
2191 if(_xdr_kind == _xdr_kind_char) free(_xdr_cvals);
2192 if(_xdr_kind == _xdr_kind_int) free(_xdr_ivals);
2193 if(_xdr_kind == _xdr_kind_double) free(_xdr_dvals);
2194 _xdr_kind = _xdr_kind_null;
2198 //================================================================================
2200 * \brief Return the current integer value
2202 //================================================================================
2204 int XDRReader::getInt() const
2206 if(_iRead < _nbToRead)
2208 return _xdr_ivals[_iRead];
2214 xdr_int((XDR*)_xdrs, &result);
2220 //================================================================================
2222 * \brief Return the current float value
2224 //================================================================================
2226 float XDRReader::getFloat() const
2230 xdr_float((XDR*)_xdrs, &result);
2235 //================================================================================
2237 * \brief Return the current double value
2239 //================================================================================
2241 double XDRReader::getDouble() const
2243 if(_iRead < _nbToRead)
2245 return _xdr_dvals[_iRead];
2251 xdr_double((XDR*)_xdrs, &result);
2257 //================================================================================
2259 * \brief Return the current string value
2261 //================================================================================
2263 std::string XDRReader::getName() const
2266 char* s = _xdr_cvals + _iRead*_width;
2267 while (( s[len-1] == ' ' || s[len-1] == 0) && len > 0 )
2269 return std::string( s, len );
2272 //================================================================================
2274 * \brief Throw an exception if not all needed data is present
2276 //================================================================================
2278 void IntermediateMED::checkDataAvailability() const
2280 if ( _spaceDim == 0 )
2281 THROW_IK_EXCEPTION("Wrong file format"); // it is the first record in the sauve file
2283 if ( _groups.empty() )
2284 THROW_IK_EXCEPTION("No elements have been read");
2286 if ( _points.empty() || _nbNodes == 0 )
2287 THROW_IK_EXCEPTION("Nodes of elements are not filled");
2289 if ( _coords.empty() )
2290 THROW_IK_EXCEPTION("Node coordinates are missing");
2292 if ( _coords.size() < _nbNodes * _spaceDim )
2293 THROW_IK_EXCEPTION("Nodes and coordinates mismatch");
2296 //================================================================================
2298 * \brief Safely adds a new Group
2300 //================================================================================
2302 Group* IntermediateMED::addNewGroup(std::vector<SauvUtilities::Group*>* groupsToFix)
2304 if ( _groups.size() == _groups.capacity() ) // re-allocation would occur
2306 std::vector<Group> newGroups( _groups.size() );
2307 newGroups.push_back( Group() );
2309 for ( size_t i = 0; i < _groups.size(); ++i )
2311 // avoid copying _cells
2312 std::vector<const Cell*> cells;
2313 cells.swap( _groups[i]._cells );
2314 newGroups[i] = _groups[i];
2315 newGroups[i]._cells.swap( cells );
2317 // correct pointers to sub-groups
2318 for ( size_t j = 0; j < _groups[i]._groups.size(); ++j )
2320 std::size_t iG = _groups[i]._groups[j] - &_groups[0];
2321 newGroups[i]._groups[j] = & newGroups[ iG ];
2327 for ( size_t i = 0; i < groupsToFix->size(); ++i )
2328 if ( (*groupsToFix)[i] )
2330 std::size_t iG = (*groupsToFix)[i] - &_groups[0];
2331 (*groupsToFix)[i] = & newGroups[ iG ];
2334 // fix field supports
2335 for ( int isNode = 0; isNode < 2; ++isNode )
2337 std::vector<DoubleField* >& fields = isNode ? _nodeFields : _cellFields;
2338 for ( size_t i = 0; i < fields.size(); ++i )
2340 if ( !fields[i] ) continue;
2341 for ( size_t j = 0; j < fields[i]->_sub.size(); ++j )
2342 if ( fields[i]->_sub[j]._support )
2344 std::size_t iG = fields[i]->_sub[j]._support - &_groups[0];
2345 fields[i]->_sub[j]._support = & newGroups[ iG ];
2347 if ( fields[i]->_group )
2349 std::size_t iG = fields[i]->_group - &_groups[0];
2350 fields[i]->_group = & newGroups[ iG ];
2355 _groups.swap( newGroups );
2359 _groups.push_back( Group() );
2361 return &_groups.back();
2364 //================================================================================
2366 * \brief Makes MEDCoupling::MEDFileData from self
2368 //================================================================================
2370 MEDCoupling::MEDFileData* IntermediateMED::convertInMEDFileDS()
2372 MCAuto< MEDFileUMesh > mesh = makeMEDFileMesh();
2373 MCAuto< MEDFileFields > fields = makeMEDFileFields(mesh);
2375 MCAuto< MEDFileMeshes > meshes = MEDFileMeshes::New();
2376 MCAuto< MEDFileData > medData = MEDFileData::New();
2377 meshes->pushMesh( mesh );
2378 medData->setMeshes( meshes );
2379 if ( fields ) medData->setFields( fields );
2381 return medData.retn();
2384 //================================================================================
2386 * \brief Creates MEDCoupling::MEDFileUMesh from its data
2388 //================================================================================
2390 MEDCoupling::MEDFileUMesh* IntermediateMED::makeMEDFileMesh()
2392 // check if all needed piles are present
2393 checkDataAvailability();
2395 decreaseHierarchicalDepthOfSubgroups();
2397 // set long names (before orienting!)
2398 setGroupLongNames();
2400 // fix element orientation
2401 if ( _spaceDim == 2 || _spaceDim == 1 )
2403 else if ( _spaceDim == 3 )
2407 eraseUselessGroups();
2408 //detectMixDimGroups();
2411 _points.numberNodes();
2414 // make the med mesh
2416 MEDFileUMesh* mesh = MEDFileUMesh::New();
2418 DataArrayDouble *coords = getCoords();
2419 setConnectivity( mesh, coords );
2424 if ( !mesh->getName().c_str() || strlen( mesh->getName().c_str() ) == 0 )
2425 mesh->setName( "MESH" );
2430 //================================================================================
2432 * \brief Set long names to groups
2434 //================================================================================
2436 void IntermediateMED::setGroupLongNames()
2438 if ( _listGIBItoMED_mail.empty() )
2441 // IMP 0023285: only keep the meshes named in the table MED_MAIL
2442 // clear all group names
2443 for ( size_t i = 0; i < _groups.size(); ++i )
2444 if ( !_groups[i]._isProfile )
2445 _groups[i]._name.clear();
2448 // IMP 0020434: mapping GIBI names to MED names
2449 // set med names to objects (mesh, fields, support, group or other)
2451 std::set<int> treatedGroups;
2453 std::list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_mail.begin();
2454 for (; itGIBItoMED != _listGIBItoMED_mail.end(); itGIBItoMED++)
2456 if ( (int)_groups.size() < itGIBItoMED->gibi_id ) continue;
2458 SauvUtilities::Group & grp = _groups[itGIBItoMED->gibi_id - 1];
2460 // if there are several names for grp then the 1st name is the name
2461 // of grp and the rest ones are names of groups referring grp (issue 0021311)
2462 const bool isRefName = !treatedGroups.insert( itGIBItoMED->gibi_id ).second;
2465 grp._name = _mapStrings[ itGIBItoMED->med_id ];
2467 else if ( !grp._refNames.empty() && grp._refNames.back().empty() )
2469 for ( unsigned i = 0; i < grp._refNames.size(); ++i )
2470 if ( grp._refNames[i].empty() )
2471 grp._refNames[i] = _mapStrings[ (*itGIBItoMED).med_id ];
2475 grp._refNames.push_back( _mapStrings[ (*itGIBItoMED).med_id ]);
2479 // IMP 0023285: only keep the meshes named in the table MED_MAIL
2480 // remove all cells belonging to non-named groups only
2482 // use Cell::_reverse to mark cells to keep
2483 for ( size_t i = 0; i < _groups.size(); ++i )
2485 SauvUtilities::Group & grp = _groups[i];
2486 if ( grp._isProfile || !grp._name.empty() )
2488 for ( size_t iC = 0; iC < grp._cells.size(); ++iC )
2489 grp._cells[iC]->_reverse = true;
2491 for (size_t j = 0; j < grp._groups.size(); ++j )
2492 for ( size_t iC = 0; iC < grp._groups[j]->_cells.size(); ++iC )
2493 grp._groups[j]->_cells[iC]->_reverse = true;
2496 // remove non-marked cells (with _reverse == false)
2497 CellsByDimIterator cellsIt( *this );
2498 while ( cellsIt.nextType() )
2500 std::set<Cell> & cells = _cellsByType[ cellsIt.type() ];
2501 std::set<Cell>::iterator cIt = cells.begin();
2502 while ( cIt != cells.end() )
2503 if ( cIt->_reverse )
2505 cIt->_reverse = false;
2510 cells.erase( cIt++ );
2515 //================================================================================
2517 * \brief Set long names to fields
2519 //================================================================================
2521 void IntermediateMED::setFieldLongNames(std::set< std::string >& usedNames)
2523 std::list<nameGIBItoMED>::iterator itGIBItoMED = _listGIBItoMED_cham.begin();
2524 for (; itGIBItoMED != _listGIBItoMED_cham.end(); itGIBItoMED++)
2526 if (itGIBItoMED->gibi_pile == PILE_FIELD)
2528 _cellFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2530 else if (itGIBItoMED->gibi_pile == PILE_NODES_FIELD)
2532 _nodeFields[itGIBItoMED->gibi_id - 1]->_name = _mapStrings[itGIBItoMED->med_id];
2534 } // iterate on _listGIBItoMED_cham
2536 for (itGIBItoMED =_listGIBItoMED_comp.begin(); itGIBItoMED != _listGIBItoMED_comp.end(); itGIBItoMED++)
2538 std::string medName = _mapStrings[itGIBItoMED->med_id];
2539 std::string gibiName = _mapStrings[itGIBItoMED->gibi_id];
2541 bool name_found = false;
2542 for ( int isNodal = 0; isNodal < 2 && !name_found; ++isNodal )
2544 std::vector<DoubleField* > & fields = isNodal ? _nodeFields : _cellFields;
2545 for ( size_t ifi = 0; ifi < fields.size() && !name_found; ifi++)
2547 if (medName.find( fields[ifi]->_name + "." ) == 0 )
2549 std::vector<DoubleField::_Sub_data>& aSubDs = fields[ifi]->_sub;
2550 std::size_t nbSub = aSubDs.size();
2551 for (std::size_t isu = 0; isu < nbSub; isu++)
2552 for (int ico = 0; ico < aSubDs[isu].nbComponents(); ico++)
2554 if (aSubDs[isu].compName(ico) == gibiName)
2556 std::string medNameCompo = medName.substr( fields[ifi]->_name.size() + 1 );
2557 fields[ifi]->_sub[isu].compName(ico) = medNameCompo;
2563 } // iterate on _listGIBItoMED_comp
2565 for ( size_t i = 0; i < _nodeFields.size() ; i++)
2566 usedNames.insert( _nodeFields[i]->_name );
2567 for ( size_t i = 0; i < _cellFields.size() ; i++)
2568 usedNames.insert( _cellFields[i]->_name );
2571 //================================================================================
2573 * \brief Decrease hierarchical depth of subgroups
2575 //================================================================================
2577 void IntermediateMED::decreaseHierarchicalDepthOfSubgroups()
2579 for (size_t i=0; i!=_groups.size(); ++i)
2581 Group& grp = _groups[i];
2582 for (size_t j = 0; j < grp._groups.size(); ++j )
2584 Group & sub_grp = *grp._groups[j];
2585 if ( !sub_grp._groups.empty() )
2587 // replace j with its 1st subgroup
2588 grp._groups[j] = sub_grp._groups[0];
2589 // push back the rest subs
2590 grp._groups.insert( grp._groups.end(), ++sub_grp._groups.begin(), sub_grp._groups.end() );
2593 // remove empty sub-_groups
2594 std::vector< Group* > newSubGroups;
2595 newSubGroups.reserve( grp._groups.size() );
2596 for (size_t j = 0; j < grp._groups.size(); ++j )
2597 if ( !grp._groups[j]->empty() )
2598 newSubGroups.push_back( grp._groups[j] );
2599 if ( newSubGroups.size() < grp._groups.size() )
2600 grp._groups.swap( newSubGroups );
2604 //================================================================================
2606 * \brief Erase _groups that won't be converted
2608 //================================================================================
2610 void IntermediateMED::eraseUselessGroups()
2612 // propagate _isProfile=true to sub-groups of composite groups
2613 // for (size_t int i=0; i!=_groups.size(); ++i)
2615 // Group* grp = _groups[i];
2616 // if ( grp->_isProfile && !grp->_groups.empty() )
2617 // for (size_t j = 0; j < grp->_groups.size(); ++j )
2618 // grp->_groups[j]->_isProfile=true;
2620 std::set<Group*> groups2convert;
2621 // keep not named sub-groups of field supports
2622 for (size_t i=0; i!=_groups.size(); ++i)
2624 Group& grp = _groups[i];
2625 if ( grp._isProfile && !grp._groups.empty() )
2626 groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2629 // keep named groups and their subgroups
2630 for (size_t i=0; i!=_groups.size(); ++i)
2632 Group& grp = _groups[i];
2633 if ( !grp._name.empty() && !grp.empty() )
2635 groups2convert.insert( &grp );
2636 groups2convert.insert( grp._groups.begin(), grp._groups.end() );
2639 // erase groups that are not in groups2convert and not _isProfile
2640 for (size_t i=0; i!=_groups.size(); ++i)
2642 Group* grp = &_groups[i];
2643 if ( !grp->_isProfile && !groups2convert.count( grp ) )
2645 grp->_cells.clear();
2646 grp->_groups.clear();
2651 //================================================================================
2653 * \brief Detect _groups of mixed dimension
2655 //================================================================================
2657 void IntermediateMED::detectMixDimGroups()
2659 //hasMixedCells = false;
2660 for ( size_t i=0; i < _groups.size(); ++i )
2662 Group& grp = _groups[i];
2663 if ( grp._groups.size() < 2 )
2666 // check if sub-groups have different dimension
2667 unsigned dim1 = getDim( &grp );
2668 for ( size_t j = 1; j < grp._groups.size(); ++j )
2670 unsigned dim2 = getDim( grp._groups[j] );
2674 grp._groups.clear();
2675 if ( !grp._name.empty() )
2676 std::cout << "Erase a group with elements of different dim |" << grp._name << "|"<< std::endl;
2683 //================================================================================
2685 * \brief Fix connectivity of elements in 2D space
2687 //================================================================================
2689 void IntermediateMED::orientElements2D()
2691 std::set<Cell>::const_iterator elemIt, elemEnd;
2692 std::vector< std::pair<int,int> > swapVec;
2694 // ------------------------------------
2695 // fix connectivity of quadratic edges
2696 // ------------------------------------
2697 std::set<Cell>& quadEdges = _cellsByType[ INTERP_KERNEL::NORM_SEG3 ];
2698 if ( !quadEdges.empty() )
2700 elemIt = quadEdges.begin(), elemEnd = quadEdges.end();
2701 for ( ; elemIt != elemEnd; ++elemIt )
2702 ConvertQuadratic( INTERP_KERNEL::NORM_SEG3, *elemIt );
2705 CellsByDimIterator faceIt( *this, 2 );
2706 while ( const std::set<Cell > * faces = faceIt.nextType() )
2708 TCellType cellType = faceIt.type();
2709 bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2711 getReverseVector( cellType, swapVec );
2713 // ------------------------------------
2714 // fix connectivity of quadratic faces
2715 // ------------------------------------
2717 for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2718 ConvertQuadratic( cellType, *elemIt );
2720 // --------------------------
2721 // orient faces clockwise
2722 // --------------------------
2723 // COMMENTED for issue 0022612 note 17739
2724 // int iQuad = isQuadratic ? 2 : 1;
2725 // for ( elemIt = faces->begin(), elemEnd = faces->end(); elemIt != elemEnd; elemIt++ )
2727 // // look for index of the most left node
2728 // int iLeft = 0, iNode, nbNodes = elemIt->_nodes.size() / iQuad;
2729 // double x, minX = nodeCoords( elemIt->_nodes[0] )[0];
2730 // for ( iNode = 1; iNode < nbNodes; ++iNode )
2731 // if (( x = nodeCoords( elemIt->_nodes[ iNode ])[ 0 ]) < minX )
2732 // minX = x, iLeft = iNode;
2734 // // indices of the nodes neighboring the most left one
2735 // int iPrev = ( iLeft - 1 < 0 ) ? nbNodes - 1 : iLeft - 1;
2736 // int iNext = ( iLeft + 1 == nbNodes ) ? 0 : iLeft + 1;
2737 // // find components of prev-left and left-next vectors
2738 // double xP = nodeCoords( elemIt->_nodes[ iPrev ])[ 0 ];
2739 // double yP = nodeCoords( elemIt->_nodes[ iPrev ])[ 1 ];
2740 // double xN = nodeCoords( elemIt->_nodes[ iNext ])[ 0 ];
2741 // double yN = nodeCoords( elemIt->_nodes[ iNext ])[ 1 ];
2742 // double xL = nodeCoords( elemIt->_nodes[ iLeft ])[ 0 ];
2743 // double yL = nodeCoords( elemIt->_nodes[ iLeft ])[ 1 ];
2744 // double xPL = xL - xP, yPL = yL - yP; // components of prev-left vector
2745 // double xLN = xN - xL, yLN = yN - yL; // components of left-next vector
2746 // // normalise y of the vectors
2747 // double modPL = sqrt ( xPL * xPL + yPL * yPL );
2748 // double modLN = sqrt ( xLN * xLN + yLN * yLN );
2749 // if ( modLN > std::numeric_limits<double>::min() &&
2750 // modPL > std::numeric_limits<double>::min() )
2754 // // summary direction of neighboring links must be positive
2755 // bool clockwise = ( yPL + yLN > 0 );
2756 // if ( !clockwise )
2757 // reverse( *elemIt, swapVec );
2763 //================================================================================
2765 * \brief Fix connectivity of elements in 3D space
2767 //================================================================================
2769 void IntermediateMED::orientElements3D()
2771 // set _reverse flags of faces
2772 // COMMENTED for issue 0022612 note 17739
2775 // -----------------
2777 // -----------------
2779 std::set<Cell>::const_iterator elemIt, elemEnd;
2780 std::vector< std::pair<int,int> > swapVec;
2782 for ( int dim = 1; dim <= 3; ++dim )
2784 CellsByDimIterator cellsIt( *this, dim );
2785 while ( const std::set<Cell > * elems = cellsIt.nextType() )
2787 TCellType cellType = cellsIt.type();
2788 bool isQuadratic = getGibi2MedQuadraticInterlace( cellType );
2789 getReverseVector( cellType, swapVec );
2791 elemIt = elems->begin(), elemEnd = elems->end();
2792 for ( ; elemIt != elemEnd; elemIt++ )
2794 // GIBI connectivity -> MED one
2796 ConvertQuadratic( cellType, *elemIt );
2799 if ( elemIt->_reverse )
2800 reverse ( *elemIt, swapVec );
2805 // COMMENTED for issue 0022612 note 17739
2809 //================================================================================
2811 * \brief Orient equally (by setting _reverse flag) all connected faces in 3D space
2813 //================================================================================
2815 void IntermediateMED::orientFaces3D()
2817 // fill map of links and their faces
2818 std::set<const Cell*> faces;
2819 std::map<const Cell*, Group*> fgm;
2820 std::map<Link, std::list<const Cell*> > linkFacesMap;
2821 std::map<Link, std::list<const Cell*> >::iterator lfIt, lfIt2;
2823 for (size_t i=0; i!=_groups.size(); ++i)
2825 Group& grp = _groups[i];
2826 if ( !grp._cells.empty() && getDimension( grp._cellType ) == 2 )
2827 for ( size_t j = 0; j < grp._cells.size(); ++j )
2828 if ( faces.insert( grp._cells[j] ).second )
2830 for ( unsigned int k = 0; k < grp._cells[j]->_nodes.size(); ++k )
2831 linkFacesMap[ grp._cells[j]->link( k ) ].push_back( grp._cells[j] );
2832 fgm.insert( std::make_pair( grp._cells[j], &grp ));
2835 // dump linkFacesMap
2836 // for ( lfIt = linkFacesMap.begin(); lfIt!=linkFacesMap.end(); lfIt++) {
2837 // cout<< "LINK: " << lfIt->first.first << "-" << lfIt->first.second << std::endl;
2838 // std::list<const Cell*> & fList = lfIt->second;
2839 // std::list<const Cell*>::iterator fIt = fList.begin();
2840 // for ( ; fIt != fList.end(); fIt++ )
2841 // cout << "\t" << **fIt << fgm[*fIt]->nom << std::endl;
2844 // Each oriented link must appear in one face only, else a face is reversed.
2846 std::queue<const Cell*> faceQueue; /* the queue contains well oriented faces
2847 whose neighbors orientation is to be checked */
2848 bool manifold = true;
2849 while ( !linkFacesMap.empty() )
2851 if ( faceQueue.empty() )
2853 assert( !linkFacesMap.begin()->second.empty() );
2854 faceQueue.push( linkFacesMap.begin()->second.front() );
2856 while ( !faceQueue.empty() )
2858 const Cell* face = faceQueue.front();
2861 // loop on links of <face>
2862 for ( int i = 0; i < (int)face->_nodes.size(); ++i )
2864 Link link = face->link( i );
2865 // find the neighbor faces
2866 lfIt = linkFacesMap.find( link );
2867 int nbFaceByLink = 0;
2868 std::list< const Cell* > ml;
2869 if ( lfIt != linkFacesMap.end() )
2871 std::list<const Cell*> & fList = lfIt->second;
2872 std::list<const Cell*>::iterator fIt = fList.begin();
2873 assert( fIt != fList.end() );
2874 for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2876 ml.push_back( *fIt );
2877 if ( *fIt != face ) // wrongly oriented neighbor face
2879 const Cell* badFace = *fIt;
2880 // reverse and remove badFace from linkFacesMap
2881 for ( int j = 0; j < (int)badFace->_nodes.size(); ++j )
2883 Link badlink = badFace->link( j );
2884 if ( badlink == link ) continue;
2885 lfIt2 = linkFacesMap.find( badlink );
2886 if ( lfIt2 != linkFacesMap.end() )
2888 std::list<const Cell*> & ff = lfIt2->second;
2889 std::list<const Cell*>::iterator lfIt3 = find( ff.begin(), ff.end(), badFace );
2890 // check if badFace has been found,
2891 // else we can't erase it
2892 // case of degenerated face in edge
2893 if (lfIt3 != ff.end())
2897 linkFacesMap.erase( lfIt2 );
2901 badFace->_reverse = true; // reverse
2902 //INFOS_MED( "REVERSE " << *badFace );
2903 faceQueue.push( badFace );
2906 linkFacesMap.erase( lfIt );
2908 // add good neighbors to the queue
2909 Link revLink( link.second, link.first );
2910 lfIt = linkFacesMap.find( revLink );
2911 if ( lfIt != linkFacesMap.end() )
2913 std::list<const Cell*> & fList = lfIt->second;
2914 std::list<const Cell*>::iterator fIt = fList.begin();
2915 for ( ; fIt != fList.end(); fIt++, nbFaceByLink++ )
2917 ml.push_back( *fIt );
2919 faceQueue.push( *fIt );
2921 linkFacesMap.erase( lfIt );
2923 if ( nbFaceByLink > 2 )
2927 std::list<const Cell*>::iterator ii = ml.begin();
2928 std::cout << nbFaceByLink << " faces by 1 link:" << std::endl;
2929 for( ; ii!= ml.end(); ii++ )
2930 std::cout << "in sub-mesh <" << fgm[ *ii ]->_name << "> " << **ii << std::endl;
2934 } // loop on links of the being checked face
2935 } // loop on the face queue
2936 } // while ( !linkFacesMap.empty() )
2939 std::cout << " -> Non manifold mesh, faces orientation may be incorrect" << std::endl;
2942 //================================================================================
2944 * \brief Orient volumes according to MED conventions:
2945 * normal of a bottom (first) face should be outside
2947 //================================================================================
2949 void IntermediateMED::orientVolumes()
2951 std::set<Cell>::const_iterator elemIt, elemEnd;
2952 std::vector< std::pair<int,int> > swapVec;
2954 CellsByDimIterator cellsIt( *this, 3 );
2955 while ( const std::set<Cell > * elems = cellsIt.nextType() )
2957 TCellType cellType = cellsIt.type();
2958 elemIt = elems->begin(), elemEnd = elems->end();
2959 int nbBottomNodes = 0;
2966 nbBottomNodes = 3; break;
2971 nbBottomNodes = 4; break;
2974 getReverseVector( cellType, swapVec );
2976 for ( ; elemIt != elemEnd; elemIt++ )
2978 // find a normal to the bottom face
2980 n[0] = nodeCoords( elemIt->_nodes[0]); // 3 bottom nodes
2981 n[1] = nodeCoords( elemIt->_nodes[1]);
2982 n[2] = nodeCoords( elemIt->_nodes[2]);
2983 n[3] = nodeCoords( elemIt->_nodes[nbBottomNodes]); // a top node
2984 double vec01[3]; // vector n[0]-n[1]
2985 vec01[0] = n[1][0] - n[0][0];
2986 vec01[1] = n[1][1] - n[0][1];
2987 vec01[2] = n[1][2] - n[0][2];
2988 double vec02 [3]; // vector n[0]-n[2]
2989 vec02[0] = n[2][0] - n[0][0];
2990 vec02[1] = n[2][1] - n[0][1];
2991 vec02[2] = n[2][2] - n[0][2];
2992 double normal [3]; // vec01 ^ vec02
2993 normal[0] = vec01[1] * vec02[2] - vec01[2] * vec02[1];
2994 normal[1] = vec01[2] * vec02[0] - vec01[0] * vec02[2];
2995 normal[2] = vec01[0] * vec02[1] - vec01[1] * vec02[0];
2996 // check if the 102 angle is convex
2997 if ( nbBottomNodes > 3 )
2999 const double* n3 = nodeCoords( elemIt->_nodes[nbBottomNodes-1] );// last bottom node
3000 double vec03 [3]; // vector n[0]-n3
3001 vec03[0] = n3[0] - n[0][0];
3002 vec03[1] = n3[1] - n[0][1];
3003 vec03[2] = n3[2] - n[0][2];
3004 if ( fabs( normal[0]+normal[1]+normal[2] ) <= std::numeric_limits<double>::max() ) // vec01 || vec02
3006 normal[0] = vec01[1] * vec03[2] - vec01[2] * vec03[1]; // vec01 ^ vec03
3007 normal[1] = vec01[2] * vec03[0] - vec01[0] * vec03[2];
3008 normal[2] = vec01[0] * vec03[1] - vec01[1] * vec03[0];
3012 double vec [3]; // normal ^ vec01
3013 vec[0] = normal[1] * vec01[2] - normal[2] * vec01[1];
3014 vec[1] = normal[2] * vec01[0] - normal[0] * vec01[2];
3015 vec[2] = normal[0] * vec01[1] - normal[1] * vec01[0];
3016 double dot2 = vec[0]*vec03[0] + vec[1]*vec03[1] + vec[2]*vec03[2]; // vec*vec03
3017 if ( dot2 < 0 ) // concave -> reverse normal
3025 // direction from top to bottom
3027 tbDir[0] = n[0][0] - n[3][0];
3028 tbDir[1] = n[0][1] - n[3][1];
3029 tbDir[2] = n[0][2] - n[3][2];
3031 // compare 2 directions: normal and top-bottom
3032 double dot = normal[0]*tbDir[0] + normal[1]*tbDir[1] + normal[2]*tbDir[2];
3033 if ( dot < 0. ) // need reverse
3034 reverse( *elemIt, swapVec );
3036 } // loop on volumes of one geometry
3037 } // loop on 3D geometry types
3041 //================================================================================
3043 * \brief Assign new IDs to nodes by skipping not used nodes and return their number
3045 //================================================================================
3047 int NodeContainer::numberNodes()
3050 for ( size_t i = 0; i < _nodes.size(); ++i )
3051 for ( size_t j = 0; j < _nodes[i].size(); ++j )
3052 if ( _nodes[i][j].isUsed() )
3053 _nodes[i][j]._number = id++;
3058 //================================================================================
3060 * \brief Assign new IDs to elements
3062 //================================================================================
3064 void IntermediateMED::numberElements()
3066 std::set<Cell>::const_iterator elemIt, elemEnd;
3068 // numbering _cells of type NORM_POINT1 by node number
3070 const std::set<Cell>& points = _cellsByType[ INTERP_KERNEL::NORM_POINT1 ];
3071 elemIt = points.begin(), elemEnd = points.end();
3072 for ( ; elemIt != elemEnd; ++elemIt )
3073 elemIt->_number = elemIt->_nodes[0]->_number;
3076 // numbering 1D-3D _cells
3077 for ( int dim = 1; dim <= 3; ++dim )
3079 // check if re-numeration is needed (to try to keep elem oreder as in sauve file )
3080 bool ok = true, renumEntity = false;
3081 CellsByDimIterator cellsIt( *this, dim );
3082 mcIdType prevNbElems = 0;
3083 while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3085 TID minNumber = std::numeric_limits<TID>::max(), maxNumber = 0;
3086 for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3088 if ( elemIt->_number < minNumber ) minNumber = elemIt->_number;
3089 if ( elemIt->_number > maxNumber ) maxNumber = elemIt->_number;
3091 mcIdType typeSize = ToIdType( typeCells->size() );
3092 if ( typeSize != maxNumber - minNumber + 1 )
3094 if ( prevNbElems+1 != (int)minNumber )
3096 if ( prevNbElems != 0 && minNumber == 1 )
3099 prevNbElems += typeSize;
3102 if ( ok && renumEntity ) // each geom type was numerated separately
3104 cellsIt.init( dim );
3105 prevNbElems = ToIdType( cellsIt.nextType()->size()); // no need to renumber the first type
3106 while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3108 for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3109 elemIt->_number += prevNbElems;
3110 prevNbElems += ToIdType( typeCells->size() );
3116 cellsIt.init( dim );
3117 while ( const std::set<Cell> * typeCells = cellsIt.nextType() )
3118 for ( elemIt = typeCells->begin(), elemEnd = typeCells->end(); elemIt!=elemEnd; ++elemIt)
3119 elemIt->_number = cellID++;
3124 //================================================================================
3126 * \brief Creates coord array
3128 //================================================================================
3130 MEDCoupling::DataArrayDouble * IntermediateMED::getCoords()
3132 DataArrayDouble* coordArray = DataArrayDouble::New();
3133 coordArray->alloc( _nbNodes, _spaceDim );
3134 double * coordPrt = coordArray->getPointer();
3135 for ( unsigned int i = 0; i < _points.size(); ++i )
3137 Node* n = getNode( i+1 );
3140 const double* nCoords = nodeCoords( n );
3141 std::copy( nCoords, nCoords+_spaceDim, coordPrt );
3142 coordPrt += _spaceDim;
3148 //================================================================================
3150 * \brief Sets connectivity of elements to the mesh
3151 * \param mesh - mesh to fill in
3152 * \param coords - coordinates that must be shared by all meshes of different dim
3154 //================================================================================
3156 void IntermediateMED::setConnectivity( MEDCoupling::MEDFileUMesh* mesh,
3157 MEDCoupling::DataArrayDouble* coords )
3161 mesh->setCoords( coords );
3163 std::set<Cell>::const_iterator elemIt, elemEnd;
3164 for ( int dim = 3; dim > 0; --dim )
3166 CellsByDimIterator dimCells( *this, dim );
3168 mcIdType nbOfCells = 0;
3169 while ( const std::set<Cell > * cells = dimCells.nextType() )
3170 nbOfCells += ToIdType( cells->size() );
3171 if ( nbOfCells == 0 )
3174 if ( !meshDim ) meshDim = dim;
3176 MEDCouplingUMesh* dimMesh = MEDCouplingUMesh::New();
3177 dimMesh->setCoords( coords );
3178 dimMesh->setMeshDimension( dim );
3179 dimMesh->allocateCells( nbOfCells );
3181 mcIdType prevNbCells = 0;
3182 dimCells.init( dim );
3183 while ( const std::set<Cell > * cells = dimCells.nextType() )
3185 // fill connectivity array to take into account order of elements in the sauv file
3186 const mcIdType nbCellNodes = ToIdType( cells->begin()->_nodes.size() );
3187 std::vector< TID > connectivity( cells->size() * nbCellNodes );
3188 TID * nodalConnOfCell;
3189 for ( elemIt = cells->begin(), elemEnd = cells->end(); elemIt != elemEnd; ++elemIt )
3191 const Cell& cell = *elemIt;
3192 const TID index = cell._number - 1 - prevNbCells;
3193 nodalConnOfCell = &connectivity[ index * nbCellNodes ];
3194 if ( cell._reverse )
3195 for ( mcIdType i = nbCellNodes-1; i >= 0; --i )
3196 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3198 for ( mcIdType i = 0; i < nbCellNodes; ++i )
3199 *nodalConnOfCell++ = cell._nodes[i]->_number - 1;
3201 prevNbCells += ToIdType( cells->size() );
3204 TCellType cellType = dimCells.type();
3205 nodalConnOfCell = &connectivity[0];
3206 for ( size_t i = 0; i < cells->size(); ++i, nodalConnOfCell += nbCellNodes )
3207 dimMesh->insertNextCell( cellType, nbCellNodes, nodalConnOfCell );
3209 dimMesh->finishInsertingCells();
3210 mesh->setMeshAtLevel( dim - meshDim, dimMesh );
3215 //================================================================================
3217 * \brief Fill in the mesh with groups
3218 * \param mesh - mesh to fill in
3220 //================================================================================
3222 void IntermediateMED::setGroups( MEDCoupling::MEDFileUMesh* mesh )
3224 bool isMeshNameSet = false;
3225 const int meshDim = mesh->getMeshDimension();
3226 for ( int dim = 0; dim <= meshDim; ++dim )
3228 const int meshDimRelToMaxExt = ( dim == 0 ? 1 : dim - meshDim );
3230 std::vector<const DataArrayIdType *> medGroups;
3231 std::vector<MCAuto<DataArrayIdType> > refGroups;
3232 for ( size_t i = 0; i < _groups.size(); ++i )
3234 Group& grp = _groups[i];
3235 if ( (int)getDim( &grp ) != dim &&
3236 grp._groups.empty() ) // to allow groups on diff dims
3238 // convert only named groups or field supports
3239 if ( grp.empty() || (grp._name.empty() && !grp._isProfile ))
3241 //if ( grp._medGroup ) continue; // already converted
3243 // sort cells by ID and remember their initial order in the group
3244 TCellToOrderMap cell2order;
3245 unsigned orderInGroup = 0;
3246 std::vector< Group* > groupVec;
3247 if ( grp._groups.empty() ) groupVec.push_back( & grp );
3248 else groupVec = grp._groups;
3249 for ( size_t iG = 0; iG < groupVec.size(); ++iG )
3251 Group* aG = groupVec[ iG ];
3252 if ( (int)getDim( aG ) != dim )
3254 for ( size_t iC = 0; iC < aG->_cells.size(); ++iC )
3255 cell2order.insert( cell2order.end(), std::make_pair( aG->_cells[iC], orderInGroup++ ));
3257 if ( cell2order.empty() )
3259 bool isSelfIntersect = ( orderInGroup != cell2order.size() );
3260 if ( isSelfIntersect ) // self intersecting group
3262 std::ostringstream msg;
3263 msg << "Self intersecting sub-mesh: id = " << i+1
3264 << ", name = |" << grp._name << "|" << std::endl
3265 << " nb unique elements = " << cell2order.size() << std::endl
3266 << " total nb elements = " << orderInGroup;
3267 if ( grp._isProfile )
3269 THROW_IK_EXCEPTION( msg.str() );
3273 std::cout << msg.str() << std::endl;
3276 // create a med group
3277 grp._medGroup = DataArrayIdType::New();
3278 grp._medGroup->setName( grp._name.c_str() );
3279 grp._medGroup->alloc( cell2order.size(), /*nbOfCompo=*/1 );
3280 TID * idsPtr = grp._medGroup->getPointer();
3281 TCellToOrderMap::iterator cell2orderIt, cell2orderEnd = cell2order.end();
3282 for ( cell2orderIt = cell2order.begin(); cell2orderIt != cell2orderEnd; ++cell2orderIt )
3283 *idsPtr++ = (*cell2orderIt).first->_number - 1;
3285 // try to set the mesh name
3286 if ( !isMeshNameSet &&
3288 !grp._name.empty() &&
3289 grp.size() == mesh->getSizeAtLevel( meshDimRelToMaxExt ))
3291 mesh->setName( grp._name.c_str() );
3292 isMeshNameSet = true;
3294 if ( !grp._name.empty() )
3296 medGroups.push_back( grp._medGroup );
3298 // set relocation table
3299 setRelocationTable( &grp, cell2order );
3301 // Issue 0021311. Use case: a gibi group has references (recorded in pile 1)
3302 // and several names (pile 27) refer (pile 10) to this group.
3303 // We create a copy of this group per each named reference
3304 std::set<std::string> uniqueNames;
3305 uniqueNames.insert( grp._name );
3306 for ( unsigned iRef = 0 ; iRef < grp._refNames.size(); ++iRef )
3307 if ( !grp._refNames[ iRef ].empty() &&
3308 uniqueNames.insert( grp._refNames[ iRef ]).second ) // for name uniqueness (23155)
3310 refGroups.push_back( grp._medGroup->deepCopy() );
3311 refGroups.back()->setName( grp._refNames[ iRef ].c_str() );
3312 medGroups.push_back( refGroups.back() );
3315 mesh->setGroupsAtLevel( meshDimRelToMaxExt, medGroups );
3319 //================================================================================
3321 * \brief Return true if the group is on all elements and return its relative dimension
3323 //================================================================================
3325 bool IntermediateMED::isOnAll( const Group* grp, int & dimRel ) const
3327 int dim = getDim( grp );
3329 mcIdType nbElems = 0;
3337 CellsByDimIterator dimCells( *this, dim );
3338 while ( const std::set<Cell > * cells = dimCells.nextType() )
3339 nbElems += ToIdType( cells->size() );
3342 for ( ; meshDim > 0; --meshDim )
3344 dimCells.init( meshDim );
3345 if ( dimCells.nextType() )
3348 dimRel = dim - meshDim;
3351 bool onAll = ( nbElems == grp->size() );
3355 //================================================================================
3357 * \brief Makes fields from own data
3359 //================================================================================
3361 MEDCoupling::MEDFileFields * IntermediateMED::makeMEDFileFields(MEDCoupling::MEDFileUMesh* mesh)
3363 if ( _nodeFields.empty() && _cellFields.empty() ) return 0;
3366 std::set< std::string > usedFieldNames;
3367 setFieldLongNames(usedFieldNames);
3369 MEDFileFields* fields = MEDFileFields::New();
3371 for ( unsigned int i = 0; i < _nodeFields.size(); ++i )
3372 setFields( _nodeFields[i], fields, mesh, i+1, usedFieldNames );
3374 for ( unsigned int i = 0; i < _cellFields.size(); ++i )
3375 setFields( _cellFields[i], fields, mesh, i+1, usedFieldNames );
3380 //================================================================================
3382 * \brief Make med fields from a SauvUtilities::DoubleField
3384 //================================================================================
3386 void IntermediateMED::setFields( SauvUtilities::DoubleField* fld,
3387 MEDCoupling::MEDFileFields* medFields,
3388 MEDCoupling::MEDFileUMesh* mesh,
3390 std::set< std::string >& usedFieldNames)
3392 bool sameNbGauss = true;
3393 if ( !fld || !fld->isMedCompatible( sameNbGauss )) return;
3396 fld->splitSubWithDiffNbGauss();
3398 // if ( !fld->hasCommonSupport() ):
3399 // each sub makes MEDFileFieldMultiTS
3401 // unite several subs into a MEDCouplingFieldDouble
3403 const bool uniteSubs = fld->hasCommonSupport() && sameNbGauss;
3405 std::cout << "Castem field #" << castemID << " <" << fld->_name
3406 << "> is incompatible with MED format, so we split it into several fields:" << std::endl;
3408 for ( unsigned int iSub = 0; iSub < fld->_sub.size(); )
3411 if ( !uniteSubs || fld->_name.empty() )
3412 makeFieldNewName( usedFieldNames, fld );
3415 DataArrayDouble * values = DataArrayDouble::New();
3416 values->alloc( fld->getNbTuples(iSub), fld->_sub[iSub].nbComponents() );
3419 double * valPtr = values->getPointer();
3422 mcIdType nbElems = fld->_group->size();
3423 for ( mcIdType elemShift = 0; elemShift < nbElems && iSub < fld->_sub.size(); )
3424 elemShift += fld->setValues( valPtr, iSub++, elemShift );
3425 setTS( fld, values, medFields, mesh );
3429 fld->setValues( valPtr, iSub );
3430 setTS( fld, values, medFields, mesh, iSub++ );
3432 std::cout << fld->_name << " with components";
3433 for ( size_t i = 0; i < (size_t)fld->_sub[iSub-1].nbComponents(); ++i )
3434 std::cout << " " << fld->_sub[iSub-1]._comp_names[ i ];
3435 std::cout << std::endl;
3440 //================================================================================
3442 * \brief Store value array of a field into med fields
3444 //================================================================================
3446 void IntermediateMED::setTS( SauvUtilities::DoubleField* fld,
3447 MEDCoupling::DataArrayDouble* values,
3448 MEDCoupling::MEDFileFields* medFields,
3449 MEDCoupling::MEDFileUMesh* mesh,
3452 // treat a field support
3453 const Group* support = fld->getSupport( iSub );
3455 const bool onAll = isOnAll( support, dimRel );
3456 if ( !onAll && support->_name.empty() )
3458 const_cast<Group*>(support)->_name += "PFL_" + fld->_name;
3459 support->_medGroup->setName( support->_name.c_str() );
3462 // make and fill a time-stamp
3464 MEDCouplingFieldDouble * timeStamp = MEDCouplingFieldDouble::New( fld->getMedType( iSub ),
3465 fld->getMedTimeDisc() );
3466 timeStamp->setName( fld->_name.c_str() );
3467 timeStamp->setDescription( fld->_description.c_str() );
3472 < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3473 timeStamp->setMesh( dimMesh );
3475 else if ( timeStamp->getTypeOfField() == MEDCoupling::ON_NODES )
3477 DataArrayDouble * coo = mesh->getCoords();
3479 <DataArrayDouble> subCoo = coo->selectByTupleId(support->_medGroup->begin(),
3480 support->_medGroup->end());
3481 MCAuto< MEDCouplingUMesh > nodeSubMesh =
3482 MEDCouplingUMesh::Build0DMeshFromCoords( subCoo );
3483 timeStamp->setMesh( nodeSubMesh );
3488 < MEDCouplingUMesh > dimMesh = mesh->getMeshAtLevel( dimRel );
3490 <MEDCouplingMesh> subMesh = dimMesh->buildPart(support->_medGroup->begin(),
3491 support->_medGroup->end());
3492 timeStamp->setMesh( subMesh);
3495 for ( size_t i = 0; i < (size_t)fld->_sub[iSub].nbComponents(); ++i )
3496 values->setInfoOnComponent( i, fld->_sub[iSub]._comp_names[ i ].c_str() );
3497 timeStamp->setArray( values );
3500 if ( timeStamp->getTypeOfField() == MEDCoupling::ON_GAUSS_PT )
3502 TGaussDef gaussDef( fld->_sub[iSub]._support->_cellType,
3503 fld->_sub[iSub].nbGauss() );
3504 timeStamp->setGaussLocalizationOnType( fld->_sub[iSub]._support->_cellType,
3505 gaussDef.myRefCoords,
3507 gaussDef.myWeights );
3509 // get a field to add the time-stamp
3510 bool isNewMedField = false;
3511 if ( !fld->_curMedField || fld->_name != fld->_curMedField->getName() )
3513 fld->_curMedField = MEDFileFieldMultiTS::New();
3514 isNewMedField = true;
3518 const int nbTS = fld->_curMedField->getNumberOfTS();
3520 timeStamp->setOrder( nbTS );
3522 // add the time-stamp
3523 timeStamp->checkConsistencyLight();
3525 fld->_curMedField->appendFieldNoProfileSBT( timeStamp );
3527 fld->_curMedField->appendFieldProfile( timeStamp, mesh, dimRel, support->_medGroup );
3528 timeStamp->decrRef();
3530 if ( isNewMedField ) // timeStamp must be added before this
3532 medFields->pushField( fld->_curMedField );
3536 //================================================================================
3538 * \brief Make a new unique name for a field
3540 //================================================================================
3542 void IntermediateMED::makeFieldNewName(std::set< std::string >& usedNames,
3543 SauvUtilities::DoubleField* fld )
3545 std::string base = fld->_name;
3552 std::string::size_type pos = base.rfind('_');
3553 if ( pos != std::string::npos )
3554 base = base.substr( 0, pos+1 );
3562 fld->_name = base + SauvUtilities::toString( i++ );
3564 while( !usedNames.insert( fld->_name ).second );
3567 //================================================================================
3569 * \brief Split sub-components with different nb of gauss points into several sub-components
3571 //================================================================================
3573 void DoubleField::splitSubWithDiffNbGauss()
3575 for ( size_t iSub = 0; iSub < _sub.size(); ++iSub )
3577 if ( _sub[iSub].isSameNbGauss() ) continue;
3579 _sub.insert( _sub.begin() + iSub + 1, 1, _Sub_data() );
3580 _Sub_data & subToSplit = _sub[iSub];
3581 _Sub_data & subNew = _sub[iSub+1];
3583 while ( subToSplit._nb_gauss[ 0 ] == subToSplit._nb_gauss[ iDiff ] )
3585 subNew._support = subToSplit._support;
3586 subNew._comp_names.assign( subToSplit._comp_names.begin() + iDiff,
3587 subToSplit._comp_names.end() );
3588 subNew._nb_gauss.assign ( subToSplit._nb_gauss.begin() + iDiff,
3589 subToSplit._nb_gauss.end() );
3590 subToSplit._comp_names.resize( iDiff );
3591 subToSplit._nb_gauss.resize ( iDiff );
3595 //================================================================================
3597 * \brief Return a vector ready to fill in
3599 //================================================================================
3601 std::vector< double >& DoubleField::addComponent( int nb_values )
3603 _comp_values.push_back( std::vector< double >() );
3604 std::vector< double >& res = _comp_values.back();
3605 res.resize( nb_values );
3609 DoubleField::~DoubleField()
3612 _curMedField->decrRef();
3615 //================================================================================
3617 * \brief Returns a supporting group
3619 //================================================================================
3621 const Group* DoubleField::getSupport( const int iSub ) const
3623 return _group ? _group : _sub[iSub]._support;
3626 //================================================================================
3628 * \brief Return true if each sub-component is a time stamp
3630 //================================================================================
3632 bool DoubleField::isMultiTimeStamps() const
3634 if ( _sub.size() < 2 )
3636 bool sameSupports = true;
3637 Group* grpp1 = _sub[0]._support;// grpp NOT grp because XDR under Windows defines grp...
3638 for ( size_t i = 1; i < _sub.size() && sameSupports; ++i )
3639 sameSupports = ( grpp1 == _sub[i]._support );
3641 return sameSupports;
3644 //================================================================================
3646 * \brief True if the field can be converted into the med field
3648 //================================================================================
3650 bool DoubleField::isMedCompatible(bool& sameNbGauss) const
3652 for ( unsigned int iSub = 0; iSub < _sub.size(); ++iSub )
3654 if ( !getSupport(iSub) || !getSupport(iSub)->_medGroup )
3655 THROW_IK_EXCEPTION("SauvReader INTERNAL ERROR: NULL field support");
3658 if ( !_sub[iSub].isSameNbGauss() )
3660 std::cout << "Field <" << _name << "> : different nb of gauss points in components" << std::endl;
3661 sameNbGauss = false;
3668 //================================================================================
3670 * \brief return true if all sub-components has same components and same nbGauss
3672 //================================================================================
3674 bool DoubleField::hasSameComponentsBySupport() const
3676 std::vector< _Sub_data >::const_iterator sub_data = _sub.begin();
3677 const _Sub_data& first_sub_data = *sub_data;
3678 for ( ++sub_data ; sub_data != _sub.end(); ++sub_data )
3680 if ( first_sub_data._comp_names != sub_data->_comp_names )
3681 return false; // diff names of components
3683 if ( first_sub_data._nb_gauss != sub_data->_nb_gauss &&
3684 first_sub_data._support->_cellType == sub_data->_support->_cellType)
3685 return false; // diff nb of gauss points on same cell type
3690 //================================================================================
3692 * \brief Return type of MEDCouplingFieldDouble
3694 //================================================================================
3696 MEDCoupling::TypeOfField DoubleField::getMedType( const int iSub ) const
3698 using namespace INTERP_KERNEL;
3700 const Group* grp = hasCommonSupport() ? _group : _sub[iSub]._support;
3701 if ( _sub[iSub].nbGauss() > 1 )
3703 const CellModel& cm = CellModel::GetCellModel( _sub[iSub]._support->_cellType );
3704 return (int) cm.getNumberOfNodes() == _sub[iSub].nbGauss() ? ON_GAUSS_NE : ON_GAUSS_PT;
3708 return getDim( grp ) == 0 ? ON_NODES : ON_CELLS;
3712 //================================================================================
3714 * \brief Return TypeOfTimeDiscretization
3716 //================================================================================
3718 MEDCoupling::TypeOfTimeDiscretization DoubleField::getMedTimeDisc() const
3724 // CONST_ON_TIME_INTERVAL = 7
3727 //================================================================================
3729 * \brief Return nb tuples to be used to allocate DataArrayDouble
3731 //================================================================================
3733 mcIdType DoubleField::getNbTuples( const int iSub ) const
3736 if ( hasCommonSupport() && !_group->_groups.empty() )
3737 for ( size_t i = 0; i < _group->_groups.size(); ++i )
3738 nb += _sub[i].nbGauss() * _sub[i]._support->size();
3740 nb = _sub[iSub].nbGauss() * getSupport(iSub)->size();
3744 //================================================================================
3746 * \brief Store values of a sub-component and return nb of elements in the iSub
3748 //================================================================================
3750 mcIdType DoubleField::setValues( double * valPtr, const int iSub, const mcIdType elemShift ) const
3752 // find values for iSub
3754 for ( int iS = 0; iS < iSub; ++iS )
3755 iComp += _sub[iS].nbComponents();
3756 const std::vector< double > * compValues = &_comp_values[ iComp ];
3760 const std::vector< mcIdType >& relocTable = getSupport( iSub )->_relocTable;
3762 const mcIdType nbElems = _sub[iSub]._support->size();
3763 const int nbGauss = _sub[iSub].nbGauss();
3764 const int nbComponents = _sub[iSub].nbComponents();
3765 const int nbValsByElem = nbComponents * nbGauss;
3768 mcIdType nbVals = 0;
3769 for ( iComp = 0; iComp < nbComponents; ++iComp )
3770 nbVals += ToIdType( compValues[iComp].size() );
3771 const bool isConstField = ( nbVals == nbComponents ); // one value per component (issue 22321)
3772 if ( !isConstField && nbVals != nbElems * nbValsByElem )
3773 THROW_IK_EXCEPTION("SauvMedConvertor.cxx: support size mismatches field size");
3775 // compute nb values in previous subs
3776 mcIdType valsShift = 0;
3777 for ( mcIdType iS = iSub-1, shift = elemShift; shift > 0; --iS)
3779 mcIdType nbE = _sub[iS]._support->size();
3781 valsShift += nbE * _sub[iS].nbComponents() * _sub[iS].nbGauss();
3785 for ( mcIdType iE = 0; iE < nbElems; ++iE )
3787 mcIdType iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3788 for ( iComp = 0; iComp < nbComponents; ++iComp )
3789 valPtr[ iMed + iComp ] = compValues[iComp][ 0 ];
3792 for ( mcIdType iE = 0; iE < nbElems; ++iE )
3794 mcIdType iMed = valsShift + nbValsByElem * ( relocTable.empty() ? iE : relocTable[iE+elemShift]-elemShift );
3795 for ( iComp = 0; iComp < nbComponents; ++iComp )
3796 for ( int iG = 0; iG < nbGauss; ++iG )
3797 valPtr[ iMed + iG * nbComponents + iComp ] = compValues[iComp][ iE * nbGauss + iG ];
3802 //================================================================================
3804 * \brief Destructor of IntermediateMED
3806 //================================================================================
3808 IntermediateMED::~IntermediateMED()
3810 for ( size_t i = 0; i < _nodeFields.size(); ++i )
3811 if ( _nodeFields[i] )
3812 delete _nodeFields[i];
3813 _nodeFields.clear();
3815 for ( size_t i = 0; i < _cellFields.size(); ++i )
3816 if ( _cellFields[i] )
3817 delete _cellFields[i];
3818 _cellFields.clear();
3820 for ( size_t i = 0; i < _groups.size(); ++i )
3821 if ( _groups[i]._medGroup )
3822 _groups[i]._medGroup->decrRef();
3825 //================================================================================
3827 * \brief CellsByDimIterator constructor
3829 CellsByDimIterator::CellsByDimIterator( const IntermediateMED & medi, int dimm)
3835 * \brief Initialize iteration on cells of given dimension
3837 void CellsByDimIterator::init(const int dimm)
3840 myTypeEnd = INTERP_KERNEL::NORM_HEXA20 + 1;
3844 * \brief return next set of Cell's of required dimension
3846 const std::set< Cell > * CellsByDimIterator::nextType()
3848 while ( ++myCurType < myTypeEnd )
3849 if ( !myImed->_cellsByType[myCurType].empty() && ( myDim < 0 || dim(false) == myDim ))
3850 return & myImed->_cellsByType[myCurType];
3854 * \brief return dimension of cells returned by the last or further next()
3856 int CellsByDimIterator::dim(const bool last) const
3858 int typp = myCurType;
3860 while ( typp < myTypeEnd && myImed->_cellsByType[typp].empty() )
3862 return typp < myTypeEnd ? getDimension( TCellType( typp )) : 4;
3864 // END CellsByDimIterator ========================================================