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()
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 = 2;
272 _my_local_nb_ref = 9;
274 aSatify = isSatisfy();
279 _my_local_ref_dim = 3;
280 _my_local_nb_ref = 4;
282 aSatify = isSatisfy();
287 aSatify = isSatisfy();
293 _my_local_ref_dim = 3;
294 _my_local_nb_ref = 10;
296 aSatify = isSatisfy();
301 aSatify = isSatisfy();
307 _my_local_ref_dim = 3;
308 _my_local_nb_ref = 5;
310 aSatify = isSatisfy();
315 aSatify = isSatisfy();
321 _my_local_ref_dim = 3;
322 _my_local_nb_ref = 13;
324 aSatify = isSatisfy();
329 aSatify = isSatisfy();
335 _my_local_ref_dim = 3;
336 _my_local_nb_ref = 6;
338 aSatify = isSatisfy();
343 aSatify = isSatisfy();
349 _my_local_ref_dim = 3;
350 _my_local_nb_ref = 15;
352 aSatify = isSatisfy();
357 aSatify = isSatisfy();
363 _my_local_ref_dim = 3;
364 _my_local_nb_ref = 8;
366 aSatify = isSatisfy();
371 aSatify = isSatisfy();
377 _my_local_ref_dim = 3;
378 _my_local_nb_ref = 20;
380 aSatify = isSatisfy();
385 aSatify = isSatisfy();
391 throw INTERP_KERNEL::Exception("Not manged cell type !");
397 * Return shape function value by node id
399 const double* GaussInfo::getFunctionValues( const int theGaussId ) const
401 return &_my_function_value[ _my_nb_ref*theGaussId ];
405 * Init Segment 2 Reference coordinates ans Shape function.
407 void GaussInfo::seg2Init()
409 LOCAL_COORD_MACRO_BEGIN;
416 LOCAL_COORD_MACRO_END;
418 SHAPE_FUN_MACRO_BEGIN;
419 funValue[0] = 0.5*(1.0 - gc[0]);
420 funValue[1] = 0.5*(1.0 + gc[0]);
425 * Init Segment 3 Reference coordinates ans Shape function.
427 void GaussInfo::seg3Init()
429 LOCAL_COORD_MACRO_BEGIN;
439 LOCAL_COORD_MACRO_END;
441 SHAPE_FUN_MACRO_BEGIN;
442 funValue[0] = -0.5*(1.0 - gc[0])*gc[0];
443 funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
444 funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
449 * Init Triangle Reference coordinates ans Shape function.
452 void GaussInfo::tria3aInit()
454 LOCAL_COORD_MACRO_BEGIN;
467 LOCAL_COORD_MACRO_END;
469 SHAPE_FUN_MACRO_BEGIN;
470 funValue[0] = 0.5*(1.0 + gc[1]);
471 funValue[1] = -0.5*(gc[0] + gc[1]);
472 funValue[2] = 0.5*(1.0 + gc[0]);
477 * Init Triangle Reference coordinates ans Shape function.
480 void GaussInfo::tria3bInit()
482 LOCAL_COORD_MACRO_BEGIN;
495 LOCAL_COORD_MACRO_END;
497 SHAPE_FUN_MACRO_BEGIN;
498 funValue[0] = 1.0 - gc[0] - gc[1];
505 * Init Quadratic Triangle Reference coordinates ans Shape function.
508 void GaussInfo::tria6aInit()
510 LOCAL_COORD_MACRO_BEGIN;
535 LOCAL_COORD_MACRO_END;
537 SHAPE_FUN_MACRO_BEGIN;
538 funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
539 funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
540 funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
541 funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
542 funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
543 funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
548 * Init Quadratic Triangle Reference coordinates ans Shape function.
551 void GaussInfo::tria6bInit()
553 LOCAL_COORD_MACRO_BEGIN;
584 LOCAL_COORD_MACRO_END;
586 SHAPE_FUN_MACRO_BEGIN;
587 funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
588 funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
589 funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
590 funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
591 funValue[4] = 4.0*gc[0]*gc[1];
592 funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
597 * Init Quadrangle Reference coordinates ans Shape function.
600 void GaussInfo::quad4aInit()
602 LOCAL_COORD_MACRO_BEGIN;
620 LOCAL_COORD_MACRO_END;
622 SHAPE_FUN_MACRO_BEGIN;
623 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
624 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
625 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
626 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
631 * Init Quadrangle Reference coordinates ans Shape function.
634 void GaussInfo::quad4bInit()
636 LOCAL_COORD_MACRO_BEGIN;
654 LOCAL_COORD_MACRO_END;
656 SHAPE_FUN_MACRO_BEGIN;
657 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
658 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
659 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
660 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
666 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
669 void GaussInfo::quad8aInit()
671 LOCAL_COORD_MACRO_BEGIN;
704 LOCAL_COORD_MACRO_END;
706 SHAPE_FUN_MACRO_BEGIN;
707 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
708 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
709 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
710 funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
711 funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
712 funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
713 funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
714 funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
719 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
722 void GaussInfo::quad8bInit()
724 LOCAL_COORD_MACRO_BEGIN;
757 LOCAL_COORD_MACRO_END;
759 SHAPE_FUN_MACRO_BEGIN;
760 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
761 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
762 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
763 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
764 funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
765 funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
766 funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
767 funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
771 void GaussInfo::quad9aInit()
773 LOCAL_COORD_MACRO_BEGIN;
810 LOCAL_COORD_MACRO_END;
812 SHAPE_FUN_MACRO_BEGIN;
813 funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
814 funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
815 funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
816 funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
817 funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
818 funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
819 funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
820 funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
821 funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
826 * Init Tetrahedron Reference coordinates ans Shape function.
829 void GaussInfo::tetra4aInit()
831 LOCAL_COORD_MACRO_BEGIN;
852 LOCAL_COORD_MACRO_END;
854 SHAPE_FUN_MACRO_BEGIN;
857 funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
863 * Init Tetrahedron Reference coordinates ans Shape function.
866 void GaussInfo::tetra4bInit()
868 LOCAL_COORD_MACRO_BEGIN;
889 LOCAL_COORD_MACRO_END;
891 SHAPE_FUN_MACRO_BEGIN;
894 funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
901 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
904 void GaussInfo::tetra10aInit()
906 LOCAL_COORD_MACRO_BEGIN;
957 LOCAL_COORD_MACRO_END;
959 SHAPE_FUN_MACRO_BEGIN;
960 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
961 funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
962 funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
963 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
964 funValue[4] = 4.0*gc[1]*gc[2];
965 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
966 funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
967 funValue[7] = 4.0*gc[0]*gc[1];
968 funValue[8] = 4.0*gc[0]*gc[2];
969 funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
974 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
977 void GaussInfo::tetra10bInit()
979 LOCAL_COORD_MACRO_BEGIN;
1030 LOCAL_COORD_MACRO_END;
1032 SHAPE_FUN_MACRO_BEGIN;
1033 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1034 funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1035 funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1036 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1037 funValue[6] = 4.0*gc[1]*gc[2];
1038 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1039 funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1040 funValue[7] = 4.0*gc[0]*gc[1];
1041 funValue[9] = 4.0*gc[0]*gc[2];
1042 funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1043 SHAPE_FUN_MACRO_END;
1047 * Init Pyramid Reference coordinates ans Shape function.
1050 void GaussInfo::pyra5aInit()
1052 LOCAL_COORD_MACRO_BEGIN;
1078 LOCAL_COORD_MACRO_END;
1080 SHAPE_FUN_MACRO_BEGIN;
1081 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1082 funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1083 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1084 funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1085 funValue[4] = gc[2];
1086 SHAPE_FUN_MACRO_END;
1089 * Init Pyramid Reference coordinates ans Shape function.
1092 void GaussInfo::pyra5bInit()
1094 LOCAL_COORD_MACRO_BEGIN;
1120 LOCAL_COORD_MACRO_END;
1122 SHAPE_FUN_MACRO_BEGIN;
1123 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1124 funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1125 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1126 funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1127 funValue[4] = gc[2];
1128 SHAPE_FUN_MACRO_END;
1132 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1135 void GaussInfo::pyra13aInit()
1137 LOCAL_COORD_MACRO_BEGIN;
1204 LOCAL_COORD_MACRO_END;
1206 SHAPE_FUN_MACRO_BEGIN;
1207 funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1208 (gc[0] - 0.5)/(1.0 - gc[2]);
1209 funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1210 (gc[1] - 0.5)/(1.0 - gc[2]);
1211 funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1212 (-gc[0] - 0.5)/(1.0 - gc[2]);
1213 funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1214 (-gc[1] - 0.5)/(1.0 - gc[2]);
1216 funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1218 funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1219 (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1220 funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1221 (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1222 funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1223 (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1224 funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1225 (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1227 funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1229 funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1231 funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1233 funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1235 SHAPE_FUN_MACRO_END;
1239 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1242 void GaussInfo::pyra13bInit()
1244 LOCAL_COORD_MACRO_BEGIN;
1310 LOCAL_COORD_MACRO_END;
1312 SHAPE_FUN_MACRO_BEGIN;
1313 funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1314 (gc[0] - 0.5)/(1.0 - gc[2]);
1315 funValue[3] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1316 (gc[1] - 0.5)/(1.0 - gc[2]);
1317 funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1318 (-gc[0] - 0.5)/(1.0 - gc[2]);
1319 funValue[1] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1320 (-gc[1] - 0.5)/(1.0 - gc[2]);
1322 funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1324 funValue[8] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1325 (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1326 funValue[7] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1327 (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1328 funValue[6] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1329 (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1330 funValue[5] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1331 (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1333 funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1335 funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1337 funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1339 funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1341 SHAPE_FUN_MACRO_END;
1346 * Init Pentahedron Reference coordinates and Shape function.
1349 void GaussInfo::penta6aInit()
1351 LOCAL_COORD_MACRO_BEGIN;
1382 LOCAL_COORD_MACRO_END;
1384 SHAPE_FUN_MACRO_BEGIN;
1385 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1386 funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
1387 funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1389 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1390 funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
1391 funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1392 SHAPE_FUN_MACRO_END;
1396 * Init Pentahedron Reference coordinates and Shape function.
1399 void GaussInfo::penta6bInit()
1401 LOCAL_COORD_MACRO_BEGIN;
1432 LOCAL_COORD_MACRO_END;
1434 SHAPE_FUN_MACRO_BEGIN;
1435 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1436 funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
1437 funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1438 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1439 funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
1440 funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1441 SHAPE_FUN_MACRO_END;
1444 * Init Pentahedron Reference coordinates and Shape function.
1447 void GaussInfo::penta15aInit()
1449 LOCAL_COORD_MACRO_BEGIN;
1526 LOCAL_COORD_MACRO_END;
1528 SHAPE_FUN_MACRO_BEGIN;
1529 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1530 funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1531 funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1533 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1534 funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1535 funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1537 funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1538 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1539 funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1541 funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
1542 funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
1543 funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1545 funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1546 funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1547 funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1548 SHAPE_FUN_MACRO_END;
1552 * Init Qaudratic Pentahedron Reference coordinates and Shape function.
1555 void GaussInfo::penta15bInit()
1557 LOCAL_COORD_MACRO_BEGIN;
1634 LOCAL_COORD_MACRO_END;
1636 SHAPE_FUN_MACRO_BEGIN;
1637 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1638 funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1639 funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1641 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1642 funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1643 funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1645 funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1646 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1647 funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1649 funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
1650 funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
1651 funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1653 funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1654 funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1655 funValue[9] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1656 SHAPE_FUN_MACRO_END;
1660 * Init Hehahedron Reference coordinates and Shape function.
1663 void GaussInfo::hexa8aInit()
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[1] = 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[3] = 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[5] = 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[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1718 SHAPE_FUN_MACRO_END;
1722 * Init Hehahedron Reference coordinates and Shape function.
1725 void GaussInfo::hexa8bInit()
1727 LOCAL_COORD_MACRO_BEGIN;
1768 LOCAL_COORD_MACRO_END;
1770 SHAPE_FUN_MACRO_BEGIN;
1771 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1772 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1773 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1774 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1776 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1777 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1778 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1779 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1780 SHAPE_FUN_MACRO_END;
1784 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1787 void GaussInfo::hexa20aInit()
1789 LOCAL_COORD_MACRO_BEGIN;
1891 LOCAL_COORD_MACRO_END;
1893 SHAPE_FUN_MACRO_BEGIN;
1894 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1895 (-2.0 - gc[0] - gc[1] - gc[2]);
1896 funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1897 (-2.0 + gc[0] - gc[1] - gc[2]);
1898 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1899 (-2.0 + gc[0] + gc[1] - gc[2]);
1900 funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1901 (-2.0 - gc[0] + gc[1] - gc[2]);
1902 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1903 (-2.0 - gc[0] - gc[1] + gc[2]);
1904 funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1905 (-2.0 + gc[0] - gc[1] + gc[2]);
1906 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1907 (-2.0 + gc[0] + gc[1] + gc[2]);
1908 funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1909 (-2.0 - gc[0] + gc[1] + gc[2]);
1911 funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1912 funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
1913 funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1914 funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
1915 funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
1916 funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
1917 funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
1918 funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
1919 funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1920 funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
1921 funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1922 funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
1923 SHAPE_FUN_MACRO_END;
1927 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1930 void GaussInfo::hexa20bInit()
1932 LOCAL_COORD_MACRO_BEGIN;
2034 LOCAL_COORD_MACRO_END;
2036 SHAPE_FUN_MACRO_BEGIN;
2038 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2039 (-2.0 - gc[0] - gc[1] - gc[2]);
2040 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2041 (-2.0 + gc[0] - gc[1] - gc[2]);
2042 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2043 (-2.0 + gc[0] + gc[1] - gc[2]);
2044 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2045 (-2.0 - gc[0] + gc[1] - gc[2]);
2046 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2047 (-2.0 - gc[0] - gc[1] + gc[2]);
2048 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2049 (-2.0 + gc[0] - gc[1] + gc[2]);
2050 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2051 (-2.0 + gc[0] + gc[1] + gc[2]);
2052 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2053 (-2.0 - gc[0] + gc[1] + gc[2]);
2055 funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2056 funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
2057 funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2058 funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
2059 funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
2060 funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
2061 funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
2062 funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
2063 funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2064 funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
2065 funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2066 funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
2067 SHAPE_FUN_MACRO_END;
2072 ////////////////////////////////////////////////////////////////////////////////////////////////
2073 // GAUSS COORD CLASS //
2074 ////////////////////////////////////////////////////////////////////////////////////////////////
2078 GaussCoords::GaussCoords()
2085 GaussCoords::~GaussCoords()
2087 GaussInfoVector::iterator it = _my_gauss_info.begin();
2088 for( ; it != _my_gauss_info.end(); it++ )
2096 * Add Gauss localization info
2098 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
2100 const double* theGaussCoord,
2102 const double* theReferenceCoord,
2105 GaussInfoVector::iterator it = _my_gauss_info.begin();
2106 for( ; it != _my_gauss_info.end(); it++ )
2108 if( (*it)->getCellType() == theGeometry )
2114 DataVector aGaussCoord;
2115 for(int i = 0 ; i < theNbGauss*coordDim; i++ )
2116 aGaussCoord.push_back(theGaussCoord[i]);
2118 DataVector aReferenceCoord;
2119 for(int i = 0 ; i < theNbRef*coordDim; i++ )
2120 aReferenceCoord.push_back(theReferenceCoord[i]);
2123 GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
2124 info->initLocalInfo();
2126 //If info with cell type doesn't exist add it
2127 if( it == _my_gauss_info.end() )
2129 _my_gauss_info.push_back(info);
2131 // If information exists, update it
2135 int index = std::distance(_my_gauss_info.begin(),it);
2137 _my_gauss_info[index] = info;
2143 * Calculate gauss points coordinates
2145 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
2146 const double *theNodeCoords,
2147 const int theSpaceDim,
2148 const int *theIndex)
2150 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2151 int nbCoords = theSpaceDim * info->getNbGauss();
2152 double *aCoords = new double[nbCoords];
2153 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
2158 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2160 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2161 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
2164 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2166 int aConn = info->getNbRef();
2168 int nbCoords = theSpaceDim * info->getNbGauss();
2169 std::fill(result,result+nbCoords,0.);
2171 for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
2173 double *coord=result+gaussId*theSpaceDim;
2174 const double *function=info->getFunctionValues(gaussId);
2175 for ( int connId = 0; connId < aConn ; connId++ )
2177 const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
2178 for( int dimId = 0; dimId < theSpaceDim; dimId++ )
2179 coord[dimId] += nodeCoord[dimId]*function[connId];
2184 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
2186 GaussInfoVector::const_iterator it = _my_gauss_info.begin();
2187 //Try to find gauss localization info
2188 for( ; it != _my_gauss_info.end() ; it++ )
2189 if( (*it)->getCellType()==cellType)
2191 throw INTERP_KERNEL::Exception("Can't find gauss localization information !");