1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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
21 #include "InterpKernelGaussCoords.hxx"
22 #include "CellModel.hxx"
29 using namespace INTERP_KERNEL;
31 //Define common part of the code in the MACRO
32 //---------------------------------------------------------------
33 #define LOCAL_COORD_MACRO_BEGIN \
34 _my_local_reference_coord.resize( _my_local_ref_dim*_my_local_nb_ref ); \
35 for( int refId = 0; refId < _my_local_nb_ref; refId++ ) \
37 double* coords = &_my_local_reference_coord[ refId*_my_local_ref_dim ]; \
41 //---------------------------------------------------------------
42 #define LOCAL_COORD_MACRO_END \
46 //---------------------------------------------------------------
47 #define SHAPE_FUN_MACRO_BEGIN \
48 for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ ) \
50 double* funValue = &_my_function_value[ gaussId * _my_nb_ref ]; \
51 const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ];
53 //---------------------------------------------------------------
54 #define SHAPE_FUN_MACRO_END \
60 std::ostringstream stream; \
61 stream << "Error in the gauss localization for the cell with type "; \
62 stream << cellModel.getRepr(); \
64 throw INTERP_KERNEL::Exception(stream.str().c_str()); \
68 //---------------------------------------------------------------
69 static bool IsEqual(double theLeft, double theRight)
71 static double EPS = 1.0E-3;
72 if(fabs(theLeft) + fabs(theRight) > EPS)
73 return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
78 ////////////////////////////////////////////////////////////////////////////////////////////////
79 // GAUSS INFO CLASS //
80 ////////////////////////////////////////////////////////////////////////////////////////////////
83 * Constructor of the GaussInfo
85 GaussInfo::GaussInfo( NormalizedCellType theGeometry,
86 const DataVector& theGaussCoord,
88 const DataVector& theReferenceCoord,
90 _my_geometry(theGeometry),
91 _my_nb_gauss(theNbGauss),
92 _my_gauss_coord(theGaussCoord),
94 _my_reference_coord(theReferenceCoord)
97 //Allocate shape function values
98 _my_function_value.resize( _my_nb_gauss * _my_nb_ref );
104 GaussInfo::~GaussInfo()
109 * Return dimension of the gauss coordinates
111 int GaussInfo::getGaussCoordDim() const
115 return _my_gauss_coord.size()/_my_nb_gauss;
124 * Return dimension of the reference coordinates
126 int GaussInfo::getReferenceCoordDim() const
130 return _my_reference_coord.size()/_my_nb_ref;
139 * Return type of the cell.
141 NormalizedCellType GaussInfo::getCellType() const
147 * Return Nb of the gauss points.
149 int GaussInfo::getNbGauss() const
155 * Return Nb of the reference coordinates.
157 int GaussInfo::getNbRef() const
165 bool GaussInfo::isSatisfy()
168 bool anIsSatisfy = ((_my_local_nb_ref == _my_nb_ref) && (_my_local_ref_dim == getReferenceCoordDim()));
172 for( int refId = 0; refId < _my_local_nb_ref; refId++ )
174 double* refCoord = &_my_reference_coord[ refId*_my_local_ref_dim ];
175 double* localRefCoord = &_my_local_reference_coord[ refId*_my_local_ref_dim ];
176 bool anIsEqual = false;
177 for( int dimId = 0; dimId < _my_local_ref_dim; dimId++ )
179 anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
191 * Initialize the internal vectors
193 void GaussInfo::initLocalInfo() throw (INTERP_KERNEL::Exception)
195 bool aSatify = false;
196 const CellModel& cellModel=CellModel::GetCellModel(_my_geometry);
197 switch( _my_geometry )
200 _my_local_ref_dim = 1;
201 _my_local_nb_ref = 2;
203 aSatify = isSatisfy();
208 _my_local_ref_dim = 1;
209 _my_local_nb_ref = 3;
211 aSatify = isSatisfy();
216 _my_local_ref_dim = 2;
217 _my_local_nb_ref = 3;
219 aSatify = isSatisfy();
224 aSatify = isSatisfy();
230 _my_local_ref_dim = 2;
231 _my_local_nb_ref = 6;
233 aSatify = isSatisfy();
237 aSatify = isSatisfy();
243 _my_local_ref_dim = 2;
244 _my_local_nb_ref = 4;
246 aSatify = isSatisfy();
251 aSatify = isSatisfy();
257 _my_local_ref_dim = 2;
258 _my_local_nb_ref = 8;
260 aSatify = isSatisfy();
265 aSatify = isSatisfy();
271 _my_local_ref_dim = 3;
272 _my_local_nb_ref = 4;
274 aSatify = isSatisfy();
279 aSatify = isSatisfy();
285 _my_local_ref_dim = 3;
286 _my_local_nb_ref = 10;
288 aSatify = isSatisfy();
293 aSatify = isSatisfy();
299 _my_local_ref_dim = 3;
300 _my_local_nb_ref = 5;
302 aSatify = isSatisfy();
307 aSatify = isSatisfy();
313 _my_local_ref_dim = 3;
314 _my_local_nb_ref = 13;
316 aSatify = isSatisfy();
321 aSatify = isSatisfy();
327 _my_local_ref_dim = 3;
328 _my_local_nb_ref = 6;
330 aSatify = isSatisfy();
335 aSatify = isSatisfy();
341 _my_local_ref_dim = 3;
342 _my_local_nb_ref = 15;
344 aSatify = isSatisfy();
349 aSatify = isSatisfy();
355 _my_local_ref_dim = 3;
356 _my_local_nb_ref = 8;
358 aSatify = isSatisfy();
363 aSatify = isSatisfy();
369 _my_local_ref_dim = 3;
370 _my_local_nb_ref = 20;
372 aSatify = isSatisfy();
377 aSatify = isSatisfy();
383 throw INTERP_KERNEL::Exception("Not manged cell type !");
389 * Return shape function value by node id
391 const double* GaussInfo::getFunctionValues( const int theGaussId ) const
393 return &_my_function_value[ _my_nb_ref*theGaussId ];
397 * Init Segment 2 Reference coordinates ans Shape function.
399 void GaussInfo::seg2Init()
401 LOCAL_COORD_MACRO_BEGIN;
408 LOCAL_COORD_MACRO_END;
410 SHAPE_FUN_MACRO_BEGIN;
411 funValue[0] = 0.5*(1.0 - gc[0]);
412 funValue[1] = 0.5*(1.0 + gc[0]);
417 * Init Segment 3 Reference coordinates ans Shape function.
419 void GaussInfo::seg3Init()
421 LOCAL_COORD_MACRO_BEGIN;
431 LOCAL_COORD_MACRO_END;
433 SHAPE_FUN_MACRO_BEGIN;
434 funValue[0] = 0.5*(1.0 - gc[0])*gc[0];
435 funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
436 funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
441 * Init Triangle Reference coordinates ans Shape function.
444 void GaussInfo::tria3aInit()
446 LOCAL_COORD_MACRO_BEGIN;
459 LOCAL_COORD_MACRO_END;
461 SHAPE_FUN_MACRO_BEGIN;
462 funValue[0] = 0.5*(1.0 + gc[1]);
463 funValue[1] = -0.5*(gc[0] + gc[1]);
464 funValue[2] = 0.5*(1.0 + gc[0]);
469 * Init Triangle Reference coordinates ans Shape function.
472 void GaussInfo::tria3bInit()
474 LOCAL_COORD_MACRO_BEGIN;
487 LOCAL_COORD_MACRO_END;
489 SHAPE_FUN_MACRO_BEGIN;
490 funValue[0] = 1.0 - gc[0] - gc[1];
497 * Init Quadratic Triangle Reference coordinates ans Shape function.
500 void GaussInfo::tria6aInit()
502 LOCAL_COORD_MACRO_BEGIN;
527 LOCAL_COORD_MACRO_END;
529 SHAPE_FUN_MACRO_BEGIN;
530 funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
531 funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
532 funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
533 funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
534 funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
535 funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
540 * Init Quadratic Triangle Reference coordinates ans Shape function.
543 void GaussInfo::tria6bInit()
545 LOCAL_COORD_MACRO_BEGIN;
576 LOCAL_COORD_MACRO_END;
578 SHAPE_FUN_MACRO_BEGIN;
579 funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
580 funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
581 funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
582 funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
583 funValue[4] = 4.0*gc[0]*gc[1];
584 funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
589 * Init Quadrangle Reference coordinates ans Shape function.
592 void GaussInfo::quad4aInit()
594 LOCAL_COORD_MACRO_BEGIN;
612 LOCAL_COORD_MACRO_END;
614 SHAPE_FUN_MACRO_BEGIN;
615 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
616 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
617 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
618 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
623 * Init Quadrangle Reference coordinates ans Shape function.
626 void GaussInfo::quad4bInit()
628 LOCAL_COORD_MACRO_BEGIN;
646 LOCAL_COORD_MACRO_END;
648 SHAPE_FUN_MACRO_BEGIN;
649 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
650 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
651 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
652 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
658 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
661 void GaussInfo::quad8aInit()
663 LOCAL_COORD_MACRO_BEGIN;
696 LOCAL_COORD_MACRO_END;
698 SHAPE_FUN_MACRO_BEGIN;
699 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
700 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
701 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
702 funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
703 funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
704 funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
705 funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
706 funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
711 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
714 void GaussInfo::quad8bInit()
716 LOCAL_COORD_MACRO_BEGIN;
749 LOCAL_COORD_MACRO_END;
751 SHAPE_FUN_MACRO_BEGIN;
752 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
753 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
754 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
755 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
756 funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
757 funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
758 funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
759 funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
764 * Init Tetrahedron Reference coordinates ans Shape function.
767 void GaussInfo::tetra4aInit()
769 LOCAL_COORD_MACRO_BEGIN;
790 LOCAL_COORD_MACRO_END;
792 SHAPE_FUN_MACRO_BEGIN;
795 funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
801 * Init Tetrahedron Reference coordinates ans Shape function.
804 void GaussInfo::tetra4bInit()
806 LOCAL_COORD_MACRO_BEGIN;
827 LOCAL_COORD_MACRO_END;
829 SHAPE_FUN_MACRO_BEGIN;
832 funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
839 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
842 void GaussInfo::tetra10aInit()
844 LOCAL_COORD_MACRO_BEGIN;
895 LOCAL_COORD_MACRO_END;
897 SHAPE_FUN_MACRO_BEGIN;
898 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
899 funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
900 funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
901 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
902 funValue[4] = 4.0*gc[1]*gc[2];
903 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
904 funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
905 funValue[7] = 4.0*gc[0]*gc[1];
906 funValue[8] = 4.0*gc[0]*gc[2];
907 funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
912 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
915 void GaussInfo::tetra10bInit()
917 LOCAL_COORD_MACRO_BEGIN;
968 LOCAL_COORD_MACRO_END;
970 SHAPE_FUN_MACRO_BEGIN;
971 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
972 funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
973 funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
974 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
975 funValue[6] = 4.0*gc[1]*gc[2];
976 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
977 funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
978 funValue[7] = 4.0*gc[0]*gc[1];
979 funValue[9] = 4.0*gc[0]*gc[2];
980 funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
985 * Init Pyramid Reference coordinates ans Shape function.
988 void GaussInfo::pyra5aInit()
990 LOCAL_COORD_MACRO_BEGIN;
1016 LOCAL_COORD_MACRO_END;
1018 SHAPE_FUN_MACRO_BEGIN;
1019 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1020 funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1021 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1022 funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1023 funValue[4] = gc[2];
1024 SHAPE_FUN_MACRO_END;
1027 * Init Pyramid Reference coordinates ans Shape function.
1030 void GaussInfo::pyra5bInit()
1032 LOCAL_COORD_MACRO_BEGIN;
1058 LOCAL_COORD_MACRO_END;
1060 SHAPE_FUN_MACRO_BEGIN;
1061 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1062 funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1063 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1064 funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1065 funValue[4] = gc[2];
1066 SHAPE_FUN_MACRO_END;
1070 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1073 void GaussInfo::pyra13aInit()
1075 LOCAL_COORD_MACRO_BEGIN;
1142 LOCAL_COORD_MACRO_END;
1144 SHAPE_FUN_MACRO_BEGIN;
1145 funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1146 (gc[0] - 0.5)/(1.0 - gc[2]);
1147 funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1148 (gc[1] - 0.5)/(1.0 - gc[2]);
1149 funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1150 (-gc[0] - 0.5)/(1.0 - gc[2]);
1151 funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1152 (-gc[1] - 0.5)/(1.0 - gc[2]);
1154 funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1156 funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1157 (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1158 funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1159 (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1160 funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1161 (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1162 funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1163 (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1165 funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1167 funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1169 funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1171 funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1173 SHAPE_FUN_MACRO_END;
1177 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1180 void GaussInfo::pyra13bInit()
1182 LOCAL_COORD_MACRO_BEGIN;
1248 LOCAL_COORD_MACRO_END;
1250 SHAPE_FUN_MACRO_BEGIN;
1251 funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1252 (gc[0] - 0.5)/(1.0 - gc[2]);
1253 funValue[3] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1254 (gc[1] - 0.5)/(1.0 - gc[2]);
1255 funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1256 (-gc[0] - 0.5)/(1.0 - gc[2]);
1257 funValue[1] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1258 (-gc[1] - 0.5)/(1.0 - gc[2]);
1260 funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1262 funValue[8] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1263 (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1264 funValue[7] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1265 (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1266 funValue[6] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1267 (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1268 funValue[5] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1269 (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1271 funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1273 funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1275 funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1277 funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1279 SHAPE_FUN_MACRO_END;
1284 * Init Pentahedron Reference coordinates and Shape function.
1287 void GaussInfo::penta6aInit()
1289 LOCAL_COORD_MACRO_BEGIN;
1320 LOCAL_COORD_MACRO_END;
1322 SHAPE_FUN_MACRO_BEGIN;
1323 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1324 funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
1325 funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1327 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1328 funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
1329 funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1330 SHAPE_FUN_MACRO_END;
1334 * Init Pentahedron Reference coordinates and Shape function.
1337 void GaussInfo::penta6bInit()
1339 LOCAL_COORD_MACRO_BEGIN;
1370 LOCAL_COORD_MACRO_END;
1372 SHAPE_FUN_MACRO_BEGIN;
1373 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1374 funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
1375 funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1376 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1377 funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
1378 funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1379 SHAPE_FUN_MACRO_END;
1382 * Init Pentahedron Reference coordinates and Shape function.
1385 void GaussInfo::penta15aInit()
1387 LOCAL_COORD_MACRO_BEGIN;
1464 LOCAL_COORD_MACRO_END;
1466 SHAPE_FUN_MACRO_BEGIN;
1467 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1468 funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1469 funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1471 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1472 funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1473 funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1475 funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1476 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1477 funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1479 funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
1480 funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
1481 funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1483 funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1484 funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1485 funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1486 SHAPE_FUN_MACRO_END;
1490 * Init Qaudratic Pentahedron Reference coordinates and Shape function.
1493 void GaussInfo::penta15bInit()
1495 LOCAL_COORD_MACRO_BEGIN;
1572 LOCAL_COORD_MACRO_END;
1574 SHAPE_FUN_MACRO_BEGIN;
1575 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1576 funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1577 funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1579 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1580 funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1581 funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1583 funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1584 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1585 funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1587 funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
1588 funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
1589 funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1591 funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1592 funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1593 funValue[9] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1594 SHAPE_FUN_MACRO_END;
1598 * Init Hehahedron Reference coordinates and Shape function.
1601 void GaussInfo::hexa8aInit()
1603 LOCAL_COORD_MACRO_BEGIN;
1644 LOCAL_COORD_MACRO_END;
1646 SHAPE_FUN_MACRO_BEGIN;
1647 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1648 funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1649 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1650 funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1652 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1653 funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1654 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1655 funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1656 SHAPE_FUN_MACRO_END;
1660 * Init Hehahedron Reference coordinates and Shape function.
1663 void GaussInfo::hexa8bInit()
1665 LOCAL_COORD_MACRO_BEGIN;
1706 LOCAL_COORD_MACRO_END;
1708 SHAPE_FUN_MACRO_BEGIN;
1709 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1710 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1711 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1712 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1714 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1715 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1716 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1717 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1718 SHAPE_FUN_MACRO_END;
1722 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1725 void GaussInfo::hexa20aInit()
1727 LOCAL_COORD_MACRO_BEGIN;
1829 LOCAL_COORD_MACRO_END;
1831 SHAPE_FUN_MACRO_BEGIN;
1832 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1833 (-2.0 - gc[0] - gc[1] - gc[2]);
1834 funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1835 (-2.0 + gc[0] - gc[1] - gc[2]);
1836 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1837 (-2.0 + gc[0] + gc[1] - gc[2]);
1838 funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1839 (-2.0 - gc[0] + gc[1] - gc[2]);
1840 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1841 (-2.0 - gc[0] - gc[1] + gc[2]);
1842 funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1843 (-2.0 + gc[0] - gc[1] + gc[2]);
1844 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1845 (-2.0 + gc[0] + gc[1] + gc[2]);
1846 funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1847 (-2.0 - gc[0] + gc[1] + gc[2]);
1849 funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1850 funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
1851 funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1852 funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
1853 funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
1854 funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
1855 funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
1856 funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
1857 funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1858 funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
1859 funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1860 funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
1861 SHAPE_FUN_MACRO_END;
1865 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1868 void GaussInfo::hexa20bInit()
1870 LOCAL_COORD_MACRO_BEGIN;
1972 LOCAL_COORD_MACRO_END;
1974 SHAPE_FUN_MACRO_BEGIN;
1976 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1977 (-2.0 - gc[0] - gc[1] - gc[2]);
1978 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1979 (-2.0 + gc[0] - gc[1] - gc[2]);
1980 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1981 (-2.0 + gc[0] + gc[1] - gc[2]);
1982 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1983 (-2.0 - gc[0] + gc[1] - gc[2]);
1984 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1985 (-2.0 - gc[0] - gc[1] + gc[2]);
1986 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1987 (-2.0 + gc[0] - gc[1] + gc[2]);
1988 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1989 (-2.0 + gc[0] + gc[1] + gc[2]);
1990 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1991 (-2.0 - gc[0] + gc[1] + gc[2]);
1993 funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1994 funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
1995 funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1996 funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
1997 funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
1998 funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
1999 funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
2000 funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
2001 funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2002 funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
2003 funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2004 funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
2005 SHAPE_FUN_MACRO_END;
2010 ////////////////////////////////////////////////////////////////////////////////////////////////
2011 // GAUSS COORD CLASS //
2012 ////////////////////////////////////////////////////////////////////////////////////////////////
2016 GaussCoords::GaussCoords()
2023 GaussCoords::~GaussCoords()
2025 GaussInfoVector::iterator it = _my_gauss_info.begin();
2026 for( ; it != _my_gauss_info.end(); it++ )
2034 * Add Gauss localization info
2036 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
2038 const double* theGaussCoord,
2040 const double* theReferenceCoord,
2041 int theNbRef) throw (INTERP_KERNEL::Exception)
2043 GaussInfoVector::iterator it = _my_gauss_info.begin();
2044 for( ; it != _my_gauss_info.end(); it++ )
2046 if( (*it)->getCellType() == theGeometry )
2052 DataVector aGaussCoord;
2053 for(int i = 0 ; i < theNbGauss*coordDim; i++ )
2054 aGaussCoord.push_back(theGaussCoord[i]);
2056 DataVector aReferenceCoord;
2057 for(int i = 0 ; i < theNbRef*coordDim; i++ )
2058 aReferenceCoord.push_back(theReferenceCoord[i]);
2061 GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
2062 info->initLocalInfo();
2064 //If info with cell type doesn't exist add it
2065 if( it == _my_gauss_info.end() )
2067 _my_gauss_info.push_back(info);
2069 // If information exists, update it
2073 int index = std::distance(_my_gauss_info.begin(),it);
2075 _my_gauss_info[index] = info;
2081 * Calculate gauss points coordinates
2083 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
2084 const double *theNodeCoords,
2085 const int theSpaceDim,
2086 const int *theIndex) throw (INTERP_KERNEL::Exception)
2088 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2089 int nbCoords = theSpaceDim * info->getNbGauss();
2090 double *aCoords = new double[nbCoords];
2091 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
2096 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const int *theIndex, double *result) throw(INTERP_KERNEL::Exception)
2098 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2099 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
2102 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2104 int aConn = info->getNbRef();
2106 int nbCoords = theSpaceDim * info->getNbGauss();
2107 std::fill(result,result+nbCoords,0.);
2109 for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
2111 double *coord=result+gaussId*theSpaceDim;
2112 const double *function=info->getFunctionValues(gaussId);
2113 for ( int connId = 0; connId < aConn ; connId++ )
2115 const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
2116 for( int dimId = 0; dimId < theSpaceDim; dimId++ )
2117 coord[dimId] += nodeCoord[dimId]*function[connId];
2122 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
2124 GaussInfoVector::const_iterator it = _my_gauss_info.begin();
2125 //Try to find gauss localization info
2126 for( ; it != _my_gauss_info.end() ; it++ )
2127 if( (*it)->getCellType()==cellType)
2129 throw INTERP_KERNEL::Exception("Can't find gauss localization information !");