1 // Copyright (C) 2007-2016 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, 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
21 #include "InterpKernelGaussCoords.hxx"
22 #include "CellModel.hxx"
29 using namespace INTERP_KERNEL;
31 const double GaussInfo::TETRA4A_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0};
33 const double GaussInfo::TETRA4B_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0};
35 const double GaussInfo::TETRA10A_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0};
37 const double GaussInfo::TETRA10B_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5};
40 //Define common part of the code in the MACRO
41 //---------------------------------------------------------------
42 #define LOCAL_COORD_MACRO_BEGIN \
43 _my_local_reference_coord.resize( _my_local_ref_dim*_my_local_nb_ref ); \
44 for( int refId = 0; refId < _my_local_nb_ref; refId++ ) \
46 double* coords = &_my_local_reference_coord[ refId*_my_local_ref_dim ]; \
50 //---------------------------------------------------------------
51 #define LOCAL_COORD_MACRO_END \
55 //---------------------------------------------------------------
56 #define SHAPE_FUN_MACRO_BEGIN \
57 for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ ) \
59 double* funValue = &_my_function_value[ gaussId * _my_nb_ref ]; \
60 const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ];
62 //---------------------------------------------------------------
63 #define SHAPE_FUN_MACRO_END \
69 std::ostringstream stream; \
70 stream << "Error in the gauss localization for the cell with type "; \
71 stream << cellModel.getRepr(); \
73 throw INTERP_KERNEL::Exception(stream.str().c_str()); \
77 //---------------------------------------------------------------
78 static bool IsEqual(double theLeft, double theRight)
80 static double EPS = 1.0E-3;
81 if(fabs(theLeft) + fabs(theRight) > EPS)
82 return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
87 ////////////////////////////////////////////////////////////////////////////////////////////////
88 // GAUSS INFO CLASS //
89 ////////////////////////////////////////////////////////////////////////////////////////////////
92 * Constructor of the GaussInfo
94 GaussInfo::GaussInfo( NormalizedCellType theGeometry,
95 const DataVector& theGaussCoord,
97 const DataVector& theReferenceCoord,
99 _my_geometry(theGeometry),
100 _my_nb_gauss(theNbGauss),
101 _my_gauss_coord(theGaussCoord),
102 _my_nb_ref(theNbRef),
103 _my_reference_coord(theReferenceCoord)
106 //Allocate shape function values
107 _my_function_value.resize( _my_nb_gauss * _my_nb_ref );
113 GaussInfo::~GaussInfo()
118 * Return dimension of the gauss coordinates
120 int GaussInfo::getGaussCoordDim() const
124 return _my_gauss_coord.size()/_my_nb_gauss;
133 * Return dimension of the reference coordinates
135 int GaussInfo::getReferenceCoordDim() const
139 return _my_reference_coord.size()/_my_nb_ref;
148 * Return type of the cell.
150 NormalizedCellType GaussInfo::getCellType() const
156 * Return Nb of the gauss points.
158 int GaussInfo::getNbGauss() const
164 * Return Nb of the reference coordinates.
166 int GaussInfo::getNbRef() const
171 GaussInfo GaussInfo::convertToLinear() const
177 std::vector<double> a(TETRA10A_REF,TETRA10A_REF+30),b(TETRA10B_REF,TETRA10B_REF+30);
178 if(IsSatisfy(a,_my_reference_coord))
180 std::vector<double> c(TETRA4A_REF,TETRA4A_REF+12);
181 return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
183 if(IsSatisfy(b,_my_reference_coord))
185 std::vector<double> c(TETRA4B_REF,TETRA4B_REF+12);
186 return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
188 throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TETRA10 !");
191 throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not implemented yet for other types than TETRA10 !");
196 bool GaussInfo::IsSatisfy(const std::vector<double>& ref1, const std::vector<double>& ref2)
198 std::size_t sz(ref1.size());
201 for(std::size_t i=0;i<sz;i++)
202 if(!IsEqual(ref1[i],ref2[i]))
210 bool GaussInfo::isSatisfy()
213 bool anIsSatisfy = ((_my_local_nb_ref == _my_nb_ref) && (_my_local_ref_dim == getReferenceCoordDim()));
217 for( int refId = 0; refId < _my_local_nb_ref; refId++ )
219 double* refCoord = &_my_reference_coord[ refId*_my_local_ref_dim ];
220 double* localRefCoord = &_my_local_reference_coord[ refId*_my_local_ref_dim ];
221 bool anIsEqual = false;
222 for( int dimId = 0; dimId < _my_local_ref_dim; dimId++ )
224 anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
235 std::vector<double> GaussInfo::NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray)
237 std::size_t sz(inputArray.size()),dim((std::size_t)inputDim);
239 throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid dimension ! Must be !=0 !");
241 throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid input array ! Inconsistent with the given dimension !");
242 const CellModel& cm(CellModel::GetCellModel(ct));
243 std::size_t baseDim((std::size_t)cm.getDimension());
246 std::size_t nbOfItems(sz/dim);
247 std::vector<double> ret(nbOfItems*baseDim);
250 for(std::size_t i=0;i<nbOfItems;i++)
254 ret[i*baseDim+j]=inputArray[i*dim+j];
261 for(std::size_t i=0;i<nbOfItems;i++)
265 ret[i*baseDim+j]=inputArray[i*dim+j];
271 typedef void (*MapToShapeFunction)(GaussInfo& obj);
274 * Initialize the internal vectors
276 void GaussInfo::initLocalInfo()
278 bool aSatify = false;
279 const CellModel& cellModel(CellModel::GetCellModel(_my_geometry));
280 switch( _my_geometry )
283 _my_local_ref_dim = 0;
284 _my_local_nb_ref = 1;
286 aSatify = isSatisfy();
291 _my_local_ref_dim = 1;
292 _my_local_nb_ref = 2;
294 aSatify = isSatisfy();
299 _my_local_ref_dim = 1;
300 _my_local_nb_ref = 3;
302 aSatify = isSatisfy();
307 _my_local_ref_dim = 2;
308 _my_local_nb_ref = 3;
310 aSatify = isSatisfy();
315 aSatify = isSatisfy();
321 _my_local_ref_dim = 2;
322 _my_local_nb_ref = 6;
324 aSatify = isSatisfy();
328 aSatify = isSatisfy();
334 _my_local_ref_dim = 2;
335 _my_local_nb_ref = 7;
337 aSatify = isSatisfy();
343 _my_local_ref_dim = 2;
344 _my_local_nb_ref = 4;
345 MapToShapeFunction QUAD4PTR[]={Quad4aInit,Quad4bInit,Quad4cInit,Quad4DegSeg2Init};
346 std::size_t NB_OF_QUAD4PTR(sizeof(QUAD4PTR)/sizeof(MapToShapeFunction));
347 for(std::size_t i=0;i<NB_OF_QUAD4PTR && !aSatify;i++)
349 (QUAD4PTR[i])(*this);
350 aSatify = isSatisfy();
358 _my_local_ref_dim = 2;
359 _my_local_nb_ref = 8;
361 aSatify = isSatisfy();
366 aSatify = isSatisfy();
372 _my_local_ref_dim = 2;
373 _my_local_nb_ref = 9;
375 aSatify = isSatisfy();
380 _my_local_ref_dim = 3;
381 _my_local_nb_ref = 4;
383 aSatify = isSatisfy();
388 aSatify = isSatisfy();
394 _my_local_ref_dim = 3;
395 _my_local_nb_ref = 10;
397 aSatify = isSatisfy();
402 aSatify = isSatisfy();
408 _my_local_ref_dim = 3;
409 _my_local_nb_ref = 5;
411 aSatify = isSatisfy();
416 aSatify = isSatisfy();
422 _my_local_ref_dim = 3;
423 _my_local_nb_ref = 13;
425 aSatify = isSatisfy();
430 aSatify = isSatisfy();
437 _my_local_ref_dim = 3;
438 _my_local_nb_ref = 6;
439 MapToShapeFunction PENTA6PTR[]={Penta6aInit,Penta6bInit,Penta6DegTria3aInit,Penta6DegTria3bInit};
440 std::size_t NB_OF_PENTA6PTR(sizeof(PENTA6PTR)/sizeof(MapToShapeFunction));
441 for(std::size_t i=0;i<NB_OF_PENTA6PTR && !aSatify;i++)
443 (PENTA6PTR[i])(*this);
444 aSatify = isSatisfy();
451 _my_local_ref_dim = 3;
452 _my_local_nb_ref = 6;
454 aSatify = isSatisfy();
459 aSatify = isSatisfy();
466 _my_local_ref_dim = 3;
467 _my_local_nb_ref = 15;
468 MapToShapeFunction PENTA15PTR[]={Penta15aInit,Penta15bInit};
469 std::size_t NB_OF_PENTA15PTR(sizeof(PENTA15PTR)/sizeof(MapToShapeFunction));
470 for(std::size_t i=0;i<NB_OF_PENTA15PTR && !aSatify;i++)
472 (PENTA15PTR[i])(*this);
473 aSatify = isSatisfy();
481 _my_local_ref_dim = 3;
482 _my_local_nb_ref = 8;
483 MapToShapeFunction HEXA8PTR[]={Hexa8aInit,Hexa8bInit,Hexa8DegQuad4aInit,Hexa8DegQuad4bInit,Hexa8DegQuad4cInit};
484 std::size_t NB_OF_HEXA8PTR(sizeof(HEXA8PTR)/sizeof(MapToShapeFunction));
485 for(std::size_t i=0;i<NB_OF_HEXA8PTR && !aSatify;i++)
487 (HEXA8PTR[i])(*this);
488 aSatify = isSatisfy();
495 _my_local_ref_dim = 3;
496 _my_local_nb_ref = 20;
498 aSatify = isSatisfy();
503 aSatify = isSatisfy();
509 _my_local_ref_dim = 3;
510 _my_local_nb_ref = 27;
512 aSatify = isSatisfy();
517 throw INTERP_KERNEL::Exception("Not managed cell type !");
523 * Return shape function value by node id
525 const double* GaussInfo::getFunctionValues( const int theGaussId ) const
527 return &_my_function_value[ _my_nb_ref*theGaussId ];
530 void GaussInfo::point1Init()
532 double *funValue(&_my_function_value[0]);
537 * Init Segment 2 Reference coordinates ans Shape function.
539 void GaussInfo::seg2Init()
541 LOCAL_COORD_MACRO_BEGIN;
548 LOCAL_COORD_MACRO_END;
550 SHAPE_FUN_MACRO_BEGIN;
551 funValue[0] = 0.5*(1.0 - gc[0]);
552 funValue[1] = 0.5*(1.0 + gc[0]);
557 * Init Segment 3 Reference coordinates ans Shape function.
559 void GaussInfo::seg3Init()
561 LOCAL_COORD_MACRO_BEGIN;
571 LOCAL_COORD_MACRO_END;
573 SHAPE_FUN_MACRO_BEGIN;
574 funValue[0] = -0.5*(1.0 - gc[0])*gc[0];
575 funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
576 funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
581 * Init Triangle Reference coordinates ans Shape function.
584 void GaussInfo::tria3aInit()
586 LOCAL_COORD_MACRO_BEGIN;
599 LOCAL_COORD_MACRO_END;
601 SHAPE_FUN_MACRO_BEGIN;
602 funValue[0] = 0.5*(1.0 + gc[1]);
603 funValue[1] = -0.5*(gc[0] + gc[1]);
604 funValue[2] = 0.5*(1.0 + gc[0]);
609 * Init Triangle Reference coordinates ans Shape function.
612 void GaussInfo::tria3bInit()
614 LOCAL_COORD_MACRO_BEGIN;
627 LOCAL_COORD_MACRO_END;
629 SHAPE_FUN_MACRO_BEGIN;
630 funValue[0] = 1.0 - gc[0] - gc[1];
637 * Init Quadratic Triangle Reference coordinates ans Shape function.
640 void GaussInfo::tria6aInit()
642 LOCAL_COORD_MACRO_BEGIN;
667 LOCAL_COORD_MACRO_END;
669 SHAPE_FUN_MACRO_BEGIN;
670 funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
671 funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
672 funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
673 funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
674 funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
675 funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
680 * Init Quadratic Triangle Reference coordinates ans Shape function.
683 void GaussInfo::tria6bInit()
685 LOCAL_COORD_MACRO_BEGIN;
716 LOCAL_COORD_MACRO_END;
718 SHAPE_FUN_MACRO_BEGIN;
719 funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
720 funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
721 funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
722 funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
723 funValue[4] = 4.0*gc[0]*gc[1];
724 funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
728 void GaussInfo::tria7aInit()
730 LOCAL_COORD_MACRO_BEGIN;
756 coords[0] = 0.3333333333333333;
757 coords[1] = 0.3333333333333333;
760 LOCAL_COORD_MACRO_END;
762 SHAPE_FUN_MACRO_BEGIN;
763 funValue[0]=1-3*(gc[0]+gc[1])+2*(gc[0]*gc[0]+gc[1]*gc[1])+7*gc[0]*gc[1]-3*gc[0]*gc[1]*(gc[0]+gc[1]);
764 funValue[1]=gc[0]*(-1+2*gc[0]+3*gc[1]-3*gc[1]*(gc[0]+gc[1]));
765 funValue[2]=gc[1]*(-1.+3.*gc[0]+2.*gc[1]-3.*gc[0]*(gc[0]+gc[1]));
766 funValue[3]=4*gc[0]*(1-gc[0]-4*gc[1]+3*gc[1]*(gc[0]+gc[1]));
767 funValue[4]=4*gc[0]*gc[1]*(-2+3*(gc[0]+gc[1]));
768 funValue[5]=4*gc[1]*(1-4*gc[0]-gc[1]+3*gc[0]*(gc[0]+gc[1]));
769 funValue[6]=27*gc[0]*gc[1]*(1-gc[0]-gc[1]);
774 * Init Quadrangle Reference coordinates ans Shape function.
777 void GaussInfo::quad4aInit()
779 LOCAL_COORD_MACRO_BEGIN;
797 LOCAL_COORD_MACRO_END;
799 SHAPE_FUN_MACRO_BEGIN;
800 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
801 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
802 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
803 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
808 * Init Quadrangle Reference coordinates ans Shape function.
811 void GaussInfo::quad4bInit()
813 LOCAL_COORD_MACRO_BEGIN;
831 LOCAL_COORD_MACRO_END;
833 SHAPE_FUN_MACRO_BEGIN;
834 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
835 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
836 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
837 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
841 void GaussInfo::quad4cInit()
843 LOCAL_COORD_MACRO_BEGIN;
861 LOCAL_COORD_MACRO_END;
863 SHAPE_FUN_MACRO_BEGIN;
864 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
865 funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
866 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
867 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
872 * This shapefunc map is same as degenerated seg2Init
874 void GaussInfo::quad4DegSeg2Init()
876 LOCAL_COORD_MACRO_BEGIN;
893 LOCAL_COORD_MACRO_END;
895 SHAPE_FUN_MACRO_BEGIN;
896 funValue[0] = 0.5*(1.0 - gc[0]);
897 funValue[1] = 0.5*(1.0 + gc[0]);
904 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
907 void GaussInfo::quad8aInit()
909 LOCAL_COORD_MACRO_BEGIN;
942 LOCAL_COORD_MACRO_END;
944 SHAPE_FUN_MACRO_BEGIN;
945 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
946 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
947 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
948 funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
949 funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
950 funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
951 funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
952 funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
957 * Init Quadratic Quadrangle Reference coordinates ans Shape function.
960 void GaussInfo::quad8bInit()
962 LOCAL_COORD_MACRO_BEGIN;
995 LOCAL_COORD_MACRO_END;
997 SHAPE_FUN_MACRO_BEGIN;
998 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
999 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
1000 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
1001 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
1002 funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
1003 funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
1004 funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
1005 funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
1006 SHAPE_FUN_MACRO_END;
1009 void GaussInfo::quad9aInit()
1011 LOCAL_COORD_MACRO_BEGIN;
1048 LOCAL_COORD_MACRO_END;
1050 SHAPE_FUN_MACRO_BEGIN;
1051 funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
1052 funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
1053 funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
1054 funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
1055 funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
1056 funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
1057 funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
1058 funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
1059 funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
1060 SHAPE_FUN_MACRO_END;
1064 * Init Tetrahedron Reference coordinates ans Shape function.
1067 void GaussInfo::tetra4aInit()
1069 LOCAL_COORD_MACRO_BEGIN;
1071 coords[0] = TETRA4A_REF[0];
1072 coords[1] = TETRA4A_REF[1];
1073 coords[2] = TETRA4A_REF[2];
1076 coords[0] = TETRA4A_REF[3];
1077 coords[1] = TETRA4A_REF[4];
1078 coords[2] = TETRA4A_REF[5];
1081 coords[0] = TETRA4A_REF[6];
1082 coords[1] = TETRA4A_REF[7];
1083 coords[2] = TETRA4A_REF[8];
1086 coords[0] = TETRA4A_REF[9];
1087 coords[1] = TETRA4A_REF[10];
1088 coords[2] = TETRA4A_REF[11];
1090 LOCAL_COORD_MACRO_END;
1092 SHAPE_FUN_MACRO_BEGIN;
1093 funValue[0] = gc[1];
1094 funValue[1] = gc[2];
1095 funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
1096 funValue[3] = gc[0];
1097 SHAPE_FUN_MACRO_END;
1101 * Init Tetrahedron Reference coordinates ans Shape function.
1104 void GaussInfo::tetra4bInit()
1106 LOCAL_COORD_MACRO_BEGIN;
1108 coords[0] = TETRA4B_REF[0];
1109 coords[1] = TETRA4B_REF[1];
1110 coords[2] = TETRA4B_REF[2];
1113 coords[0] = TETRA4B_REF[3];
1114 coords[1] = TETRA4B_REF[4];
1115 coords[2] = TETRA4B_REF[5];
1118 coords[0] = TETRA4B_REF[6];
1119 coords[1] = TETRA4B_REF[7];
1120 coords[2] = TETRA4B_REF[8];
1123 coords[0] = TETRA4B_REF[9];
1124 coords[1] = TETRA4B_REF[10];
1125 coords[2] = TETRA4B_REF[11];
1127 LOCAL_COORD_MACRO_END;
1129 SHAPE_FUN_MACRO_BEGIN;
1130 funValue[0] = gc[1];
1131 funValue[2] = gc[2];
1132 funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
1133 funValue[3] = gc[0];
1134 SHAPE_FUN_MACRO_END;
1138 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1141 void GaussInfo::tetra10aInit()
1143 LOCAL_COORD_MACRO_BEGIN;
1145 coords[0] = TETRA10A_REF[0];
1146 coords[1] = TETRA10A_REF[1];
1147 coords[2] = TETRA10A_REF[2];
1150 coords[0] = TETRA10A_REF[3];
1151 coords[1] = TETRA10A_REF[4];
1152 coords[2] = TETRA10A_REF[5];
1155 coords[0] = TETRA10A_REF[6];
1156 coords[1] = TETRA10A_REF[7];
1157 coords[2] = TETRA10A_REF[8];
1160 coords[0] = TETRA10A_REF[9];
1161 coords[1] = TETRA10A_REF[10];
1162 coords[2] = TETRA10A_REF[11];
1165 coords[0] = TETRA10A_REF[12];
1166 coords[1] = TETRA10A_REF[13];
1167 coords[2] = TETRA10A_REF[14];
1170 coords[0] = TETRA10A_REF[15];
1171 coords[1] = TETRA10A_REF[16];
1172 coords[2] = TETRA10A_REF[17];
1175 coords[0] = TETRA10A_REF[18];
1176 coords[1] = TETRA10A_REF[19];
1177 coords[2] = TETRA10A_REF[20];
1180 coords[0] = TETRA10A_REF[21];
1181 coords[1] = TETRA10A_REF[22];
1182 coords[2] = TETRA10A_REF[23];
1185 coords[0] = TETRA10A_REF[24];
1186 coords[1] = TETRA10A_REF[25];
1187 coords[2] = TETRA10A_REF[26];
1190 coords[0] = TETRA10A_REF[27];
1191 coords[1] = TETRA10A_REF[28];
1192 coords[2] = TETRA10A_REF[29];
1194 LOCAL_COORD_MACRO_END;
1196 SHAPE_FUN_MACRO_BEGIN;
1197 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1198 funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
1199 funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1200 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1201 funValue[4] = 4.0*gc[1]*gc[2];
1202 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1203 funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1204 funValue[7] = 4.0*gc[0]*gc[1];
1205 funValue[8] = 4.0*gc[0]*gc[2];
1206 funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1207 SHAPE_FUN_MACRO_END;
1211 * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1214 void GaussInfo::tetra10bInit()
1216 LOCAL_COORD_MACRO_BEGIN;
1218 coords[0] = TETRA10B_REF[0];
1219 coords[1] = TETRA10B_REF[1];
1220 coords[2] = TETRA10B_REF[2];
1223 coords[0] = TETRA10B_REF[3];
1224 coords[1] = TETRA10B_REF[4];
1225 coords[2] = TETRA10B_REF[5];
1228 coords[0] = TETRA10B_REF[6];
1229 coords[1] = TETRA10B_REF[7];
1230 coords[2] = TETRA10B_REF[8];
1233 coords[0] = TETRA10B_REF[9];
1234 coords[1] = TETRA10B_REF[10];
1235 coords[2] = TETRA10B_REF[11];
1238 coords[0] = TETRA10B_REF[12];
1239 coords[1] = TETRA10B_REF[13];
1240 coords[2] = TETRA10B_REF[14];
1243 coords[0] = TETRA10B_REF[15];
1244 coords[1] = TETRA10B_REF[16];
1245 coords[2] = TETRA10B_REF[17];
1248 coords[0] = TETRA10B_REF[18];
1249 coords[1] = TETRA10B_REF[19];
1250 coords[2] = TETRA10B_REF[20];
1253 coords[0] = TETRA10B_REF[21];
1254 coords[1] = TETRA10B_REF[22];
1255 coords[2] = TETRA10B_REF[23];
1258 coords[0] = TETRA10B_REF[24];
1259 coords[1] = TETRA10B_REF[25];
1260 coords[2] = TETRA10B_REF[26];
1263 coords[0] = TETRA10B_REF[27];
1264 coords[1] = TETRA10B_REF[28];
1265 coords[2] = TETRA10B_REF[29];
1267 LOCAL_COORD_MACRO_END;
1268 SHAPE_FUN_MACRO_BEGIN;
1269 funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1270 funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1271 funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1272 funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1273 funValue[6] = 4.0*gc[1]*gc[2];
1274 funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1275 funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1276 funValue[7] = 4.0*gc[0]*gc[1];
1277 funValue[9] = 4.0*gc[0]*gc[2];
1278 funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1279 SHAPE_FUN_MACRO_END;
1283 * Init Pyramid Reference coordinates ans Shape function.
1286 void GaussInfo::pyra5aInit()
1288 LOCAL_COORD_MACRO_BEGIN;
1314 LOCAL_COORD_MACRO_END;
1316 SHAPE_FUN_MACRO_BEGIN;
1317 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1318 funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1319 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1320 funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1321 funValue[4] = gc[2];
1322 SHAPE_FUN_MACRO_END;
1325 * Init Pyramid Reference coordinates ans Shape function.
1328 void GaussInfo::pyra5bInit()
1330 LOCAL_COORD_MACRO_BEGIN;
1356 LOCAL_COORD_MACRO_END;
1358 SHAPE_FUN_MACRO_BEGIN;
1359 funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1360 funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1361 funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1362 funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1363 funValue[4] = gc[2];
1364 SHAPE_FUN_MACRO_END;
1368 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1371 void GaussInfo::pyra13aInit()
1373 LOCAL_COORD_MACRO_BEGIN;
1440 LOCAL_COORD_MACRO_END;
1442 SHAPE_FUN_MACRO_BEGIN;
1443 funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1444 (gc[0] - 0.5)/(1.0 - gc[2]);
1445 funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1446 (gc[1] - 0.5)/(1.0 - gc[2]);
1447 funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1448 (-gc[0] - 0.5)/(1.0 - gc[2]);
1449 funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1450 (-gc[1] - 0.5)/(1.0 - gc[2]);
1452 funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1454 funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1455 (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1456 funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1457 (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1458 funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1459 (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1460 funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1461 (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1463 funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1465 funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1467 funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1469 funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1471 SHAPE_FUN_MACRO_END;
1475 * Init Quadratic Pyramid Reference coordinates ans Shape function.
1478 void GaussInfo::pyra13bInit()
1480 LOCAL_COORD_MACRO_BEGIN;
1546 LOCAL_COORD_MACRO_END;
1548 SHAPE_FUN_MACRO_BEGIN;
1549 funValue[0] =0.5*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-0.5)/(1.0-gc[2]);
1550 funValue[1] =0.5*(+gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[1]-0.5)/(1.0-gc[2]);
1551 funValue[2] =0.5*(+gc[0]-gc[1]+gc[2]-1.0)*(+gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-0.5)/(1.0-gc[2]);
1552 funValue[3] =0.5*(-gc[0]-gc[1]+gc[2]-1.0)*(+gc[0]-gc[1]+gc[2]-1.0)*(gc[1]-0.5)/(1.0-gc[2]);
1554 funValue[4] =2.*gc[2]*(gc[2]-0.5);
1556 funValue[5] =-0.5*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1557 funValue[6] =-0.5*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1558 funValue[7] =-0.5*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1559 funValue[8] =-0.5*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1561 funValue[9] =gc[2]*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1562 funValue[10]=gc[2]*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1563 funValue[11]=gc[2]*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1564 funValue[12]=gc[2]*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1566 SHAPE_FUN_MACRO_END;
1571 * Init Pentahedron Reference coordinates and Shape function.
1574 void GaussInfo::penta6aInit()
1576 LOCAL_COORD_MACRO_BEGIN;
1607 LOCAL_COORD_MACRO_END;
1609 SHAPE_FUN_MACRO_BEGIN;
1610 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1611 funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
1612 funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1614 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1615 funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
1616 funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1617 SHAPE_FUN_MACRO_END;
1621 * Init Pentahedron Reference coordinates and Shape function.
1624 void GaussInfo::penta6bInit()
1626 LOCAL_COORD_MACRO_BEGIN;
1657 LOCAL_COORD_MACRO_END;
1659 SHAPE_FUN_MACRO_BEGIN;
1660 funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1661 funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
1662 funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1663 funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1664 funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
1665 funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1666 SHAPE_FUN_MACRO_END;
1670 * This shapefunc map is same as degenerated tria3aInit
1672 void GaussInfo::penta6DegTria3aInit()
1674 LOCAL_COORD_MACRO_BEGIN;
1705 LOCAL_COORD_MACRO_END;
1707 SHAPE_FUN_MACRO_BEGIN;
1708 funValue[0] = 0.5*(1.0 + gc[1]);
1709 funValue[1] = -0.5*(gc[0] + gc[1]);
1710 funValue[2] = 0.5*(1.0 + gc[0]);
1714 SHAPE_FUN_MACRO_END;
1718 * This shapefunc map is same as degenerated tria3bInit
1720 void GaussInfo::penta6DegTria3bInit()
1722 LOCAL_COORD_MACRO_BEGIN;
1753 LOCAL_COORD_MACRO_END;
1755 SHAPE_FUN_MACRO_BEGIN;
1756 funValue[0] = 1.0 - gc[0] - gc[1];
1757 funValue[1] = gc[0];
1758 funValue[2] = gc[1];
1762 SHAPE_FUN_MACRO_END;
1766 * Init Pentahedron Reference coordinates and Shape function.
1769 void GaussInfo::penta15aInit()
1771 LOCAL_COORD_MACRO_BEGIN;
1848 LOCAL_COORD_MACRO_END;
1850 SHAPE_FUN_MACRO_BEGIN;
1851 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1852 funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1853 funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1855 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1856 funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1857 funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1859 funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1860 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1861 funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1863 funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
1864 funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
1865 funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1867 funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1868 funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1869 funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1870 SHAPE_FUN_MACRO_END;
1874 * Init Qaudratic Pentahedron Reference coordinates and Shape function.
1877 void GaussInfo::penta15bInit()
1879 LOCAL_COORD_MACRO_BEGIN;
1956 LOCAL_COORD_MACRO_END;
1958 SHAPE_FUN_MACRO_BEGIN;
1959 funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1960 funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1961 funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1963 funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1964 funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1965 funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1967 funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1968 funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1969 funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1971 funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
1972 funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
1973 funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1975 funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1976 funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1977 funValue[9] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1978 SHAPE_FUN_MACRO_END;
1982 * Init Hehahedron Reference coordinates and Shape function.
1985 void GaussInfo::hexa8aInit()
1987 LOCAL_COORD_MACRO_BEGIN;
2028 LOCAL_COORD_MACRO_END;
2030 SHAPE_FUN_MACRO_BEGIN;
2031 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2032 funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2033 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2034 funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2036 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2037 funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2038 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2039 funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2040 SHAPE_FUN_MACRO_END;
2044 * Init Hehahedron Reference coordinates and Shape function.
2047 void GaussInfo::hexa8bInit()
2049 LOCAL_COORD_MACRO_BEGIN;
2090 LOCAL_COORD_MACRO_END;
2092 SHAPE_FUN_MACRO_BEGIN;
2093 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2094 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2095 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2096 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2098 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2099 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2100 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2101 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2102 SHAPE_FUN_MACRO_END;
2106 * This shapefunc map is same as degenerated quad4bInit
2108 void GaussInfo::hexa8DegQuad4aInit()
2110 LOCAL_COORD_MACRO_BEGIN;
2151 LOCAL_COORD_MACRO_END;
2153 SHAPE_FUN_MACRO_BEGIN;
2154 funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
2155 funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
2156 funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
2157 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2162 SHAPE_FUN_MACRO_END;
2166 * This shapefunc map is same as degenerated quad4bInit
2168 void GaussInfo::hexa8DegQuad4bInit()
2170 LOCAL_COORD_MACRO_BEGIN;
2211 LOCAL_COORD_MACRO_END;
2213 SHAPE_FUN_MACRO_BEGIN;
2214 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
2215 funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
2216 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2217 funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
2222 SHAPE_FUN_MACRO_END;
2226 * This shapefunc map is same as degenerated quad4cInit
2228 void GaussInfo::hexa8DegQuad4cInit()
2230 LOCAL_COORD_MACRO_BEGIN;
2272 LOCAL_COORD_MACRO_END;
2274 SHAPE_FUN_MACRO_BEGIN;
2275 funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
2276 funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
2277 funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2278 funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
2283 SHAPE_FUN_MACRO_END;
2287 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
2290 void GaussInfo::hexa20aInit()
2292 LOCAL_COORD_MACRO_BEGIN;
2394 LOCAL_COORD_MACRO_END;
2396 SHAPE_FUN_MACRO_BEGIN;
2397 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2398 (-2.0 - gc[0] - gc[1] - gc[2]);
2399 funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2400 (-2.0 + gc[0] - gc[1] - gc[2]);
2401 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2402 (-2.0 + gc[0] + gc[1] - gc[2]);
2403 funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2404 (-2.0 - gc[0] + gc[1] - gc[2]);
2405 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2406 (-2.0 - gc[0] - gc[1] + gc[2]);
2407 funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2408 (-2.0 + gc[0] - gc[1] + gc[2]);
2409 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2410 (-2.0 + gc[0] + gc[1] + gc[2]);
2411 funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2412 (-2.0 - gc[0] + gc[1] + gc[2]);
2414 funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2415 funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
2416 funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2417 funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
2418 funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
2419 funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
2420 funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
2421 funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
2422 funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2423 funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
2424 funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2425 funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
2426 SHAPE_FUN_MACRO_END;
2430 * Init Qaudratic Hehahedron Reference coordinates and Shape function.
2433 void GaussInfo::hexa20bInit()
2435 LOCAL_COORD_MACRO_BEGIN;
2537 LOCAL_COORD_MACRO_END;
2539 SHAPE_FUN_MACRO_BEGIN;
2541 funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2542 (-2.0 - gc[0] - gc[1] - gc[2]);
2543 funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2544 (-2.0 + gc[0] - gc[1] - gc[2]);
2545 funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2546 (-2.0 + gc[0] + gc[1] - gc[2]);
2547 funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2548 (-2.0 - gc[0] + gc[1] - gc[2]);
2549 funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2550 (-2.0 - gc[0] - gc[1] + gc[2]);
2551 funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2552 (-2.0 + gc[0] - gc[1] + gc[2]);
2553 funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2554 (-2.0 + gc[0] + gc[1] + gc[2]);
2555 funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2556 (-2.0 - gc[0] + gc[1] + gc[2]);
2558 funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2559 funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
2560 funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2561 funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
2562 funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
2563 funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
2564 funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
2565 funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
2566 funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2567 funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
2568 funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2569 funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
2570 SHAPE_FUN_MACRO_END;
2573 void GaussInfo::hexa27aInit()
2575 LOCAL_COORD_MACRO_BEGIN;
2711 LOCAL_COORD_MACRO_END;
2713 SHAPE_FUN_MACRO_BEGIN;
2715 funValue[0] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
2716 funValue[1] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
2717 funValue[2] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
2718 funValue[3] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
2719 funValue[4] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
2720 funValue[5] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
2721 funValue[6] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
2722 funValue[7] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
2723 funValue[8] =0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
2724 funValue[9] =0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
2725 funValue[10]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
2726 funValue[11]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
2727 funValue[12]=0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
2728 funValue[13]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
2729 funValue[14]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
2730 funValue[15]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
2731 funValue[16]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
2732 funValue[17]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
2733 funValue[18]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
2734 funValue[19]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
2735 funValue[20]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
2736 funValue[21]=0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
2737 funValue[22]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
2738 funValue[23]=0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
2739 funValue[24]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
2740 funValue[25]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
2741 funValue[26]=(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
2743 SHAPE_FUN_MACRO_END;
2746 ////////////////////////////////////////////////////////////////////////////////////////////////
2747 // GAUSS COORD CLASS //
2748 ////////////////////////////////////////////////////////////////////////////////////////////////
2752 GaussCoords::GaussCoords()
2759 GaussCoords::~GaussCoords()
2761 GaussInfoVector::iterator it = _my_gauss_info.begin();
2762 for( ; it != _my_gauss_info.end(); it++ )
2770 * Add Gauss localization info
2772 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
2774 const double* theGaussCoord,
2776 const double* theReferenceCoord,
2779 GaussInfoVector::iterator it = _my_gauss_info.begin();
2780 for( ; it != _my_gauss_info.end(); it++ )
2782 if( (*it)->getCellType() == theGeometry )
2788 DataVector aGaussCoord;
2789 for(int i = 0 ; i < theNbGauss*coordDim; i++ )
2790 aGaussCoord.push_back(theGaussCoord[i]);
2792 DataVector aReferenceCoord;
2793 for(int i = 0 ; i < theNbRef*coordDim; i++ )
2794 aReferenceCoord.push_back(theReferenceCoord[i]);
2797 GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
2798 info->initLocalInfo();
2800 //If info with cell type doesn't exist add it
2801 if( it == _my_gauss_info.end() )
2803 _my_gauss_info.push_back(info);
2805 // If information exists, update it
2809 int index = std::distance(_my_gauss_info.begin(),it);
2811 _my_gauss_info[index] = info;
2817 * Calculate gauss points coordinates
2819 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
2820 const double *theNodeCoords,
2821 const int theSpaceDim,
2822 const int *theIndex)
2824 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2825 int nbCoords = theSpaceDim * info->getNbGauss();
2826 double *aCoords = new double[nbCoords];
2827 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
2832 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2834 const GaussInfo *info = getInfoGivenCellType(theGeometry);
2835 calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
2838 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2840 int aConn = info->getNbRef();
2842 int nbCoords = theSpaceDim * info->getNbGauss();
2843 std::fill(result,result+nbCoords,0.);
2845 for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
2847 double *coord=result+gaussId*theSpaceDim;
2848 const double *function=info->getFunctionValues(gaussId);
2849 for ( int connId = 0; connId < aConn ; connId++ )
2851 const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
2852 for( int dimId = 0; dimId < theSpaceDim; dimId++ )
2853 coord[dimId] += nodeCoord[dimId]*function[connId];
2858 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
2860 GaussInfoVector::const_iterator it = _my_gauss_info.begin();
2861 //Try to find gauss localization info
2862 for( ; it != _my_gauss_info.end() ; it++ )
2863 if( (*it)->getCellType()==cellType)
2865 throw INTERP_KERNEL::Exception("Can't find gauss localization information !");