src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 //Local includes
  21 #include "InterpKernelGaussCoords.hxx"
  22 #include "CellModel.hxx"
  23
  24 //STL includes
  25 #include <math.h>
  26 #include <algorithm>
  27 #include <sstream>
  28 #include <cmath>
  29
  30 using namespace INTERP_KERNEL;
  31
  32 const double GaussInfo::SEG2A_REF[2]={-1.0, 1.0};
  33
  34 const double GaussInfo::SEG2B_REF[2]={0., 1.0};
  35
  36 const double GaussInfo::SEG3_REF[3]={-1.0, 1.0, 0.0};
  37
  38 const double GaussInfo::TRIA3A_REF[6]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0};
  39
  40 const double GaussInfo::TRIA3B_REF[6]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
  41
  42 const double GaussInfo::TRIA6A_REF[12]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0};
  43
  44 const double GaussInfo::TRIA6B_REF[12]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5};
  45
  46 const double GaussInfo::TRIA7A_REF[14]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.3333333333333333, 0.3333333333333333};
  47
  48 const double GaussInfo::QUAD4A_REF[8]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
  49
  50 const double GaussInfo::QUAD4B_REF[8]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
  51
  52 const double GaussInfo::QUAD8A_REF[16]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0};
  53
  54 const double GaussInfo::QUAD8B_REF[16]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0};
  55
  56 const double GaussInfo::QUAD9A_REF[18]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0};
  57
  58 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};
  59
  60 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};
  61
  62 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};
  63
  64 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};
  65
  66 const double GaussInfo::PYRA5A_REF[15]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0};
  67
  68 const double GaussInfo::PYRA5B_REF[15]={1.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
  69
  70 const double GaussInfo::PYRA13A_REF[39]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, -0.5, 0.5, 0.0, -0.5, -0.5, 0.0, 0.5, -0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, -0.5, 0.0, 0.5, 0.0, -0.5, 0.5};
  71
  72 const double GaussInfo::PYRA13B_REF[39]={1.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, -0.5, 0.0, -0.5, -0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, -0.5, 0.5, -0.5, 0.0, 0.5, 0.0, 0.5, 0.5};
  73
  74 const double GaussInfo::PENTA6A_REF[18]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0};
  75
  76 const double GaussInfo::PENTA6B_REF[18]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0};
  77
  78 const double GaussInfo::PENTA15A_REF[45]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, 0.5, 0.5, -1.0, 0.0, 0.5, -1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0};
  79
  80 const double GaussInfo::PENTA15B_REF[45]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -1.0, 0.5, 0.0, -1.0, 0.0, 0.5, -1.0, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0};
  81
  82 const double GaussInfo::PENTA18A_REF[54]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, 0.5, 0.5, -1.0, 0.0, 0.5, -1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0};
  83
  84 const double GaussInfo::PENTA18B_REF[54]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -1.0, 0.5, 0.0, -1.0, 0.0, 0.5, -1.0, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5};
  85
  86 const double GaussInfo::HEXA8A_REF[24]={-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
  87
  88 const double GaussInfo::HEXA8B_REF[24]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0};
  89
  90 const double GaussInfo::HEXA20A_REF[60]={-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 0.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, -1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, -1.0, 0.0, 1.0};
  91
  92 const double GaussInfo::HEXA20B_REF[60]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 0.0, -1.0, 0.0, 1.0, -1.0, 1.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, -1.0, 1.0, -1.0, -1.0, 0.0, -1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, -1.0, 0.0};
  93
  94 const double GaussInfo::HEXA27A_REF[81]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 0.0, -1.0, 0.0, 1.0, -1.0, 1.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, -1.0, 1.0, -1.0, -1.0, 0.0, -1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 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};
  95
  96 //Define common part of the code in the MACRO
  97 //---------------------------------------------------------------
  98 #define LOCAL_COORD_MACRO_BEGIN                                         \
  99   _my_local_reference_coord.resize( _my_local_ref_dim*_my_local_nb_ref );           \
 100   for( int refId = 0; refId < _my_local_nb_ref; refId++ )                   \
 101     {                                                                   \
 102       double* coords = &_my_local_reference_coord[ refId*_my_local_ref_dim ];   \
 103       switch(refId)                                                     \
 104         {
 105
 106 //---------------------------------------------------------------
 107 #define LOCAL_COORD_MACRO_END                   \
 108   }                                             \
 109 }
 110
 111 //---------------------------------------------------------------
 112 #define SHAPE_FUN_MACRO_BEGIN                                              \
 113   for( int gaussId     = 0 ; gaussId < _my_nb_gauss ; gaussId++ )          \
 114     {                                                                      \
 115       double* funValue =  &_my_function_value[ gaussId * _my_nb_ref ];     \
 116       const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ];
 117
 118 //---------------------------------------------------------------
 119 #define SHAPE_FUN_MACRO_END                                                \
 120   }
 121
 122 #define DEV_SHAPE_FUN_MACRO_BEGIN                                                                               \
 123   for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ )                                                   \
 124     {                                                                                                           \
 125       double *devFunValue = _my_derivative_func_value.data() + gaussId * getReferenceCoordDim() * _my_nb_ref;   \
 126       const double *gc = _my_gauss_coord.data() + gaussId * getGaussCoordDim();
 127
 128 #define DEV_SHAPE_FUN_MACRO_END                                                                            \
 129     }
 130
 131 #define CHECK_MACRO                                                        \
 132   if( ! aSatify )                                                          \
 133     {                                                                      \
 134       std::ostringstream stream;                                           \
 135       stream << "Error in the gauss localization for the cell with type "; \
 136       stream << cellModel.getRepr();                                       \
 137       stream << " !!!";                                                    \
 138       throw INTERP_KERNEL::Exception(stream.str().c_str());                \
 139     }
 140
 141
 142 //---------------------------------------------------------------
 143 static bool IsEqual(double theLeft, double theRight)
 144 {
 145   static double EPS = 1.0E-3;
 146   if(fabs(theLeft) + fabs(theRight) > EPS)
 147     return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
 148   return true;
 149 }
 150
 151
 152 ////////////////////////////////////////////////////////////////////////////////////////////////
 153 //                                GAUSS INFO CLASS                                            //
 154 ////////////////////////////////////////////////////////////////////////////////////////////////
 155
 156 /*!
 157  * Constructor of the GaussInfo
 158  */
 159 GaussInfo::GaussInfo( NormalizedCellType theGeometry,
 160                       const DataVector& theGaussCoord,
 161                       int theNbGauss,
 162                       const DataVector& theReferenceCoord,
 163                       int theNbRef ) :
 164   _my_geometry(theGeometry),
 165   _my_nb_gauss(theNbGauss),
 166   _my_gauss_coord(theGaussCoord),
 167   _my_nb_ref(theNbRef),
 168   _my_reference_coord(theReferenceCoord)
 169 {
 170   //Allocate shape function values
 171   _my_function_value.resize( _my_nb_gauss * _my_nb_ref );
 172   _my_derivative_func_value.resize( _my_nb_gauss * _my_nb_ref * getReferenceCoordDim() );
 173 }
 174
 175 /*!
 176  * Destructor
 177  */
 178 GaussInfo::~GaussInfo()
 179 {
 180 }
 181
 182 /*!
 183  * Return dimension of the gauss coordinates
 184  */
 185 int GaussInfo::getGaussCoordDim() const
 186 {
 187   if( _my_nb_gauss )
 188     {
 189       return (int)_my_gauss_coord.size()/_my_nb_gauss;
 190     }
 191   else
 192     {
 193       return 0;
 194     }
 195 }
 196
 197 /*!
 198  * Return dimension of the reference coordinates
 199  */
 200 int GaussInfo::getReferenceCoordDim() const
 201 {
 202   if( _my_nb_ref )
 203     {
 204       return (int)(_my_reference_coord.size()/_my_nb_ref);
 205     }
 206   else
 207     {
 208       return 0;
 209     }
 210 }
 211
 212 /*!
 213  * Return type of the cell.
 214  */
 215 NormalizedCellType GaussInfo::getCellType() const
 216 {
 217   return _my_geometry;
 218 }
 219
 220 /*!
 221  * Return Nb of the gauss points.
 222  */
 223 int GaussInfo::getNbGauss() const
 224 {
 225   return _my_nb_gauss;
 226 }
 227
 228 /*!
 229  * Return Nb of the reference coordinates.
 230  */
 231 int GaussInfo::getNbRef() const
 232 {
 233   return _my_nb_ref;
 234 }
 235
 236 GaussInfo GaussInfo::convertToLinear() const
 237 {
 238   switch(_my_geometry)
 239     {
 240     case NORM_SEG3:
 241       {
 242         std::vector<double> a(SEG3_REF,SEG3_REF+3);
 243         if(IsSatisfy(a,_my_reference_coord))
 244           {
 245             std::vector<double> c(SEG2A_REF,SEG2A_REF+2);
 246             return GaussInfo(NORM_SEG2,_my_gauss_coord,getNbGauss(),c,2);
 247           }
 248         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for SEG3 !");
 249       }
 250     case NORM_TRI6:
 251       {
 252         std::vector<double> a(TRIA6A_REF,TRIA6A_REF+12),b(TRIA6B_REF,TRIA6B_REF+12);
 253         if(IsSatisfy(a,_my_reference_coord))
 254           {
 255             std::vector<double> c(TRIA3A_REF,TRIA3A_REF+6);
 256             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
 257           }
 258         if(IsSatisfy(b,_my_reference_coord))
 259           {
 260             std::vector<double> c(TRIA3B_REF,TRIA3B_REF+6);
 261             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
 262           }
 263         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TRI6 !");
 264       }
 265     case NORM_TRI7:
 266       {
 267         std::vector<double> a(TRIA7A_REF,TRIA7A_REF+14);
 268         if(IsSatisfy(a,_my_reference_coord))
 269           {
 270             std::vector<double> c(TRIA3B_REF,TRIA3B_REF+6);
 271             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
 272           }
 273         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TRI7 !");
 274       }
 275     case NORM_QUAD8:
 276       {
 277         std::vector<double> a(QUAD8A_REF,QUAD8A_REF+16),b(QUAD8B_REF,QUAD8B_REF+16);
 278         if(IsSatisfy(a,_my_reference_coord))
 279           {
 280             std::vector<double> c(QUAD4A_REF,QUAD4A_REF+8);
 281             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
 282           }
 283         if(IsSatisfy(b,_my_reference_coord))
 284           {
 285             std::vector<double> c(QUAD4B_REF,QUAD4B_REF+8);
 286             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
 287           }
 288         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for QUAD8 !");
 289       }
 290     case NORM_QUAD9:
 291       {
 292         std::vector<double> a(QUAD9A_REF,QUAD9A_REF+18);
 293         if(IsSatisfy(a,_my_reference_coord))
 294           {
 295             std::vector<double> c(QUAD4B_REF,QUAD4B_REF+8);
 296             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
 297           }
 298         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for QUAD9 !");
 299       }
 300     case NORM_TETRA10:
 301       {
 302         std::vector<double> a(TETRA10A_REF,TETRA10A_REF+30),b(TETRA10B_REF,TETRA10B_REF+30);
 303         if(IsSatisfy(a,_my_reference_coord))
 304           {
 305             std::vector<double> c(TETRA4A_REF,TETRA4A_REF+12);
 306             return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
 307           }
 308         if(IsSatisfy(b,_my_reference_coord))
 309           {
 310             std::vector<double> c(TETRA4B_REF,TETRA4B_REF+12);
 311             return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
 312           }
 313         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TETRA10 !");
 314       }
 315     case NORM_PYRA13:
 316       {
 317         std::vector<double> a(PYRA13A_REF,PYRA13A_REF+39),b(PYRA13B_REF,PYRA13B_REF+39);
 318         if(IsSatisfy(a,_my_reference_coord))
 319           {
 320             std::vector<double> c(PYRA5A_REF,PYRA5A_REF+15);
 321             return GaussInfo(NORM_PYRA5,_my_gauss_coord,getNbGauss(),c,5);
 322           }
 323         if(IsSatisfy(b,_my_reference_coord))
 324           {
 325             std::vector<double> c(PYRA5B_REF,PYRA5B_REF+15);
 326             return GaussInfo(NORM_PYRA5,_my_gauss_coord,getNbGauss(),c,5);
 327           }
 328         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PYRA13 !");
 329       }
 330     case NORM_PENTA15:
 331       {
 332         std::vector<double> a(PENTA15A_REF,PENTA15A_REF+45),b(PENTA15B_REF,PENTA15B_REF+45);
 333         if(IsSatisfy(a,_my_reference_coord))
 334           {
 335             std::vector<double> c(PENTA6A_REF,PENTA6A_REF+18);
 336             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
 337           }
 338         if(IsSatisfy(b,_my_reference_coord))
 339           {
 340             std::vector<double> c(PENTA6B_REF,PENTA6B_REF+18);
 341             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
 342           }
 343         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PENTA15 !");
 344       }
 345     case NORM_PENTA18:
 346       {
 347         std::vector<double> a(PENTA18A_REF,PENTA18A_REF+54),b(PENTA18B_REF,PENTA18B_REF+54);
 348         if(IsSatisfy(a,_my_reference_coord))
 349           {
 350             std::vector<double> c(PENTA6A_REF,PENTA6A_REF+18);
 351             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
 352           }
 353         if(IsSatisfy(b,_my_reference_coord))
 354           {
 355             std::vector<double> c(PENTA6B_REF,PENTA6B_REF+18);
 356             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
 357           }
 358         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PENTA18 !");
 359       }
 360     case NORM_HEXA20:
 361       {
 362         std::vector<double> a(HEXA20A_REF,HEXA20A_REF+60),b(HEXA20B_REF,HEXA20B_REF+60);
 363         if(IsSatisfy(a,_my_reference_coord))
 364           {
 365             std::vector<double> c(HEXA8A_REF,HEXA8A_REF+24);
 366             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
 367           }
 368         if(IsSatisfy(b,_my_reference_coord))
 369           {
 370             std::vector<double> c(HEXA8B_REF,HEXA8B_REF+24);
 371             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
 372           }
 373         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for HEXA20 !");
 374       }
 375     case NORM_HEXA27:
 376       {
 377         std::vector<double> a(HEXA27A_REF,HEXA27A_REF+81);
 378         if(IsSatisfy(a,_my_reference_coord))
 379           {
 380             std::vector<double> c(HEXA8B_REF,HEXA8B_REF+24);
 381             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
 382           }
 383         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for HEXA27 !");
 384       }
 385     default:
 386       throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not implemented yet for other types than TETRA10, HEXA20, HEXA27, TRI3, QUAD8, QUAD8, PYRA13, PENTA15 !");
 387     }
 388 }
 389
 390
 391 bool GaussInfo::IsSatisfy(const std::vector<double>& ref1, const std::vector<double>& ref2)
 392 {
 393   std::size_t sz(ref1.size());
 394   if(sz!=ref2.size())
 395     return false;
 396   for(std::size_t i=0;i<sz;i++)
 397     if(!IsEqual(ref1[i],ref2[i]))
 398       return false;
 399   return true;
 400 }
 401
 402 /*!
 403  * Check coordinates
 404  */
 405 bool GaussInfo::isSatisfy()
 406 {
 407
 408   bool anIsSatisfy = ((_my_local_nb_ref == _my_nb_ref) && (_my_local_ref_dim == getReferenceCoordDim()));
 409   //Check coordinates
 410   if(anIsSatisfy)
 411     {
 412       for( int refId = 0; refId < _my_local_nb_ref; refId++ )
 413         {
 414           double* refCoord = &_my_reference_coord[ refId*_my_local_ref_dim ];
 415           double* localRefCoord = &_my_local_reference_coord[ refId*_my_local_ref_dim ];
 416           bool anIsEqual = false;
 417           for( int dimId = 0; dimId < _my_local_ref_dim; dimId++ )
 418             {
 419               anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
 420               if(!anIsEqual )
 421                 {
 422                   return false;
 423                 }
 424             }
 425         }
 426     }
 427   return anIsSatisfy;
 428 }
 429
 430 std::vector<double> GaussInfo::NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray)
 431 {
 432   std::size_t sz(inputArray.size()),dim((std::size_t)inputDim);
 433   if(dim==0)
 434     throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid dimension ! Must be !=0 !");
 435   if(sz%dim!=0)
 436     throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid input array ! Inconsistent with the given dimension !");
 437   const CellModel& cm(CellModel::GetCellModel(ct));
 438   std::size_t baseDim((std::size_t)cm.getDimension());
 439   if(baseDim==dim)
 440     return inputArray;
 441   std::size_t nbOfItems(sz/dim);
 442   std::vector<double> ret(nbOfItems*baseDim);
 443   if(baseDim>dim)
 444     {
 445       for(std::size_t i=0;i<nbOfItems;i++)
 446         {
 447           std::size_t j=0;
 448           for(;j<dim;j++)
 449             ret[i*baseDim+j]=inputArray[i*dim+j];
 450           for(;j<baseDim;j++)
 451             ret[i*baseDim+j]=0.;
 452         }
 453     }
 454   else
 455     {
 456       for(std::size_t i=0;i<nbOfItems;i++)
 457         {
 458           std::size_t j=0;
 459           for(;j<baseDim;j++)
 460             ret[i*baseDim+j]=inputArray[i*dim+j];
 461         }
 462     }
 463   return ret;
 464 }
 465
 466 std::vector<double> GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellType ct)
 467 {
 468   switch(ct)
 469   {
 470     case INTERP_KERNEL::NORM_SEG2:
 471       return std::vector<double>(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double));
 472     case INTERP_KERNEL::NORM_SEG3:
 473       return std::vector<double>(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double));
 474     case INTERP_KERNEL::NORM_TRI3:
 475       return std::vector<double>(TRIA3A_REF,TRIA3A_REF+sizeof(TRIA3A_REF)/sizeof(double));
 476     case INTERP_KERNEL::NORM_TRI6:
 477       return std::vector<double>(TRIA6A_REF,TRIA6A_REF+sizeof(TRIA6A_REF)/sizeof(double));
 478     case INTERP_KERNEL::NORM_TRI7:
 479       return std::vector<double>(TRIA7A_REF,TRIA7A_REF+sizeof(TRIA7A_REF)/sizeof(double));
 480     case INTERP_KERNEL::NORM_QUAD4:
 481       return std::vector<double>(QUAD4A_REF,QUAD4A_REF+sizeof(QUAD4A_REF)/sizeof(double));
 482     case INTERP_KERNEL::NORM_QUAD8:
 483       return std::vector<double>(QUAD8A_REF,QUAD8A_REF+sizeof(QUAD8A_REF)/sizeof(double));
 484     case INTERP_KERNEL::NORM_QUAD9:
 485       return std::vector<double>(QUAD9A_REF,QUAD9A_REF+sizeof(QUAD9A_REF)/sizeof(double));
 486     case INTERP_KERNEL::NORM_TETRA4:
 487       return std::vector<double>(TETRA4A_REF,TETRA4A_REF+sizeof(TETRA4A_REF)/sizeof(double));
 488     case INTERP_KERNEL::NORM_TETRA10:
 489       return std::vector<double>(TETRA10A_REF,TETRA10A_REF+sizeof(TETRA10A_REF)/sizeof(double));
 490     case INTERP_KERNEL::NORM_PYRA5:
 491       return std::vector<double>(PYRA5A_REF,PYRA5A_REF+sizeof(PYRA5A_REF)/sizeof(double));
 492     case INTERP_KERNEL::NORM_PYRA13:
 493       return std::vector<double>(PYRA13A_REF,PYRA13A_REF+sizeof(PYRA13A_REF)/sizeof(double));
 494     case INTERP_KERNEL::NORM_PENTA6:
 495       return std::vector<double>(PENTA6A_REF,PENTA6A_REF+sizeof(PENTA6A_REF)/sizeof(double));
 496     case INTERP_KERNEL::NORM_PENTA15:
 497       return std::vector<double>(PENTA15A_REF,PENTA15A_REF+sizeof(PENTA15A_REF)/sizeof(double));
 498     case INTERP_KERNEL::NORM_PENTA18:
 499       return std::vector<double>(PENTA18A_REF,PENTA18A_REF+sizeof(PENTA18A_REF)/sizeof(double));
 500     case INTERP_KERNEL::NORM_HEXA8:
 501       return std::vector<double>(HEXA8A_REF,HEXA8A_REF+sizeof(HEXA8A_REF)/sizeof(double));
 502     case INTERP_KERNEL::NORM_HEXA20:
 503       return std::vector<double>(HEXA20A_REF,HEXA20A_REF+sizeof(HEXA20A_REF)/sizeof(double));
 504     case INTERP_KERNEL::NORM_HEXA27:
 505       return std::vector<double>(HEXA27A_REF,HEXA27A_REF+sizeof(HEXA27A_REF)/sizeof(double));
 506     default:
 507       THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf")
 508   }
 509 }
 510
 511 /*!
 512  * Returns true if \a ptInRefCoo is in reference cell of type \a ct (regarding GetDefaultReferenceCoordinatesOf)
 513  * \sa GetDefaultReferenceCoordinatesOf
 514  */
 515 bool GaussInfo::IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps)
 516 {
 517   switch(ct)
 518   {
 519     case INTERP_KERNEL::NORM_SEG2:
 520     case INTERP_KERNEL::NORM_SEG3:
 521     {
 522       return std::fabs(ptInRefCoo[0]) < 1.0+eps;
 523     }
 524     case INTERP_KERNEL::NORM_QUAD4:
 525     case INTERP_KERNEL::NORM_QUAD8:
 526     case INTERP_KERNEL::NORM_QUAD9:
 527     {
 528       return std::find_if(ptInRefCoo,ptInRefCoo+2,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+2;
 529     }
 530     case INTERP_KERNEL::NORM_HEXA8:
 531     case INTERP_KERNEL::NORM_HEXA20:
 532     case INTERP_KERNEL::NORM_HEXA27:
 533     {
 534       return std::find_if(ptInRefCoo,ptInRefCoo+3,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+3;
 535     }
 536     default:
 537       THROW_IK_EXCEPTION("IsInOrOutForReference : not implemented for this geo type !");
 538   }
 539 }
 540
 541 typedef void (*MapToShapeFunction)(GaussInfo& obj);
 542
 543 /*!
 544  * Initialize the internal vectors
 545  */
 546 void GaussInfo::initLocalInfo()
 547 {
 548   bool aSatify = false;
 549   const CellModel& cellModel(CellModel::GetCellModel(_my_geometry));
 550   switch( _my_geometry )
 551     {
 552     case NORM_POINT1:
 553       _my_local_ref_dim = 0;
 554       _my_local_nb_ref  = 1;
 555       point1Init();
 556       aSatify = isSatisfy();
 557       CHECK_MACRO;
 558       break;
 559
 560     case NORM_SEG2:
 561       _my_local_ref_dim = 1;
 562       _my_local_nb_ref  = 2;
 563       seg2aInit();
 564       aSatify = isSatisfy();
 565       if(!aSatify)
 566         {
 567           seg2bInit();
 568           aSatify = isSatisfy();
 569           CHECK_MACRO;
 570         }
 571       break;
 572
 573     case NORM_SEG3:
 574       _my_local_ref_dim = 1;
 575       _my_local_nb_ref  = 3;
 576       seg3Init();
 577       aSatify = isSatisfy();
 578       CHECK_MACRO;
 579       break;
 580
 581     case NORM_TRI3:
 582       _my_local_ref_dim = 2;
 583       _my_local_nb_ref  = 3;
 584       tria3aInit();
 585       aSatify = isSatisfy();
 586
 587       if(!aSatify)
 588         {
 589           tria3bInit();
 590           aSatify = isSatisfy();
 591           CHECK_MACRO;
 592         }
 593       break;
 594
 595     case NORM_TRI6:
 596       _my_local_ref_dim = 2;
 597       _my_local_nb_ref  = 6;
 598       tria6aInit();
 599       aSatify = isSatisfy();
 600       if(!aSatify)
 601         {
 602           tria6bInit();
 603           aSatify = isSatisfy();
 604           CHECK_MACRO;
 605         }
 606       break;
 607
 608     case NORM_TRI7:
 609       _my_local_ref_dim = 2;
 610       _my_local_nb_ref  = 7;
 611       tria7aInit();
 612       aSatify = isSatisfy();
 613       CHECK_MACRO;
 614       break;
 615
 616     case NORM_QUAD4:
 617       {
 618         _my_local_ref_dim = 2;
 619         _my_local_nb_ref  = 4;
 620         MapToShapeFunction QUAD4PTR[]={Quad4aInit,Quad4bInit,Quad4cInit,Quad4DegSeg2Init};
 621         std::size_t NB_OF_QUAD4PTR(sizeof(QUAD4PTR)/sizeof(MapToShapeFunction));
 622         for(std::size_t i=0;i<NB_OF_QUAD4PTR && !aSatify;i++)
 623           {
 624             (QUAD4PTR[i])(*this);
 625             aSatify = isSatisfy();
 626           }
 627         CHECK_MACRO;
 628         break;
 629       }
 630       break;
 631
 632     case NORM_QUAD8:
 633       _my_local_ref_dim = 2;
 634       _my_local_nb_ref  = 8;
 635       quad8aInit();
 636       aSatify = isSatisfy();
 637
 638       if(!aSatify)
 639         {
 640           quad8bInit();
 641           aSatify = isSatisfy();
 642           CHECK_MACRO;
 643         }
 644       break;
 645
 646     case NORM_QUAD9:
 647       _my_local_ref_dim = 2;
 648       _my_local_nb_ref  = 9;
 649       quad9aInit();
 650       aSatify = isSatisfy();
 651       CHECK_MACRO;
 652       break;
 653
 654     case NORM_TETRA4:
 655       _my_local_ref_dim = 3;
 656       _my_local_nb_ref  = 4;
 657       tetra4aInit();
 658       aSatify = isSatisfy();
 659
 660       if(!aSatify)
 661         {
 662           tetra4bInit();
 663           aSatify = isSatisfy();
 664           CHECK_MACRO;
 665         }
 666       break;
 667
 668     case NORM_TETRA10:
 669       _my_local_ref_dim = 3;
 670       _my_local_nb_ref  = 10;
 671       tetra10aInit();
 672       aSatify = isSatisfy();
 673
 674       if(!aSatify)
 675         {
 676           tetra10bInit();
 677           aSatify = isSatisfy();
 678           CHECK_MACRO;
 679         }
 680       break;
 681
 682     case NORM_PYRA5:
 683       _my_local_ref_dim = 3;
 684       _my_local_nb_ref  = 5;
 685       pyra5aInit();
 686       aSatify = isSatisfy();
 687
 688       if(!aSatify)
 689         {
 690           pyra5bInit();
 691           aSatify = isSatisfy();
 692           CHECK_MACRO;
 693         }
 694       break;
 695
 696     case NORM_PYRA13:
 697       _my_local_ref_dim = 3;
 698       _my_local_nb_ref  = 13;
 699       pyra13aInit();
 700       aSatify = isSatisfy();
 701
 702       if(!aSatify)
 703         {
 704           pyra13bInit();
 705           aSatify = isSatisfy();
 706           CHECK_MACRO;
 707         }
 708       break;
 709
 710     case NORM_PENTA6:
 711       {
 712         _my_local_ref_dim = 3;
 713         _my_local_nb_ref  = 6;
 714         MapToShapeFunction PENTA6PTR[]={Penta6aInit,Penta6bInit,Penta6DegTria3aInit,Penta6DegTria3bInit};
 715         std::size_t NB_OF_PENTA6PTR(sizeof(PENTA6PTR)/sizeof(MapToShapeFunction));
 716         for(std::size_t i=0;i<NB_OF_PENTA6PTR && !aSatify;i++)
 717           {
 718             (PENTA6PTR[i])(*this);
 719             aSatify = isSatisfy();
 720           }
 721         CHECK_MACRO;
 722         break;
 723       }
 724
 725 /*      _my_local_ref_dim = 3;
 726       _my_local_nb_ref  = 6;
 727       penta6aInit();
 728       aSatify = isSatisfy();
 729
 730       if(!aSatify)
 731         {
 732           penta6bInit();
 733           aSatify = isSatisfy();
 734           CHECK_MACRO;
 735         }
 736       break;
 737 */
 738     case NORM_PENTA15:
 739       {
 740         _my_local_ref_dim = 3;
 741         _my_local_nb_ref  = 15;
 742         MapToShapeFunction PENTA15PTR[]={Penta15aInit,Penta15bInit};
 743         std::size_t NB_OF_PENTA15PTR(sizeof(PENTA15PTR)/sizeof(MapToShapeFunction));
 744         for(std::size_t i=0;i<NB_OF_PENTA15PTR && !aSatify;i++)
 745           {
 746             (PENTA15PTR[i])(*this);
 747             aSatify = isSatisfy();
 748           }
 749         CHECK_MACRO;
 750         break;
 751       }
 752
 753     case NORM_PENTA18:
 754       {
 755         _my_local_ref_dim = 3;
 756         _my_local_nb_ref  = 18;
 757         MapToShapeFunction PENTA18PTR[]={Penta18aInit,Penta18bInit};
 758         std::size_t NB_OF_PENTA18PTR(sizeof(PENTA18PTR)/sizeof(MapToShapeFunction));
 759         for(std::size_t i=0;i<NB_OF_PENTA18PTR && !aSatify;i++)
 760           {
 761             (PENTA18PTR[i])(*this);
 762             aSatify = isSatisfy();
 763           }
 764         CHECK_MACRO;
 765         break;
 766       }
 767
 768     case NORM_HEXA8:
 769       {
 770         _my_local_ref_dim = 3;
 771         _my_local_nb_ref  = 8;
 772         MapToShapeFunction HEXA8PTR[]={Hexa8aInit,Hexa8bInit,Hexa8DegQuad4aInit,Hexa8DegQuad4bInit,Hexa8DegQuad4cInit};
 773         std::size_t NB_OF_HEXA8PTR(sizeof(HEXA8PTR)/sizeof(MapToShapeFunction));
 774         for(std::size_t i=0;i<NB_OF_HEXA8PTR && !aSatify;i++)
 775           {
 776             (HEXA8PTR[i])(*this);
 777             aSatify = isSatisfy();
 778           }
 779         CHECK_MACRO;
 780         break;
 781       }
 782
 783     case NORM_HEXA20:
 784       _my_local_ref_dim = 3;
 785       _my_local_nb_ref  = 20;
 786       hexa20aInit();
 787       aSatify = isSatisfy();
 788
 789       if(!aSatify)
 790         {
 791           hexa20bInit();
 792           aSatify = isSatisfy();
 793           CHECK_MACRO;
 794         }
 795       break;
 796
 797     case NORM_HEXA27:
 798       _my_local_ref_dim = 3;
 799       _my_local_nb_ref  = 27;
 800       hexa27aInit();
 801       aSatify = isSatisfy();
 802       CHECK_MACRO
 803       break;
 804
 805     default:
 806       throw INTERP_KERNEL::Exception("Not managed cell type !");
 807       break;
 808     }
 809 }
 810
 811 /**
 812  * Return shape function value by node id
 813  */
 814 const double *GaussInfo::getFunctionValues( const int theGaussId ) const
 815 {
 816   return _my_function_value.data() + _my_nb_ref*theGaussId ;
 817 }
 818
 819 /**
 820  * Return the derivative of shape function value by node id
 821  */
 822 const double *GaussInfo::getDerivativeOfShapeFunctionAt( const int theGaussId ) const
 823 {
 824   return _my_derivative_func_value.data() + _my_nb_ref*getReferenceCoordDim()*theGaussId;
 825 }
 826
 827 void GaussInfo::point1Init()
 828 {
 829   double *funValue(&_my_function_value[0]);
 830   funValue[0] = 1. ;
 831 }
 832
 833 /*!
 834  * Init Segment 2 Reference coordinates ans Shape function.
 835  */
 836 void GaussInfo::seg2aInit()
 837 {
 838   LOCAL_COORD_MACRO_BEGIN;
 839   case 0:
 840     coords[0] = SEG2A_REF[0];
 841     break;
 842   case 1:
 843     coords[0] = SEG2A_REF[1];
 844     break;
 845   LOCAL_COORD_MACRO_END;
 846
 847   SHAPE_FUN_MACRO_BEGIN;
 848   funValue[0] = 0.5*(1.0 - gc[0]);
 849   funValue[1] = 0.5*(1.0 + gc[0]);
 850   SHAPE_FUN_MACRO_END;
 851 }
 852
 853 void GaussInfo::seg2bInit()
 854 {
 855   LOCAL_COORD_MACRO_BEGIN;
 856   case 0:
 857     coords[0] = SEG2B_REF[0];
 858     break;
 859   case 1:
 860     coords[0] = SEG2B_REF[1];
 861     break;
 862   LOCAL_COORD_MACRO_END;
 863
 864   SHAPE_FUN_MACRO_BEGIN;
 865   funValue[0] = 1.0 - gc[0];
 866   funValue[1] = gc[0];
 867   SHAPE_FUN_MACRO_END;
 868 }
 869
 870 /*!
 871  * Init Segment 3 Reference coordinates ans Shape function.
 872  */
 873 void GaussInfo::seg3Init()
 874 {
 875   LOCAL_COORD_MACRO_BEGIN;
 876   case 0:
 877     coords[0] = SEG3_REF[0];
 878     break;
 879   case 1:
 880     coords[0] = SEG3_REF[1];
 881     break;
 882   case 2:
 883     coords[0] = SEG3_REF[2];
 884   LOCAL_COORD_MACRO_END;
 885
 886   SHAPE_FUN_MACRO_BEGIN;
 887   funValue[0] = -0.5*(1.0 - gc[0])*gc[0];
 888   funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
 889   funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
 890   SHAPE_FUN_MACRO_END;
 891 }
 892
 893 /*!
 894  * Init Triangle Reference coordinates ans Shape function.
 895  * Case A.
 896  */
 897 void GaussInfo::tria3aInit()
 898 {
 899   LOCAL_COORD_MACRO_BEGIN;
 900   case 0:
 901     coords[0] = TRIA3A_REF[0];
 902     coords[1] = TRIA3A_REF[1];
 903     break;
 904   case 1:
 905     coords[0] = TRIA3A_REF[2];
 906     coords[1] = TRIA3A_REF[3];
 907     break;
 908   case 2:
 909     coords[0] = TRIA3A_REF[4];
 910     coords[1] = TRIA3A_REF[5];
 911     break;
 912   LOCAL_COORD_MACRO_END;
 913
 914   SHAPE_FUN_MACRO_BEGIN;
 915   funValue[0] = 0.5*(1.0 + gc[1]);
 916   funValue[1] = -0.5*(gc[0] + gc[1]);
 917   funValue[2] = 0.5*(1.0 + gc[0]);
 918   SHAPE_FUN_MACRO_END;
 919 }
 920
 921 /*!
 922  * Init Triangle Reference coordinates ans Shape function.
 923  * Case B.
 924  */
 925 void GaussInfo::tria3bInit()
 926 {
 927   LOCAL_COORD_MACRO_BEGIN;
 928   case 0:
 929     coords[0] = TRIA3B_REF[0];
 930     coords[1] = TRIA3B_REF[1];
 931     break;
 932   case 1:
 933     coords[0] = TRIA3B_REF[2];
 934     coords[1] = TRIA3B_REF[3];
 935     break;
 936   case 2:
 937     coords[0] = TRIA3B_REF[4];
 938     coords[1] = TRIA3B_REF[5];
 939     break;
 940   LOCAL_COORD_MACRO_END;
 941
 942   SHAPE_FUN_MACRO_BEGIN;
 943   funValue[0] = 1.0 - gc[0] - gc[1];
 944   funValue[1] = gc[0];
 945   funValue[2] = gc[1];
 946   SHAPE_FUN_MACRO_END;
 947 }
 948
 949 /*!
 950  * Init Quadratic Triangle Reference coordinates ans Shape function.
 951  * Case A.
 952  */
 953 void GaussInfo::tria6aInit()
 954 {
 955   LOCAL_COORD_MACRO_BEGIN;
 956   case 0:
 957     coords[0] = TRIA6A_REF[0];
 958     coords[1] = TRIA6A_REF[1];
 959     break;
 960   case 1:
 961     coords[0] = TRIA6A_REF[2];
 962     coords[1] = TRIA6A_REF[3];
 963     break;
 964   case 2:
 965     coords[0] = TRIA6A_REF[4];
 966     coords[1] = TRIA6A_REF[5];
 967     break;
 968   case 3:
 969     coords[0] = TRIA6A_REF[6];
 970     coords[1] = TRIA6A_REF[7];
 971     break;
 972   case 4:
 973     coords[0] = TRIA6A_REF[8];
 974     coords[1] = TRIA6A_REF[9];
 975     break;
 976   case 5:
 977     coords[0] = TRIA6A_REF[10];
 978     coords[1] = TRIA6A_REF[11];
 979     break;
 980   LOCAL_COORD_MACRO_END;
 981
 982   SHAPE_FUN_MACRO_BEGIN;
 983   funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
 984   funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
 985   funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
 986   funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
 987   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
 988   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
 989   SHAPE_FUN_MACRO_END;
 990 }
 991
 992 /*!
 993  * Init Quadratic Triangle Reference coordinates ans Shape function.
 994  * Case B.
 995  */
 996 void GaussInfo::tria6bInit()
 997 {
 998   LOCAL_COORD_MACRO_BEGIN;
 999   case 0:
1000     coords[0] = TRIA6B_REF[0];
1001     coords[1] = TRIA6B_REF[1];
1002     break;
1003   case 1:
1004     coords[0] = TRIA6B_REF[2];
1005     coords[1] = TRIA6B_REF[3];
1006     break;
1007   case 2:
1008     coords[0] = TRIA6B_REF[4];
1009     coords[1] = TRIA6B_REF[5];
1010     break;
1011   case 3:
1012     coords[0] = TRIA6B_REF[6];
1013     coords[1] = TRIA6B_REF[7];
1014     break;
1015   case 4:
1016     coords[0] = TRIA6B_REF[8];
1017     coords[1] = TRIA6B_REF[9];
1018     break;
1019   case 5:
1020     coords[0] = TRIA6B_REF[10];
1021     coords[1] = TRIA6B_REF[11];
1022     break;
1023   LOCAL_COORD_MACRO_END;
1024
1025   SHAPE_FUN_MACRO_BEGIN;
1026   funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
1027   funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
1028   funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
1029   funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
1030   funValue[4] = 4.0*gc[0]*gc[1];
1031   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
1032   SHAPE_FUN_MACRO_END;
1033 }
1034
1035 void GaussInfo::tria7aInit()
1036 {
1037   LOCAL_COORD_MACRO_BEGIN;
1038   case 0:
1039     coords[0] = TRIA7A_REF[0];
1040     coords[1] = TRIA7A_REF[1];
1041     break;
1042   case 1:
1043     coords[0] = TRIA7A_REF[2];
1044     coords[1] = TRIA7A_REF[3];
1045     break;
1046   case 2:
1047     coords[0] = TRIA7A_REF[4];
1048     coords[1] = TRIA7A_REF[5];
1049     break;
1050   case 3:
1051     coords[0] = TRIA7A_REF[6];
1052     coords[1] = TRIA7A_REF[7];
1053     break;
1054   case 4:
1055     coords[0] = TRIA7A_REF[8];
1056     coords[1] = TRIA7A_REF[9];
1057     break;
1058   case 5:
1059     coords[0] = TRIA7A_REF[10];
1060     coords[1] = TRIA7A_REF[11];
1061     break;
1062   case 6:
1063     coords[0] = TRIA7A_REF[12];
1064     coords[1] = TRIA7A_REF[13];
1065     break;
1066   LOCAL_COORD_MACRO_END;
1067
1068   SHAPE_FUN_MACRO_BEGIN;
1069   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]);
1070   funValue[1]=gc[0]*(-1+2*gc[0]+3*gc[1]-3*gc[1]*(gc[0]+gc[1]));
1071   funValue[2]=gc[1]*(-1.+3.*gc[0]+2.*gc[1]-3.*gc[0]*(gc[0]+gc[1]));
1072   funValue[3]=4*gc[0]*(1-gc[0]-4*gc[1]+3*gc[1]*(gc[0]+gc[1]));
1073   funValue[4]=4*gc[0]*gc[1]*(-2+3*(gc[0]+gc[1]));
1074   funValue[5]=4*gc[1]*(1-4*gc[0]-gc[1]+3*gc[0]*(gc[0]+gc[1]));
1075   funValue[6]=27*gc[0]*gc[1]*(1-gc[0]-gc[1]);
1076   SHAPE_FUN_MACRO_END;
1077 }
1078
1079 /*!
1080  * Init Quadrangle Reference coordinates ans Shape function.
1081  * Case A.
1082  */
1083 void GaussInfo::quad4aInit()
1084 {
1085   LOCAL_COORD_MACRO_BEGIN;
1086   case 0:
1087     coords[0] = QUAD4A_REF[0];
1088     coords[1] = QUAD4A_REF[1];
1089     break;
1090   case 1:
1091     coords[0] = QUAD4A_REF[2];
1092     coords[1] = QUAD4A_REF[3];
1093     break;
1094   case 2:
1095     coords[0] = QUAD4A_REF[4];
1096     coords[1] = QUAD4A_REF[5];
1097     break;
1098   case 3:
1099     coords[0] = QUAD4A_REF[6];
1100     coords[1] = QUAD4A_REF[7];
1101     break;
1102   LOCAL_COORD_MACRO_END;
1103
1104   SHAPE_FUN_MACRO_BEGIN;
1105   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
1106   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
1107   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
1108   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1109   SHAPE_FUN_MACRO_END;
1110 }
1111
1112 /*!
1113  * Init Quadrangle Reference coordinates ans Shape function.
1114  * Case B.
1115  */
1116 void GaussInfo::quad4bInit()
1117 {
1118   LOCAL_COORD_MACRO_BEGIN;
1119   case 0:
1120     coords[0] = QUAD4B_REF[0];
1121     coords[1] = QUAD4B_REF[1];
1122     break;
1123   case 1:
1124     coords[0] = QUAD4B_REF[2];
1125     coords[1] = QUAD4B_REF[3];
1126     break;
1127   case 2:
1128     coords[0] = QUAD4B_REF[4];
1129     coords[1] = QUAD4B_REF[5];
1130     break;
1131   case 3:
1132     coords[0] = QUAD4B_REF[6];
1133     coords[1] = QUAD4B_REF[7];
1134     break;
1135   LOCAL_COORD_MACRO_END;
1136
1137   SHAPE_FUN_MACRO_BEGIN;
1138   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1139   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1140   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1141   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1142   SHAPE_FUN_MACRO_END;
1143 }
1144
1145 void GaussInfo::quad4cInit()
1146 {
1147   LOCAL_COORD_MACRO_BEGIN;
1148  case  0:
1149    coords[0] = -1.0;
1150    coords[1] = -1.0;
1151    break;
1152  case  1:
1153    coords[0] = -1.0;
1154    coords[1] =  1.0;
1155    break;
1156  case  2:
1157    coords[0] =  1.0;
1158    coords[1] =  1.0;
1159    break;
1160  case  3:
1161    coords[0] =  1.0;
1162    coords[1] = -1.0;
1163    break;
1164
1165    LOCAL_COORD_MACRO_END;
1166
1167    SHAPE_FUN_MACRO_BEGIN;
1168    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1169    funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1170    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1171    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1172    SHAPE_FUN_MACRO_END;
1173 }
1174
1175 /*!
1176  * This shapefunc map is same as degenerated seg2aInit
1177  */
1178 void GaussInfo::quad4DegSeg2Init()
1179 {
1180   LOCAL_COORD_MACRO_BEGIN;
1181  case  0:
1182    coords[0] = -1.0;
1183    coords[1] =  0.0;
1184    break;
1185  case  1:
1186    coords[0] =  1.0;
1187    coords[1] =  0.0;
1188    break;
1189  case  2:
1190    coords[0] =  0.0;
1191    coords[1] =  0.0;
1192    break;
1193  case  3:
1194    coords[0] =  0.0;
1195    coords[1] =  0.0;
1196    break;
1197    LOCAL_COORD_MACRO_END;
1198
1199    SHAPE_FUN_MACRO_BEGIN;
1200    funValue[0] = 0.5*(1.0 - gc[0]);
1201    funValue[1] = 0.5*(1.0 + gc[0]);
1202    funValue[2] = 0.;
1203    funValue[3] = 0.;
1204    SHAPE_FUN_MACRO_END;
1205 }
1206
1207 /*!
1208  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1209  * Case A.
1210  */
1211 void GaussInfo::quad8aInit()
1212 {
1213   LOCAL_COORD_MACRO_BEGIN;
1214   case 0:
1215     coords[0] = QUAD8A_REF[0];
1216     coords[1] = QUAD8A_REF[1];
1217     break;
1218   case 1:
1219     coords[0] = QUAD8A_REF[2];
1220     coords[1] = QUAD8A_REF[3];
1221     break;
1222   case 2:
1223     coords[0] = QUAD8A_REF[4];
1224     coords[1] = QUAD8A_REF[5];
1225     break;
1226   case 3:
1227     coords[0] = QUAD8A_REF[6];
1228     coords[1] = QUAD8A_REF[7];
1229     break;
1230   case 4:
1231     coords[0] = QUAD8A_REF[8];
1232     coords[1] = QUAD8A_REF[9];
1233     break;
1234   case 5:
1235     coords[0] = QUAD8A_REF[10];
1236     coords[1] = QUAD8A_REF[11];
1237     break;
1238   case 6:
1239     coords[0] = QUAD8A_REF[12];
1240     coords[1] = QUAD8A_REF[13];
1241     break;
1242   case 7:
1243     coords[0] = QUAD8A_REF[14];
1244     coords[1] = QUAD8A_REF[15];
1245     break;
1246   LOCAL_COORD_MACRO_END;
1247
1248   SHAPE_FUN_MACRO_BEGIN;
1249   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
1250   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
1251   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
1252   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
1253   funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1254   funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1255   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1256   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1257   SHAPE_FUN_MACRO_END;
1258 }
1259
1260 /*!
1261  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1262  * Case B.
1263  */
1264 void GaussInfo::quad8bInit()
1265 {
1266   LOCAL_COORD_MACRO_BEGIN;
1267   case 0:
1268     coords[0] = QUAD8B_REF[0];
1269     coords[1] = QUAD8B_REF[1];
1270     break;
1271   case 1:
1272     coords[0] = QUAD8B_REF[2];
1273     coords[1] = QUAD8B_REF[3];
1274     break;
1275   case 2:
1276     coords[0] = QUAD8B_REF[4];
1277     coords[1] = QUAD8B_REF[5];
1278     break;
1279   case 3:
1280     coords[0] = QUAD8B_REF[6];
1281     coords[1] = QUAD8B_REF[7];
1282     break;
1283   case 4:
1284     coords[0] = QUAD8B_REF[8];
1285     coords[1] = QUAD8B_REF[9];
1286     break;
1287   case 5:
1288     coords[0] = QUAD8B_REF[10];
1289     coords[1] = QUAD8B_REF[11];
1290     break;
1291   case 6:
1292     coords[0] = QUAD8B_REF[12];
1293     coords[1] = QUAD8B_REF[13];
1294     break;
1295   case 7:
1296     coords[0] = QUAD8B_REF[14];
1297     coords[1] = QUAD8B_REF[15];
1298     break;
1299   LOCAL_COORD_MACRO_END;
1300
1301   SHAPE_FUN_MACRO_BEGIN;
1302   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
1303   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
1304   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
1305   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
1306   funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
1307   funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
1308   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
1309   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
1310   SHAPE_FUN_MACRO_END;
1311 }
1312
1313 void GaussInfo::quad9aInit()
1314 {
1315   LOCAL_COORD_MACRO_BEGIN;
1316   case 0:
1317     coords[0] = QUAD9A_REF[0];
1318     coords[1] = QUAD9A_REF[1];
1319     break;
1320   case 1:
1321     coords[0] = QUAD9A_REF[2];
1322     coords[1] = QUAD9A_REF[3];
1323     break;
1324   case 2:
1325     coords[0] = QUAD9A_REF[4];
1326     coords[1] = QUAD9A_REF[5];
1327     break;
1328   case 3:
1329     coords[0] = QUAD9A_REF[6];
1330     coords[1] = QUAD9A_REF[7];
1331     break;
1332   case 4:
1333     coords[0] = QUAD9A_REF[8];
1334     coords[1] = QUAD9A_REF[9];
1335     break;
1336   case 5:
1337     coords[0] = QUAD9A_REF[10];
1338     coords[1] = QUAD9A_REF[11];
1339     break;
1340   case 6:
1341     coords[0] = QUAD9A_REF[12];
1342     coords[1] = QUAD9A_REF[13];
1343     break;
1344   case 7:
1345     coords[0] = QUAD9A_REF[14];
1346     coords[1] = QUAD9A_REF[15];
1347     break;
1348   case 8:
1349     coords[0] = QUAD9A_REF[16];
1350     coords[1] = QUAD9A_REF[17];
1351     break;
1352   LOCAL_COORD_MACRO_END;
1353
1354   SHAPE_FUN_MACRO_BEGIN;
1355   funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
1356   funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
1357   funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
1358   funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
1359   funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
1360   funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
1361   funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
1362   funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
1363   funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
1364   SHAPE_FUN_MACRO_END;
1365 }
1366
1367 /*!
1368  * Init Tetrahedron Reference coordinates ans Shape function.
1369  * Case A.
1370  */
1371 void GaussInfo::tetra4aInit()
1372 {
1373   LOCAL_COORD_MACRO_BEGIN;
1374   case 0:
1375     coords[0] = TETRA4A_REF[0];
1376     coords[1] = TETRA4A_REF[1];
1377     coords[2] = TETRA4A_REF[2];
1378     break;
1379   case 1:
1380     coords[0] = TETRA4A_REF[3];
1381     coords[1] = TETRA4A_REF[4];
1382     coords[2] = TETRA4A_REF[5];
1383     break;
1384   case 2:
1385     coords[0] = TETRA4A_REF[6];
1386     coords[1] = TETRA4A_REF[7];
1387     coords[2] = TETRA4A_REF[8];
1388     break;
1389   case 3:
1390     coords[0] = TETRA4A_REF[9];
1391     coords[1] = TETRA4A_REF[10];
1392     coords[2] = TETRA4A_REF[11];
1393     break;
1394   LOCAL_COORD_MACRO_END;
1395
1396   SHAPE_FUN_MACRO_BEGIN;
1397   funValue[0] = gc[1];
1398   funValue[1] = gc[2];
1399   funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
1400   funValue[3] = gc[0];
1401   SHAPE_FUN_MACRO_END;
1402 }
1403
1404 /*!
1405  * Init Tetrahedron Reference coordinates ans Shape function.
1406  * Case B.
1407  */
1408 void GaussInfo::tetra4bInit()
1409 {
1410   LOCAL_COORD_MACRO_BEGIN;
1411   case 0:
1412     coords[0] = TETRA4B_REF[0];
1413     coords[1] = TETRA4B_REF[1];
1414     coords[2] = TETRA4B_REF[2];
1415     break;
1416   case 1:
1417     coords[0] = TETRA4B_REF[3];
1418     coords[1] = TETRA4B_REF[4];
1419     coords[2] = TETRA4B_REF[5];
1420     break;
1421   case 2:
1422     coords[0] = TETRA4B_REF[6];
1423     coords[1] = TETRA4B_REF[7];
1424     coords[2] = TETRA4B_REF[8];
1425     break;
1426   case 3:
1427     coords[0] = TETRA4B_REF[9];
1428     coords[1] = TETRA4B_REF[10];
1429     coords[2] = TETRA4B_REF[11];
1430     break;
1431   LOCAL_COORD_MACRO_END;
1432
1433   SHAPE_FUN_MACRO_BEGIN;
1434   funValue[0] = gc[1];
1435   funValue[2] = gc[2];
1436   funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
1437   funValue[3] = gc[0];
1438   SHAPE_FUN_MACRO_END;
1439 }
1440
1441 /*!
1442  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1443  * Case A.
1444  */
1445 void GaussInfo::tetra10aInit()
1446 {
1447   LOCAL_COORD_MACRO_BEGIN;
1448   case 0:
1449     coords[0] = TETRA10A_REF[0];
1450     coords[1] = TETRA10A_REF[1];
1451     coords[2] = TETRA10A_REF[2];
1452     break;
1453   case 1:
1454     coords[0] = TETRA10A_REF[3];
1455     coords[1] = TETRA10A_REF[4];
1456     coords[2] = TETRA10A_REF[5];
1457     break;
1458   case 2:
1459     coords[0] = TETRA10A_REF[6];
1460     coords[1] = TETRA10A_REF[7];
1461     coords[2] = TETRA10A_REF[8];
1462     break;
1463   case 3:
1464     coords[0] = TETRA10A_REF[9];
1465     coords[1] = TETRA10A_REF[10];
1466     coords[2] = TETRA10A_REF[11];
1467     break;
1468   case 4:
1469     coords[0] = TETRA10A_REF[12];
1470     coords[1] = TETRA10A_REF[13];
1471     coords[2] = TETRA10A_REF[14];
1472     break;
1473   case 5:
1474     coords[0] = TETRA10A_REF[15];
1475     coords[1] = TETRA10A_REF[16];
1476     coords[2] = TETRA10A_REF[17];
1477     break;
1478   case 6:
1479     coords[0] = TETRA10A_REF[18];
1480     coords[1] = TETRA10A_REF[19];
1481     coords[2] = TETRA10A_REF[20];
1482     break;
1483   case 7:
1484     coords[0] = TETRA10A_REF[21];
1485     coords[1] = TETRA10A_REF[22];
1486     coords[2] = TETRA10A_REF[23];
1487     break;
1488   case 8:
1489     coords[0] = TETRA10A_REF[24];
1490     coords[1] = TETRA10A_REF[25];
1491     coords[2] = TETRA10A_REF[26];
1492     break;
1493   case 9:
1494     coords[0] = TETRA10A_REF[27];
1495     coords[1] = TETRA10A_REF[28];
1496     coords[2] = TETRA10A_REF[29];
1497     break;
1498   LOCAL_COORD_MACRO_END;
1499
1500   SHAPE_FUN_MACRO_BEGIN;
1501   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1502   funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
1503   funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1504   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1505   funValue[4] = 4.0*gc[1]*gc[2];
1506   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1507   funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1508   funValue[7] = 4.0*gc[0]*gc[1];
1509   funValue[8] = 4.0*gc[0]*gc[2];
1510   funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1511   SHAPE_FUN_MACRO_END;
1512 }
1513
1514 /*!
1515  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1516  * Case B.
1517  */
1518 void GaussInfo::tetra10bInit()
1519 {
1520   LOCAL_COORD_MACRO_BEGIN;
1521   case 0:
1522     coords[0] = TETRA10B_REF[0];
1523     coords[1] = TETRA10B_REF[1];
1524     coords[2] = TETRA10B_REF[2];
1525     break;
1526   case 1:
1527     coords[0] = TETRA10B_REF[3];
1528     coords[1] = TETRA10B_REF[4];
1529     coords[2] = TETRA10B_REF[5];
1530     break;
1531   case 2:
1532     coords[0] = TETRA10B_REF[6];
1533     coords[1] = TETRA10B_REF[7];
1534     coords[2] = TETRA10B_REF[8];
1535     break;
1536   case 3:
1537     coords[0] = TETRA10B_REF[9];
1538     coords[1] = TETRA10B_REF[10];
1539     coords[2] = TETRA10B_REF[11];
1540     break;
1541   case 4:
1542     coords[0] = TETRA10B_REF[12];
1543     coords[1] = TETRA10B_REF[13];
1544     coords[2] = TETRA10B_REF[14];
1545     break;
1546   case 5:
1547     coords[0] = TETRA10B_REF[15];
1548     coords[1] = TETRA10B_REF[16];
1549     coords[2] = TETRA10B_REF[17];
1550     break;
1551   case 6:
1552     coords[0] = TETRA10B_REF[18];
1553     coords[1] = TETRA10B_REF[19];
1554     coords[2] = TETRA10B_REF[20];
1555     break;
1556   case 7:
1557     coords[0] = TETRA10B_REF[21];
1558     coords[1] = TETRA10B_REF[22];
1559     coords[2] = TETRA10B_REF[23];
1560     break;
1561   case 8:
1562     coords[0] = TETRA10B_REF[24];
1563     coords[1] = TETRA10B_REF[25];
1564     coords[2] = TETRA10B_REF[26];
1565     break;
1566   case 9:
1567     coords[0] = TETRA10B_REF[27];
1568     coords[1] = TETRA10B_REF[28];
1569     coords[2] = TETRA10B_REF[29];
1570     break;
1571   LOCAL_COORD_MACRO_END;
1572   SHAPE_FUN_MACRO_BEGIN;
1573   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1574   funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1575   funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1576   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1577   funValue[6] = 4.0*gc[1]*gc[2];
1578   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1579   funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1580   funValue[7] = 4.0*gc[0]*gc[1];
1581   funValue[9] = 4.0*gc[0]*gc[2];
1582   funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1583   SHAPE_FUN_MACRO_END;
1584 }
1585
1586 /*!
1587  * Init Pyramid Reference coordinates ans Shape function.
1588  * Case A.
1589  */
1590 void GaussInfo::pyra5aInit()
1591 {
1592   LOCAL_COORD_MACRO_BEGIN;
1593   case 0:
1594     coords[0] = PYRA5A_REF[0];
1595     coords[1] = PYRA5A_REF[1];
1596     coords[2] = PYRA5A_REF[2];
1597     break;
1598   case 1:
1599     coords[0] = PYRA5A_REF[3];
1600     coords[1] = PYRA5A_REF[4];
1601     coords[2] = PYRA5A_REF[5];
1602     break;
1603   case 2:
1604     coords[0] = PYRA5A_REF[6];
1605     coords[1] = PYRA5A_REF[7];
1606     coords[2] = PYRA5A_REF[8];
1607     break;
1608   case 3:
1609     coords[0] = PYRA5A_REF[9];
1610     coords[1] = PYRA5A_REF[10];
1611     coords[2] = PYRA5A_REF[11];
1612     break;
1613   case 4:
1614     coords[0] = PYRA5A_REF[12];
1615     coords[1] = PYRA5A_REF[13];
1616     coords[2] = PYRA5A_REF[14];
1617     break;
1618   LOCAL_COORD_MACRO_END;
1619
1620   SHAPE_FUN_MACRO_BEGIN;
1621   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1622   funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1623   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1624   funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1625   funValue[4] = gc[2];
1626   SHAPE_FUN_MACRO_END;
1627
1628   DEV_SHAPE_FUN_MACRO_BEGIN;
1629
1630   devFunValue[0] = 0.5*(gc[0]+1)*(1.0 - gc[2]);
1631   devFunValue[1] = -0.5*gc[1]*(1.0 - gc[2]);
1632   devFunValue[2] = -0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0);
1633
1634   devFunValue[3] = -0.5*gc[0]*(1.0 - gc[2]);
1635   devFunValue[4] = 0.5*(gc[1]+1)*(1.0 - gc[2]);
1636   devFunValue[5] = -0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0);
1637
1638   devFunValue[6] = 0.5*(gc[0]-1)*(1.0 - gc[2]);
1639   devFunValue[7] = -0.5*gc[1]*(1.0 - gc[2]);
1640   devFunValue[8] = -0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0);
1641
1642   devFunValue[9] = -0.5*gc[0]*(1.0 - gc[2]);
1643   devFunValue[10] = 0.5*(1*gc[1]-1)*(1.0 - gc[2]);
1644   devFunValue[11] = -0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0);
1645
1646   devFunValue[12] = 0.0;
1647   devFunValue[13] = 0.0;
1648   devFunValue[14] = 1.0;
1649
1650   DEV_SHAPE_FUN_MACRO_END;
1651 }
1652 /*!
1653  * Init Pyramid Reference coordinates ans Shape function.
1654  * Case B.
1655  */
1656 void GaussInfo::pyra5bInit()
1657 {
1658   LOCAL_COORD_MACRO_BEGIN;
1659   case 0:
1660     coords[0] = PYRA5B_REF[0];
1661     coords[1] = PYRA5B_REF[1];
1662     coords[2] = PYRA5B_REF[2];
1663     break;
1664   case 1:
1665     coords[0] = PYRA5B_REF[3];
1666     coords[1] = PYRA5B_REF[4];
1667     coords[2] = PYRA5B_REF[5];
1668     break;
1669   case 2:
1670     coords[0] = PYRA5B_REF[6];
1671     coords[1] = PYRA5B_REF[7];
1672     coords[2] = PYRA5B_REF[8];
1673     break;
1674   case 3:
1675     coords[0] = PYRA5B_REF[9];
1676     coords[1] = PYRA5B_REF[10];
1677     coords[2] = PYRA5B_REF[11];
1678     break;
1679   case 4:
1680     coords[0] = PYRA5B_REF[12];
1681     coords[1] = PYRA5B_REF[13];
1682     coords[2] = PYRA5B_REF[14];
1683     break;
1684   LOCAL_COORD_MACRO_END;
1685
1686   SHAPE_FUN_MACRO_BEGIN;
1687   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1688   funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1689   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1690   funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1691   funValue[4] = gc[2];
1692   SHAPE_FUN_MACRO_END;
1693 }
1694
1695 /*!
1696  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1697  * Case A.
1698  */
1699 void GaussInfo::pyra13aInit()
1700 {
1701   LOCAL_COORD_MACRO_BEGIN;
1702   case 0:
1703     coords[0] = PYRA13A_REF[0];
1704     coords[1] = PYRA13A_REF[1];
1705     coords[2] = PYRA13A_REF[2];
1706     break;
1707   case 1:
1708     coords[0] = PYRA13A_REF[3];
1709     coords[1] = PYRA13A_REF[4];
1710     coords[2] = PYRA13A_REF[5];
1711     break;
1712   case 2:
1713     coords[0] = PYRA13A_REF[6];
1714     coords[1] = PYRA13A_REF[7];
1715     coords[2] = PYRA13A_REF[8];
1716     break;
1717   case 3:
1718     coords[0] = PYRA13A_REF[9];
1719     coords[1] = PYRA13A_REF[10];
1720     coords[2] = PYRA13A_REF[11];
1721     break;
1722   case 4:
1723     coords[0] = PYRA13A_REF[12];
1724     coords[1] = PYRA13A_REF[13];
1725     coords[2] = PYRA13A_REF[14];
1726     break;
1727   case 5:
1728     coords[0] = PYRA13A_REF[15];
1729     coords[1] = PYRA13A_REF[16];
1730     coords[2] = PYRA13A_REF[17];
1731     break;
1732   case 6:
1733     coords[0] = PYRA13A_REF[18];
1734     coords[1] = PYRA13A_REF[19];
1735     coords[2] = PYRA13A_REF[20];
1736     break;
1737   case 7:
1738     coords[0] = PYRA13A_REF[21];
1739     coords[1] = PYRA13A_REF[22];
1740     coords[2] = PYRA13A_REF[23];
1741     break;
1742   case 8:
1743     coords[0] = PYRA13A_REF[24];
1744     coords[1] = PYRA13A_REF[25];
1745     coords[2] = PYRA13A_REF[26];
1746     break;
1747   case 9:
1748     coords[0] = PYRA13A_REF[27];
1749     coords[1] = PYRA13A_REF[28];
1750     coords[2] = PYRA13A_REF[29];
1751     break;
1752   case 10:
1753     coords[0] = PYRA13A_REF[30];
1754     coords[1] = PYRA13A_REF[31];
1755     coords[2] = PYRA13A_REF[32];
1756     break;
1757   case 11:
1758     coords[0] = PYRA13A_REF[33];
1759     coords[1] = PYRA13A_REF[34];
1760     coords[2] = PYRA13A_REF[35];
1761     break;
1762   case 12:
1763     coords[0] = PYRA13A_REF[36];
1764     coords[1] = PYRA13A_REF[37];
1765     coords[2] = PYRA13A_REF[38];
1766     break;
1767   LOCAL_COORD_MACRO_END;
1768
1769   SHAPE_FUN_MACRO_BEGIN;
1770   funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1771     (gc[0] - 0.5)/(1.0 - gc[2]);
1772   funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1773     (gc[1] - 0.5)/(1.0 - gc[2]);
1774   funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1775     (-gc[0] - 0.5)/(1.0 - gc[2]);
1776   funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1777     (-gc[1] - 0.5)/(1.0 - gc[2]);
1778
1779   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1780
1781   funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1782     (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1783   funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1784     (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1785   funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1786     (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1787   funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1788     (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1789
1790   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1791     (1.0 - gc[2]);
1792   funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1793     (1.0 - gc[2]);
1794   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1795     (1.0 - gc[2]);
1796   funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1797     (1.0 - gc[2]);
1798   SHAPE_FUN_MACRO_END;
1799 }
1800
1801 /*!
1802  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1803  * Case B.
1804  */
1805 void GaussInfo::pyra13bInit()
1806 {
1807   LOCAL_COORD_MACRO_BEGIN;
1808   case 0:
1809     coords[0] = PYRA13B_REF[0];
1810     coords[1] = PYRA13B_REF[1];
1811     coords[2] = PYRA13B_REF[2];
1812     break;
1813   case 1:
1814     coords[0] = PYRA13B_REF[3];
1815     coords[1] = PYRA13B_REF[4];
1816     coords[2] = PYRA13B_REF[5];
1817     break;
1818   case 2:
1819     coords[0] = PYRA13B_REF[6];
1820     coords[1] = PYRA13B_REF[7];
1821     coords[2] = PYRA13B_REF[8];
1822     break;
1823   case 3:
1824     coords[0] = PYRA13B_REF[9];
1825     coords[1] = PYRA13B_REF[10];
1826     coords[2] = PYRA13B_REF[11];
1827     break;
1828   case 4:
1829     coords[0] = PYRA13B_REF[12];
1830     coords[1] = PYRA13B_REF[13];
1831     coords[2] = PYRA13B_REF[14];
1832     break;
1833   case 5:
1834     coords[0] = PYRA13B_REF[15];
1835     coords[1] = PYRA13B_REF[16];
1836     coords[2] = PYRA13B_REF[17];
1837     break;
1838   case 6:
1839     coords[0] = PYRA13B_REF[18];
1840     coords[1] = PYRA13B_REF[19];
1841     coords[2] = PYRA13B_REF[20];
1842     break;
1843   case 7:
1844     coords[0] = PYRA13B_REF[21];
1845     coords[1] = PYRA13B_REF[22];
1846     coords[2] = PYRA13B_REF[23];
1847     break;
1848   case 8:
1849     coords[0] = PYRA13B_REF[24];
1850     coords[1] = PYRA13B_REF[25];
1851     coords[2] = PYRA13B_REF[26];
1852     break;
1853   case 9:
1854     coords[0] = PYRA13B_REF[27];
1855     coords[1] = PYRA13B_REF[28];
1856     coords[2] = PYRA13B_REF[29];
1857     break;
1858   case 10:
1859     coords[0] = PYRA13B_REF[30];
1860     coords[1] = PYRA13B_REF[31];
1861     coords[2] = PYRA13B_REF[32];
1862     break;
1863   case 11:
1864     coords[0] = PYRA13B_REF[33];
1865     coords[1] = PYRA13B_REF[34];
1866     coords[2] = PYRA13B_REF[35];
1867     break;
1868   case 12:
1869     coords[0] = PYRA13B_REF[36];
1870     coords[1] = PYRA13B_REF[37];
1871     coords[2] = PYRA13B_REF[38];
1872     break;
1873   LOCAL_COORD_MACRO_END;
1874
1875   SHAPE_FUN_MACRO_BEGIN;
1876   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]);
1877   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]);
1878   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]);
1879   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]);
1880
1881   funValue[4] =2.*gc[2]*(gc[2]-0.5);
1882
1883   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]);
1884   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]);
1885   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]);
1886   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]);
1887
1888   funValue[9] =gc[2]*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1889   funValue[10]=gc[2]*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1890   funValue[11]=gc[2]*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1891   funValue[12]=gc[2]*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
1892
1893   SHAPE_FUN_MACRO_END;
1894 }
1895
1896
1897 /*!
1898  * Init Pentahedron Reference coordinates and Shape function.
1899  * Case A.
1900  */
1901 void GaussInfo::penta6aInit()
1902 {
1903   LOCAL_COORD_MACRO_BEGIN;
1904   case 0:
1905     coords[0] = PENTA6A_REF[0];
1906     coords[1] = PENTA6A_REF[1];
1907     coords[2] = PENTA6A_REF[2];
1908     break;
1909   case 1:
1910     coords[0] = PENTA6A_REF[3];
1911     coords[1] = PENTA6A_REF[4];
1912     coords[2] = PENTA6A_REF[5];
1913     break;
1914   case 2:
1915     coords[0] = PENTA6A_REF[6];
1916     coords[1] = PENTA6A_REF[7];
1917     coords[2] = PENTA6A_REF[8];
1918     break;
1919   case 3:
1920     coords[0] = PENTA6A_REF[9];
1921     coords[1] = PENTA6A_REF[10];
1922     coords[2] = PENTA6A_REF[11];
1923     break;
1924   case 4:
1925     coords[0] = PENTA6A_REF[12];
1926     coords[1] = PENTA6A_REF[13];
1927     coords[2] = PENTA6A_REF[14];
1928     break;
1929   case 5:
1930     coords[0] = PENTA6A_REF[15];
1931     coords[1] = PENTA6A_REF[16];
1932     coords[2] = PENTA6A_REF[17];
1933     break;
1934   LOCAL_COORD_MACRO_END;
1935
1936   SHAPE_FUN_MACRO_BEGIN;
1937   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1938   funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
1939   funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1940
1941   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1942   funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
1943   funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1944   SHAPE_FUN_MACRO_END;
1945
1946   DEV_SHAPE_FUN_MACRO_BEGIN;
1947
1948   devFunValue[0] = 0.5*gc[1]*(-1.0);
1949   devFunValue[1] = 0.5*(1.0 - gc[0]);
1950   devFunValue[2] = 0.0;
1951
1952   devFunValue[3] = 0.5*gc[2]*(-1.0);
1953   devFunValue[4] = 0.0;
1954   devFunValue[5] = 0.5*(1.0 - gc[0]);
1955
1956   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(-1.0);
1957   devFunValue[7] = 0.5*(-1.0)*(1.0 - gc[0]);
1958   devFunValue[8] = 0.5*(-1.0)*(1.0 - gc[0]);
1959
1960   devFunValue[9] = 0.5*gc[1];
1961   devFunValue[10] = 0.5*(gc[0] + 1.0);
1962   devFunValue[11] = 0.0;
1963
1964   devFunValue[12] = 0.5*gc[2];
1965   devFunValue[13] = 0.0;
1966   devFunValue[14] = 0.5*(gc[0] + 1.0);
1967
1968   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2]);
1969   devFunValue[16] = 0.5*(-1.0)*(1.0 + gc[0]);
1970   devFunValue[17] = 0.5*(-1.0)*(1.0 + gc[0]);
1971
1972   DEV_SHAPE_FUN_MACRO_END;
1973 }
1974
1975 /*!
1976  * Init Pentahedron Reference coordinates and Shape function.
1977  * Case B.
1978  */
1979 void GaussInfo::penta6bInit()
1980 {
1981   LOCAL_COORD_MACRO_BEGIN;
1982   case 0:
1983     coords[0] = PENTA6B_REF[0];
1984     coords[1] = PENTA6B_REF[1];
1985     coords[2] = PENTA6B_REF[2];
1986     break;
1987   case 1:
1988     coords[0] = PENTA6B_REF[3];
1989     coords[1] = PENTA6B_REF[4];
1990     coords[2] = PENTA6B_REF[5];
1991     break;
1992   case 2:
1993     coords[0] = PENTA6B_REF[6];
1994     coords[1] = PENTA6B_REF[7];
1995     coords[2] = PENTA6B_REF[8];
1996     break;
1997   case 3:
1998     coords[0] = PENTA6B_REF[9];
1999     coords[1] = PENTA6B_REF[10];
2000     coords[2] = PENTA6B_REF[11];
2001     break;
2002   case 4:
2003     coords[0] = PENTA6B_REF[12];
2004     coords[1] = PENTA6B_REF[13];
2005     coords[2] = PENTA6B_REF[14];
2006     break;
2007   case 5:
2008     coords[0] = PENTA6B_REF[15];
2009     coords[1] = PENTA6B_REF[16];
2010     coords[2] = PENTA6B_REF[17];
2011     break;
2012   LOCAL_COORD_MACRO_END;
2013
2014   SHAPE_FUN_MACRO_BEGIN;
2015   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
2016   funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
2017   funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2018   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
2019   funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
2020   funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2021   SHAPE_FUN_MACRO_END;
2022 }
2023
2024 /*!
2025  * This shapefunc map is same as degenerated tria3aInit
2026  */
2027 void GaussInfo::penta6DegTria3aInit()
2028 {
2029   LOCAL_COORD_MACRO_BEGIN;
2030  case  0:
2031    coords[0] = -1.0;
2032    coords[1] =  1.0;
2033    coords[2] =  0.0;
2034    break;
2035  case  1:
2036    coords[0] = -1.0;
2037    coords[1] = -1.0;
2038    coords[2] =  0.0;
2039    break;
2040  case  2:
2041    coords[0] =  1.0;
2042    coords[1] = -1.0;
2043    coords[2] =  0.0;
2044    break;
2045  case  3:
2046    coords[0] =  0.0;
2047    coords[1] =  0.0;
2048    coords[2] =  0.0;
2049    break;
2050  case  4:
2051    coords[0] =  0.0;
2052    coords[1] =  0.0;
2053    coords[2] =  0.0;
2054    break;
2055  case  5:
2056    coords[0] =  0.0;
2057    coords[1] =  0.0;
2058    coords[2] =  0.0;
2059    break;
2060    LOCAL_COORD_MACRO_END;
2061
2062    SHAPE_FUN_MACRO_BEGIN;
2063    funValue[0] = 0.5*(1.0 + gc[1]);
2064    funValue[1] = -0.5*(gc[0] + gc[1]);
2065    funValue[2] = 0.5*(1.0 + gc[0]);
2066    funValue[3] = 0.;
2067    funValue[4] = 0.;
2068    funValue[5] = 0.;
2069    SHAPE_FUN_MACRO_END;
2070 }
2071
2072 /*!
2073  * This shapefunc map is same as degenerated tria3bInit
2074  */
2075 void GaussInfo::penta6DegTria3bInit()
2076 {
2077   LOCAL_COORD_MACRO_BEGIN;
2078  case  0:
2079    coords[0] =  0.0;
2080    coords[1] =  0.0;
2081    coords[2] =  0.0;
2082    break;
2083  case  1:
2084    coords[0] =  1.0;
2085    coords[1] =  0.0;
2086    coords[2] =  0.0;
2087    break;
2088  case  2:
2089    coords[0] =  0.0;
2090    coords[1] =  1.0;
2091    coords[2] =  0.0;
2092    break;
2093  case  3:
2094    coords[0] =  0.0;
2095    coords[1] =  0.0;
2096    coords[2] =  0.0;
2097    break;
2098  case  4:
2099    coords[0] =  0.0;
2100    coords[1] =  0.0;
2101    coords[2] =  0.0;
2102    break;
2103  case  5:
2104    coords[0] =  0.0;
2105    coords[1] =  0.0;
2106    coords[2] =  0.0;
2107    break;
2108    LOCAL_COORD_MACRO_END;
2109
2110    SHAPE_FUN_MACRO_BEGIN;
2111    funValue[0] = 1.0 - gc[0] - gc[1];
2112    funValue[1] = gc[0];
2113    funValue[2] = gc[1];
2114    funValue[3] = 0.;
2115    funValue[4] = 0.;
2116    funValue[5] = 0.;
2117    SHAPE_FUN_MACRO_END;
2118 }
2119
2120 /*!
2121  * Init Pentahedron Reference coordinates and Shape function.
2122  * Case A.
2123  */
2124 void GaussInfo::penta15aInit()
2125 {
2126   LOCAL_COORD_MACRO_BEGIN;
2127   case 0:
2128     coords[0] = PENTA15A_REF[0];
2129     coords[1] = PENTA15A_REF[1];
2130     coords[2] = PENTA15A_REF[2];
2131     break;
2132   case 1:
2133     coords[0] = PENTA15A_REF[3];
2134     coords[1] = PENTA15A_REF[4];
2135     coords[2] = PENTA15A_REF[5];
2136     break;
2137   case 2:
2138     coords[0] = PENTA15A_REF[6];
2139     coords[1] = PENTA15A_REF[7];
2140     coords[2] = PENTA15A_REF[8];
2141     break;
2142   case 3:
2143     coords[0] = PENTA15A_REF[9];
2144     coords[1] = PENTA15A_REF[10];
2145     coords[2] = PENTA15A_REF[11];
2146     break;
2147   case 4:
2148     coords[0] = PENTA15A_REF[12];
2149     coords[1] = PENTA15A_REF[13];
2150     coords[2] = PENTA15A_REF[14];
2151     break;
2152   case 5:
2153     coords[0] = PENTA15A_REF[15];
2154     coords[1] = PENTA15A_REF[16];
2155     coords[2] = PENTA15A_REF[17];
2156     break;
2157   case 6:
2158     coords[0] = PENTA15A_REF[18];
2159     coords[1] = PENTA15A_REF[19];
2160     coords[2] = PENTA15A_REF[20];
2161     break;
2162   case 7:
2163     coords[0] = PENTA15A_REF[21];
2164     coords[1] = PENTA15A_REF[22];
2165     coords[2] = PENTA15A_REF[23];
2166     break;
2167   case 8:
2168     coords[0] = PENTA15A_REF[24];
2169     coords[1] = PENTA15A_REF[25];
2170     coords[2] = PENTA15A_REF[26];
2171     break;
2172   case 9:
2173     coords[0] = PENTA15A_REF[27];
2174     coords[1] = PENTA15A_REF[28];
2175     coords[2] = PENTA15A_REF[29];
2176     break;
2177   case 10:
2178     coords[0] = PENTA15A_REF[30];
2179     coords[1] = PENTA15A_REF[31];
2180     coords[2] = PENTA15A_REF[32];
2181     break;
2182   case 11:
2183     coords[0] = PENTA15A_REF[33];
2184     coords[1] = PENTA15A_REF[34];
2185     coords[2] = PENTA15A_REF[35];
2186     break;
2187   case 12:
2188     coords[0] = PENTA15A_REF[36];
2189     coords[1] = PENTA15A_REF[37];
2190     coords[2] = PENTA15A_REF[38];
2191     break;
2192   case 13:
2193     coords[0] = PENTA15A_REF[39];
2194     coords[1] = PENTA15A_REF[40];
2195     coords[2] = PENTA15A_REF[41];
2196     break;
2197   case 14:
2198     coords[0] = PENTA15A_REF[42];
2199     coords[1] = PENTA15A_REF[43];
2200     coords[2] = PENTA15A_REF[44];
2201     break;
2202   LOCAL_COORD_MACRO_END;
2203
2204   SHAPE_FUN_MACRO_BEGIN;
2205   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2206   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2207   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2208
2209   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2210   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2211   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2212
2213   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2214   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2215   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2216
2217   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
2218   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
2219   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2220
2221   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2222   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2223   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2224   SHAPE_FUN_MACRO_END;
2225
2226   DEV_SHAPE_FUN_MACRO_BEGIN;
2227
2228   devFunValue[0] = 0.5*gc[1]*(2 * gc[0] - 2 * gc[1] + 1.0 );
2229   devFunValue[1] = (1.0 - gc[0])*(2.0*gc[1] -1 - 0.5*gc[0]);
2230   devFunValue[2] = 0.0;
2231
2232   devFunValue[3] = 0.5*gc[2]*(2 * gc[0] - 2 * gc[2] + 1.0 );
2233   devFunValue[4] = 0.0;
2234   devFunValue[5] = (1.0 - gc[0])*(2.0*gc[2] -1 - 0.5*gc[0]);
2235
2236   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(2*gc[0] -1.0 + 2*gc[1] + 2*gc[2]);
2237   devFunValue[7] = 0.5*(gc[0] - 1.0)* ( -4*gc[1] -gc[0] -4*gc[2] + 2.0);
2238   devFunValue[8] = 0.5*(gc[0] - 1.0)* ( -4*gc[2] -4*gc[1] -gc[0] + 2.0);
2239
2240   devFunValue[9] = 0.5*gc[1]*( 2*gc[0] + 2*gc[1] -1.0 );
2241   devFunValue[10] = (1.0 + gc[0])*(2.0*gc[1] - 1.0 + 0.5*gc[0]);
2242   devFunValue[11] = 0.0;
2243
2244   devFunValue[12] = 0.5*gc[2]*( 2*gc[0] + 2*gc[2] -1.0 );
2245   devFunValue[13] = 0.0;
2246   devFunValue[14] = (1.0 + gc[0])*(2.0*gc[2] - 1.0 + 0.5*gc[0]);
2247
2248   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2])*(2 * gc[0] - 2.0 * gc[1] - 2.0 * gc[2] + 1.0 );
2249   devFunValue[16] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[1] + gc[0] - 4.0 * gc[2] + 2.0) ;
2250   devFunValue[17] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[2] - 4*gc[1] + gc[0] + 2.0 );
2251
2252   devFunValue[18] = - 2.0*gc[1]*gc[2];
2253   devFunValue[19] = 2.0*gc[2]*(1.0 - gc[0]);
2254   devFunValue[20] = 2.0*gc[1]*(1.0 - gc[0]);
2255
2256   devFunValue[21] = -2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2257   devFunValue[22] = -2.0*gc[2]*(1.0 - gc[0]);
2258   devFunValue[23] = (-4*gc[2]-2*gc[1]+2.0)*(1.0 - gc[0]);
2259
2260   devFunValue[24] = -2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2261   devFunValue[25] = (-4*gc[1]-2*gc[2]+2)*(1.0 - gc[0]);
2262   devFunValue[26] = -2.0*gc[1]*(1.0 - gc[0]);
2263
2264   devFunValue[27] = gc[1]*(- 2.0 * gc[0]);
2265   devFunValue[28] = (1.0 - gc[0]*gc[0]);
2266   devFunValue[29] = 0.0;
2267
2268   devFunValue[30] = -2.0 * gc[2] *gc[0];
2269   devFunValue[31] = 0.0;
2270   devFunValue[32] = (1.0 - gc[0]*gc[0]);
2271
2272   devFunValue[33] = -2.0 * (1.0 - gc[1] - gc[2]) * gc[0];
2273   devFunValue[34] = -1.0*(1.0 - gc[0]*gc[0]);
2274   devFunValue[35] = -1.0*(1.0 - gc[0]*gc[0]);
2275
2276   devFunValue[36] = 2.0*gc[1]*gc[2];
2277   devFunValue[37] = 2.0*gc[2]*(1.0 + gc[0]);
2278   devFunValue[38] = 2.0*gc[1]*(1.0 + gc[0]);
2279
2280   devFunValue[39] = 2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2281   devFunValue[40] = -2.0*gc[2]*(1.0 + gc[0]);
2282   devFunValue[41] = (2.0 - 2.0 * gc[1] - 4.0 * gc[2])*(1.0 + gc[0]);
2283
2284   devFunValue[42] = 2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2285   devFunValue[43] = (2.0 - 4*gc[1] - 2*gc[2])*(1.0 + gc[0]);
2286   devFunValue[44] = -2.0*gc[1]*(1.0 + gc[0]);
2287
2288   DEV_SHAPE_FUN_MACRO_END;
2289 }
2290
2291 /*!
2292  * Init Qaudratic Pentahedron Reference coordinates and Shape function.
2293  * Case B.
2294  */
2295 void GaussInfo::penta15bInit()
2296 {
2297   LOCAL_COORD_MACRO_BEGIN;
2298   case 0:
2299     coords[0] = PENTA15B_REF[0];
2300     coords[1] = PENTA15B_REF[1];
2301     coords[2] = PENTA15B_REF[2];
2302     break;
2303   case 1:
2304     coords[0] = PENTA15B_REF[3];
2305     coords[1] = PENTA15B_REF[4];
2306     coords[2] = PENTA15B_REF[5];
2307     break;
2308   case 2:
2309     coords[0] = PENTA15B_REF[6];
2310     coords[1] = PENTA15B_REF[7];
2311     coords[2] = PENTA15B_REF[8];
2312     break;
2313   case 3:
2314     coords[0] = PENTA15B_REF[9];
2315     coords[1] = PENTA15B_REF[10];
2316     coords[2] = PENTA15B_REF[11];
2317     break;
2318   case 4:
2319     coords[0] = PENTA15B_REF[12];
2320     coords[1] = PENTA15B_REF[13];
2321     coords[2] = PENTA15B_REF[14];
2322     break;
2323   case 5:
2324     coords[0] = PENTA15B_REF[15];
2325     coords[1] = PENTA15B_REF[16];
2326     coords[2] = PENTA15B_REF[17];
2327     break;
2328   case 6:
2329     coords[0] = PENTA15B_REF[18];
2330     coords[1] = PENTA15B_REF[19];
2331     coords[2] = PENTA15B_REF[20];
2332     break;
2333   case 7:
2334     coords[0] = PENTA15B_REF[21];
2335     coords[1] = PENTA15B_REF[22];
2336     coords[2] = PENTA15B_REF[23];
2337     break;
2338   case 8:
2339     coords[0] = PENTA15B_REF[24];
2340     coords[1] = PENTA15B_REF[25];
2341     coords[2] = PENTA15B_REF[26];
2342     break;
2343   case 9:
2344     coords[0] = PENTA15B_REF[27];
2345     coords[1] = PENTA15B_REF[28];
2346     coords[2] = PENTA15B_REF[29];
2347     break;
2348   case 10:
2349     coords[0] = PENTA15B_REF[30];
2350     coords[1] = PENTA15B_REF[31];
2351     coords[2] = PENTA15B_REF[32];
2352     break;
2353   case 11:
2354     coords[0] = PENTA15B_REF[33];
2355     coords[1] = PENTA15B_REF[34];
2356     coords[2] = PENTA15B_REF[35];
2357     break;
2358   case 12:
2359     coords[0] = PENTA15B_REF[36];
2360     coords[1] = PENTA15B_REF[37];
2361     coords[2] = PENTA15B_REF[38];
2362     break;
2363   case 13:
2364     coords[0] = PENTA15B_REF[39];
2365     coords[1] = PENTA15B_REF[40];
2366     coords[2] = PENTA15B_REF[41];
2367     break;
2368   case 14:
2369     coords[0] = PENTA15B_REF[42];
2370     coords[1] = PENTA15B_REF[43];
2371     coords[2] = PENTA15B_REF[44];
2372     break;
2373   LOCAL_COORD_MACRO_END;
2374
2375   SHAPE_FUN_MACRO_BEGIN;
2376   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2377   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2378   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2379
2380   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2381   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2382   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2383
2384   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2385   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2386   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2387
2388   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
2389   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
2390   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2391
2392   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2393   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2394   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2395   SHAPE_FUN_MACRO_END;
2396 }
2397
2398 void GaussInfo::penta18aInit()
2399 {
2400   LOCAL_COORD_MACRO_BEGIN;
2401   case 0:
2402     coords[0] = PENTA18A_REF[0];
2403     coords[1] = PENTA18A_REF[1];
2404     coords[2] = PENTA18A_REF[2];
2405     break;
2406   case 1:
2407     coords[0] = PENTA18A_REF[3];
2408     coords[1] = PENTA18A_REF[4];
2409     coords[2] = PENTA18A_REF[5];
2410     break;
2411   case 2:
2412     coords[0] = PENTA18A_REF[6];
2413     coords[1] = PENTA18A_REF[7];
2414     coords[2] = PENTA18A_REF[8];
2415     break;
2416   case 3:
2417     coords[0] = PENTA18A_REF[9];
2418     coords[1] = PENTA18A_REF[10];
2419     coords[2] = PENTA18A_REF[11];
2420     break;
2421   case 4:
2422     coords[0] = PENTA18A_REF[12];
2423     coords[1] = PENTA18A_REF[13];
2424     coords[2] = PENTA18A_REF[14];
2425     break;
2426   case 5:
2427     coords[0] = PENTA18A_REF[15];
2428     coords[1] = PENTA18A_REF[16];
2429     coords[2] = PENTA18A_REF[17];
2430     break;
2431   case 6:
2432     coords[0] = PENTA18A_REF[18];
2433     coords[1] = PENTA18A_REF[19];
2434     coords[2] = PENTA18A_REF[20];
2435     break;
2436   case 7:
2437     coords[0] = PENTA18A_REF[21];
2438     coords[1] = PENTA18A_REF[22];
2439     coords[2] = PENTA18A_REF[23];
2440     break;
2441   case 8:
2442     coords[0] = PENTA18A_REF[24];
2443     coords[1] = PENTA18A_REF[25];
2444     coords[2] = PENTA18A_REF[26];
2445     break;
2446   case 9:
2447     coords[0] = PENTA18A_REF[27];
2448     coords[1] = PENTA18A_REF[28];
2449     coords[2] = PENTA18A_REF[29];
2450     break;
2451   case 10:
2452     coords[0] = PENTA18A_REF[30];
2453     coords[1] = PENTA18A_REF[31];
2454     coords[2] = PENTA18A_REF[32];
2455     break;
2456   case 11:
2457     coords[0] = PENTA18A_REF[33];
2458     coords[1] = PENTA18A_REF[34];
2459     coords[2] = PENTA18A_REF[35];
2460     break;
2461   case 12:
2462     coords[0] = PENTA18A_REF[36];
2463     coords[1] = PENTA18A_REF[37];
2464     coords[2] = PENTA18A_REF[38];
2465     break;
2466   case 13:
2467     coords[0] = PENTA18A_REF[39];
2468     coords[1] = PENTA18A_REF[40];
2469     coords[2] = PENTA18A_REF[41];
2470     break;
2471   case 14:
2472     coords[0] = PENTA18A_REF[42];
2473     coords[1] = PENTA18A_REF[43];
2474     coords[2] = PENTA18A_REF[44];
2475     break;
2476   case 15:
2477     coords[0] = PENTA18A_REF[45];
2478     coords[1] = PENTA18A_REF[46];
2479     coords[2] = PENTA18A_REF[47];
2480     break;
2481   case 16:
2482     coords[0] = PENTA18A_REF[48];
2483     coords[1] = PENTA18A_REF[49];
2484     coords[2] = PENTA18A_REF[50];
2485     break;
2486   case 17:
2487     coords[0] = PENTA18A_REF[51];
2488     coords[1] = PENTA18A_REF[52];
2489     coords[2] = PENTA18A_REF[53];
2490     break;
2491   LOCAL_COORD_MACRO_END;
2492
2493   SHAPE_FUN_MACRO_BEGIN;
2494   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2495   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2496   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2497
2498   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2499   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2500   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2501
2502   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2503   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2504   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2505
2506   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
2507   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
2508   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2509
2510   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2511   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2512   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2513
2514   funValue[15] = 4.0*gc[1]*gc[2]*(1.0 - gc[0] * gc[0]);
2515   funValue[16] = 4.0*gc[2]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2516   funValue[17] = 4.0*gc[1]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2517   SHAPE_FUN_MACRO_END;
2518 }
2519
2520 void GaussInfo::penta18bInit()
2521 {
2522   LOCAL_COORD_MACRO_BEGIN;
2523   case 0:
2524     coords[0] = PENTA18B_REF[0];
2525     coords[1] = PENTA18B_REF[1];
2526     coords[2] = PENTA18B_REF[2];
2527     break;
2528   case 1:
2529     coords[0] = PENTA18B_REF[3];
2530     coords[1] = PENTA18B_REF[4];
2531     coords[2] = PENTA18B_REF[5];
2532     break;
2533   case 2:
2534     coords[0] = PENTA18B_REF[6];
2535     coords[1] = PENTA18B_REF[7];
2536     coords[2] = PENTA18B_REF[8];
2537     break;
2538   case 3:
2539     coords[0] = PENTA18B_REF[9];
2540     coords[1] = PENTA18B_REF[10];
2541     coords[2] = PENTA18B_REF[11];
2542     break;
2543   case 4:
2544     coords[0] = PENTA18B_REF[12];
2545     coords[1] = PENTA18B_REF[13];
2546     coords[2] = PENTA18B_REF[14];
2547     break;
2548   case 5:
2549     coords[0] = PENTA18B_REF[15];
2550     coords[1] = PENTA18B_REF[16];
2551     coords[2] = PENTA18B_REF[17];
2552     break;
2553   case 6:
2554     coords[0] = PENTA18B_REF[18];
2555     coords[1] = PENTA18B_REF[19];
2556     coords[2] = PENTA18B_REF[20];
2557     break;
2558   case 7:
2559     coords[0] = PENTA18B_REF[21];
2560     coords[1] = PENTA18B_REF[22];
2561     coords[2] = PENTA18B_REF[23];
2562     break;
2563   case 8:
2564     coords[0] = PENTA18B_REF[24];
2565     coords[1] = PENTA18B_REF[25];
2566     coords[2] = PENTA18B_REF[26];
2567     break;
2568   case 9:
2569     coords[0] = PENTA18B_REF[27];
2570     coords[1] = PENTA18B_REF[28];
2571     coords[2] = PENTA18B_REF[29];
2572     break;
2573   case 10:
2574     coords[0] = PENTA18B_REF[30];
2575     coords[1] = PENTA18B_REF[31];
2576     coords[2] = PENTA18B_REF[32];
2577     break;
2578   case 11:
2579     coords[0] = PENTA18B_REF[33];
2580     coords[1] = PENTA18B_REF[34];
2581     coords[2] = PENTA18B_REF[35];
2582     break;
2583   case 12:
2584     coords[0] = PENTA18B_REF[36];
2585     coords[1] = PENTA18B_REF[37];
2586     coords[2] = PENTA18B_REF[38];
2587     break;
2588   case 13:
2589     coords[0] = PENTA18B_REF[39];
2590     coords[1] = PENTA18B_REF[40];
2591     coords[2] = PENTA18B_REF[41];
2592     break;
2593   case 14:
2594     coords[0] = PENTA18B_REF[42];
2595     coords[1] = PENTA18B_REF[43];
2596     coords[2] = PENTA18B_REF[44];
2597     break;
2598   case 15:
2599     coords[0] = PENTA18B_REF[45];
2600     coords[1] = PENTA18B_REF[46];
2601     coords[2] = PENTA18B_REF[47];
2602     break;
2603   case 16:
2604     coords[0] = PENTA18B_REF[48];
2605     coords[1] = PENTA18B_REF[49];
2606     coords[2] = PENTA18B_REF[50];
2607     break;
2608   case 17:
2609     coords[0] = PENTA18B_REF[51];
2610     coords[1] = PENTA18B_REF[52];
2611     coords[2] = PENTA18B_REF[53];
2612     break;
2613   LOCAL_COORD_MACRO_END;
2614
2615   SHAPE_FUN_MACRO_BEGIN;
2616   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2617   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2618   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2619
2620   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2621   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2622   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2623
2624   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2625   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2626   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2627
2628   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
2629   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
2630   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2631
2632   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2633   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2634   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2635
2636   funValue[17] = 4.0*gc[1]*gc[2]*(1.0 - gc[0] * gc[0]);
2637   funValue[16] = 4.0*gc[2]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2638   funValue[15] = 4.0*gc[1]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2639   SHAPE_FUN_MACRO_END;
2640 }
2641
2642 /*!
2643  * Init Hehahedron Reference coordinates and Shape function.
2644  * Case A.
2645  */
2646 void GaussInfo::hexa8aInit()
2647 {
2648   LOCAL_COORD_MACRO_BEGIN;
2649   case 0:
2650     coords[0] = HEXA8A_REF[0];
2651     coords[1] = HEXA8A_REF[1];
2652     coords[2] = HEXA8A_REF[2];
2653     break;
2654   case 1:
2655     coords[0] = HEXA8A_REF[3];
2656     coords[1] = HEXA8A_REF[4];
2657     coords[2] = HEXA8A_REF[5];
2658     break;
2659   case 2:
2660     coords[0] = HEXA8A_REF[6];
2661     coords[1] = HEXA8A_REF[7];
2662     coords[2] = HEXA8A_REF[8];
2663     break;
2664   case 3:
2665     coords[0] = HEXA8A_REF[9];
2666     coords[1] = HEXA8A_REF[10];
2667     coords[2] = HEXA8A_REF[11];
2668     break;
2669   case 4:
2670     coords[0] = HEXA8A_REF[12];
2671     coords[1] = HEXA8A_REF[13];
2672     coords[2] = HEXA8A_REF[14];
2673     break;
2674   case 5:
2675     coords[0] = HEXA8A_REF[15];
2676     coords[1] = HEXA8A_REF[16];
2677     coords[2] = HEXA8A_REF[17];
2678     break;
2679   case 6:
2680     coords[0] = HEXA8A_REF[18];
2681     coords[1] = HEXA8A_REF[19];
2682     coords[2] = HEXA8A_REF[20];
2683     break;
2684   case 7:
2685     coords[0] = HEXA8A_REF[21];
2686     coords[1] = HEXA8A_REF[22];
2687     coords[2] = HEXA8A_REF[23];
2688     break;
2689   LOCAL_COORD_MACRO_END;
2690
2691   SHAPE_FUN_MACRO_BEGIN;
2692   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2693   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2694   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2695   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2696
2697   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2698   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2699   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2700   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2701   SHAPE_FUN_MACRO_END;
2702
2703   DEV_SHAPE_FUN_MACRO_BEGIN;
2704
2705   devFunValue[0] = 0.125 * (-1.0) * (1.0 - gc[1])*(1.0 - gc[2]);
2706   devFunValue[1] = 0.125 * (1.0 - gc[0]) * (-1.0) * (1.0 - gc[2]);
2707   devFunValue[2] = 0.125 * (1.0 - gc[0]) * (1.0 - gc[1]) *(-1.0);
2708
2709   devFunValue[3] = 0.125*(1.0 - gc[1])*(1.0 - gc[2]);
2710   devFunValue[4] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 - gc[2]);
2711   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0);
2712
2713   devFunValue[6] = 0.125*(1.0 + gc[1])*(1.0 - gc[2]);
2714   devFunValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[2]);
2715   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0);
2716
2717   devFunValue[9 ] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 - gc[2]);
2718   devFunValue[10] = 0.125*(1.0 - gc[0])*(1.0 - gc[2]);
2719   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0);
2720
2721   devFunValue[12] = 0.125*(-1.0)*(1.0 - gc[1])*(1.0 + gc[2]);
2722   devFunValue[13] = 0.125*(1.0 - gc[0])*(-1.0)*(1.0 + gc[2]);
2723   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1]);
2724
2725   devFunValue[15] = 0.125*(1.0 - gc[1])*(1.0 + gc[2]);
2726   devFunValue[16] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 + gc[2]);
2727   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1]);
2728
2729   devFunValue[18] = 0.125*(1.0 + gc[1])*(1.0 + gc[2]);
2730   devFunValue[19] = 0.125*(1.0 + gc[0])*(1.0 + gc[2]);
2731   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1]);
2732
2733   devFunValue[21] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 + gc[2]);
2734   devFunValue[22] = 0.125*(1.0 - gc[0])*(1.0 + gc[2]);
2735   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1]);
2736
2737   DEV_SHAPE_FUN_MACRO_END;
2738 }
2739
2740 /*!
2741  * Init Hehahedron Reference coordinates and Shape function.
2742  * Case B.
2743  */
2744 void GaussInfo::hexa8bInit()
2745 {
2746   LOCAL_COORD_MACRO_BEGIN;
2747   case 0:
2748     coords[0] = HEXA8B_REF[0];
2749     coords[1] = HEXA8B_REF[1];
2750     coords[2] = HEXA8B_REF[2];
2751     break;
2752   case 1:
2753     coords[0] = HEXA8B_REF[3];
2754     coords[1] = HEXA8B_REF[4];
2755     coords[2] = HEXA8B_REF[5];
2756     break;
2757   case 2:
2758     coords[0] = HEXA8B_REF[6];
2759     coords[1] = HEXA8B_REF[7];
2760     coords[2] = HEXA8B_REF[8];
2761     break;
2762   case 3:
2763     coords[0] = HEXA8B_REF[9];
2764     coords[1] = HEXA8B_REF[10];
2765     coords[2] = HEXA8B_REF[11];
2766     break;
2767   case 4:
2768     coords[0] = HEXA8B_REF[12];
2769     coords[1] = HEXA8B_REF[13];
2770     coords[2] = HEXA8B_REF[14];
2771     break;
2772   case 5:
2773     coords[0] = HEXA8B_REF[15];
2774     coords[1] = HEXA8B_REF[16];
2775     coords[2] = HEXA8B_REF[17];
2776     break;
2777   case 6:
2778     coords[0] = HEXA8B_REF[18];
2779     coords[1] = HEXA8B_REF[19];
2780     coords[2] = HEXA8B_REF[20];
2781     break;
2782   case 7:
2783     coords[0] = HEXA8B_REF[21];
2784     coords[1] = HEXA8B_REF[22];
2785     coords[2] = HEXA8B_REF[23];
2786     break;
2787   LOCAL_COORD_MACRO_END;
2788
2789   SHAPE_FUN_MACRO_BEGIN;
2790   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2791   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2792   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2793   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2794
2795   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2796   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2797   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2798   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2799   SHAPE_FUN_MACRO_END;
2800 }
2801
2802 /*!
2803  * This shapefunc map is same as degenerated quad4bInit
2804  */
2805 void GaussInfo::hexa8DegQuad4aInit()
2806 {
2807   LOCAL_COORD_MACRO_BEGIN;
2808  case  0:
2809    coords[0] = -1.0;
2810    coords[1] =  1.0;
2811    coords[2] =  0.0;
2812    break;
2813  case  1:
2814    coords[0] = -1.0;
2815    coords[1] = -1.0;
2816    coords[2] =  0.0;
2817    break;
2818  case  2:
2819    coords[0] =  1.0;
2820    coords[1] = -1.0;
2821    coords[2] =  0.0;
2822    break;
2823  case  3:
2824    coords[0] =  1.0;
2825    coords[1] =  1.0;
2826    coords[2] =  0.0;
2827    break;
2828  case  4:
2829    coords[0] =  0.0;
2830    coords[1] =  0.0;
2831    coords[2] =  0.0;
2832    break;
2833  case  5:
2834    coords[0] =  0.0;
2835    coords[1] =  0.0;
2836    coords[2] =  0.0;
2837    break;
2838  case  6:
2839    coords[0] =  0.0;
2840    coords[1] =  0.0;
2841    coords[2] =  0.0;
2842    break;
2843  case  7:
2844    coords[0] =  0.0;
2845    coords[1] =  0.0;
2846    coords[2] =  0.0;
2847    break;
2848   LOCAL_COORD_MACRO_END;
2849
2850   SHAPE_FUN_MACRO_BEGIN;
2851   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
2852   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
2853   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
2854   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2855   funValue[4] = 0.;
2856   funValue[5] = 0.;
2857   funValue[6] = 0.;
2858   funValue[7] = 0.;
2859   SHAPE_FUN_MACRO_END;
2860 }
2861
2862 /*!
2863  * This shapefunc map is same as degenerated quad4bInit
2864  */
2865 void GaussInfo::hexa8DegQuad4bInit()
2866 {
2867   LOCAL_COORD_MACRO_BEGIN;
2868  case  0:
2869    coords[0] = -1.0;
2870    coords[1] = -1.0;
2871    coords[2] =  0.0;
2872    break;
2873  case  1:
2874    coords[0] =  1.0;
2875    coords[1] = -1.0;
2876    coords[2] =  0.0;
2877    break;
2878  case  2:
2879    coords[0] =  1.0;
2880    coords[1] =  1.0;
2881    coords[2] =  0.0;
2882    break;
2883  case  3:
2884    coords[0] = -1.0;
2885    coords[1] =  1.0;
2886    coords[2] =  0.0;
2887    break;
2888  case  4:
2889    coords[0] =  0.0;
2890    coords[1] =  0.0;
2891    coords[2] =  0.0;
2892    break;
2893  case  5:
2894    coords[0] =  0.0;
2895    coords[1] =  0.0;
2896    coords[2] =  0.0;
2897    break;
2898  case  6:
2899    coords[0] =  0.0;
2900    coords[1] =  0.0;
2901    coords[2] =  0.0;
2902    break;
2903  case  7:
2904    coords[0] =  0.0;
2905    coords[1] =  0.0;
2906    coords[2] =  0.0;
2907    break;
2908    LOCAL_COORD_MACRO_END;
2909
2910    SHAPE_FUN_MACRO_BEGIN;
2911    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
2912    funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
2913    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2914    funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
2915    funValue[4] = 0.;
2916    funValue[5] = 0.;
2917    funValue[6] = 0.;
2918    funValue[7] = 0.;
2919    SHAPE_FUN_MACRO_END;
2920 }
2921
2922 /*!
2923  * This shapefunc map is same as degenerated quad4cInit
2924  */
2925 void GaussInfo::hexa8DegQuad4cInit()
2926 {
2927   LOCAL_COORD_MACRO_BEGIN;
2928  case  0:
2929    coords[0] = -1.0;
2930    coords[1] = -1.0;
2931    coords[2] =  0.0;
2932    break;
2933  case  1:
2934    coords[0] = -1.0;
2935    coords[1] =  1.0;
2936    coords[2] =  0.0;
2937    break;
2938  case  2:
2939    coords[0] =  1.0;
2940    coords[1] =  1.0;
2941    coords[2] =  0.0;
2942    break;
2943  case  3:
2944    coords[0] =  1.0;
2945    coords[1] = -1.0;
2946    coords[2] =  0.0;
2947    break;
2948  case  4:
2949    coords[0] =  0.0;
2950    coords[1] =  0.0;
2951    coords[2] =  0.0;
2952    break;
2953  case  5:
2954    coords[0] =  0.0;
2955    coords[1] =  0.0;
2956    coords[2] =  0.0;
2957    break;
2958  case  6:
2959    coords[0] =  0.0;
2960    coords[1] =  0.0;
2961    coords[2] =  0.0;
2962    break;
2963  case  7:
2964    coords[0] =  0.0;
2965    coords[1] =  0.0;
2966    coords[2] =  0.0;
2967    break;
2968
2969    LOCAL_COORD_MACRO_END;
2970
2971    SHAPE_FUN_MACRO_BEGIN;
2972    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
2973    funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
2974    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
2975    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
2976    funValue[4] = 0. ;
2977    funValue[5] = 0. ;
2978    funValue[6] = 0. ;
2979    funValue[7] = 0. ;
2980    SHAPE_FUN_MACRO_END;
2981 }
2982
2983 /*!
2984  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
2985  * Case A.
2986  */
2987 void GaussInfo::hexa20aInit()
2988 {
2989   LOCAL_COORD_MACRO_BEGIN;
2990   case 0:
2991     coords[0] = HEXA20A_REF[0];
2992     coords[1] = HEXA20A_REF[1];
2993     coords[2] = HEXA20A_REF[2];
2994     break;
2995   case 1:
2996     coords[0] = HEXA20A_REF[3];
2997     coords[1] = HEXA20A_REF[4];
2998     coords[2] = HEXA20A_REF[5];
2999     break;
3000   case 2:
3001     coords[0] = HEXA20A_REF[6];
3002     coords[1] = HEXA20A_REF[7];
3003     coords[2] = HEXA20A_REF[8];
3004     break;
3005   case 3:
3006     coords[0] = HEXA20A_REF[9];
3007     coords[1] = HEXA20A_REF[10];
3008     coords[2] = HEXA20A_REF[11];
3009     break;
3010   case 4:
3011     coords[0] = HEXA20A_REF[12];
3012     coords[1] = HEXA20A_REF[13];
3013     coords[2] = HEXA20A_REF[14];
3014     break;
3015   case 5:
3016     coords[0] = HEXA20A_REF[15];
3017     coords[1] = HEXA20A_REF[16];
3018     coords[2] = HEXA20A_REF[17];
3019     break;
3020   case 6:
3021     coords[0] = HEXA20A_REF[18];
3022     coords[1] = HEXA20A_REF[19];
3023     coords[2] = HEXA20A_REF[20];
3024     break;
3025   case 7:
3026     coords[0] = HEXA20A_REF[21];
3027     coords[1] = HEXA20A_REF[22];
3028     coords[2] = HEXA20A_REF[23];
3029     break;
3030   case 8:
3031     coords[0] = HEXA20A_REF[24];
3032     coords[1] = HEXA20A_REF[25];
3033     coords[2] = HEXA20A_REF[26];
3034     break;
3035   case 9:
3036     coords[0] = HEXA20A_REF[27];
3037     coords[1] = HEXA20A_REF[28];
3038     coords[2] = HEXA20A_REF[29];
3039     break;
3040   case 10:
3041     coords[0] = HEXA20A_REF[30];
3042     coords[1] = HEXA20A_REF[31];
3043     coords[2] = HEXA20A_REF[32];
3044     break;
3045   case 11:
3046     coords[0] = HEXA20A_REF[33];
3047     coords[1] = HEXA20A_REF[34];
3048     coords[2] = HEXA20A_REF[35];
3049     break;
3050   case 12:
3051     coords[0] = HEXA20A_REF[36];
3052     coords[1] = HEXA20A_REF[37];
3053     coords[2] = HEXA20A_REF[38];
3054     break;
3055   case 13:
3056     coords[0] = HEXA20A_REF[39];
3057     coords[1] = HEXA20A_REF[40];
3058     coords[2] = HEXA20A_REF[41];
3059     break;
3060   case 14:
3061     coords[0] = HEXA20A_REF[42];
3062     coords[1] = HEXA20A_REF[43];
3063     coords[2] = HEXA20A_REF[44];
3064     break;
3065   case 15:
3066     coords[0] = HEXA20A_REF[45];
3067     coords[1] = HEXA20A_REF[46];
3068     coords[2] = HEXA20A_REF[47];
3069     break;
3070   case 16:
3071     coords[0] = HEXA20A_REF[48];
3072     coords[1] = HEXA20A_REF[49];
3073     coords[2] = HEXA20A_REF[50];
3074     break;
3075   case 17:
3076     coords[0] = HEXA20A_REF[51];
3077     coords[1] = HEXA20A_REF[52];
3078     coords[2] = HEXA20A_REF[53];
3079     break;
3080   case 18:
3081     coords[0] = HEXA20A_REF[54];
3082     coords[1] = HEXA20A_REF[55];
3083     coords[2] = HEXA20A_REF[56];
3084     break;
3085   case 19:
3086     coords[0] = HEXA20A_REF[57];
3087     coords[1] = HEXA20A_REF[58];
3088     coords[2] = HEXA20A_REF[59];
3089     break;
3090   LOCAL_COORD_MACRO_END;
3091
3092   SHAPE_FUN_MACRO_BEGIN;
3093   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3094     (-2.0 - gc[0] - gc[1] - gc[2]);
3095   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3096     (-2.0 + gc[0] - gc[1] - gc[2]);
3097   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3098     (-2.0 + gc[0] + gc[1] - gc[2]);
3099   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3100     (-2.0 - gc[0] + gc[1] - gc[2]);
3101   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3102     (-2.0 - gc[0] - gc[1] + gc[2]);
3103   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3104     (-2.0 + gc[0] - gc[1] + gc[2]);
3105   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3106     (-2.0 + gc[0] + gc[1] + gc[2]);
3107   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3108     (-2.0 - gc[0] + gc[1] + gc[2]);
3109
3110   funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3111   funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3112   funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3113   funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3114   funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3115   funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3116   funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3117   funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3118   funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3119   funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3120   funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3121   funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3122   SHAPE_FUN_MACRO_END;
3123
3124   DEV_SHAPE_FUN_MACRO_BEGIN;
3125
3126   devFunValue[0] = 0.125*(1. + gc[1] + gc[2] + 2* gc[0]) * (1.0 - gc[1])*(1.0 - gc[2]);
3127   devFunValue[1] = 0.125*(1.0 - gc[0])*(1. + gc[0] + gc[2] + 2 * gc[1])*(1.0 - gc[2]);
3128   devFunValue[2] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[0] + gc[1] + 2 *gc[2]);
3129
3130   devFunValue[3] = 0.125*(-1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3131   devFunValue[4] = 0.125*(1.0 + gc[0])*  (1 - gc[0] + gc[2] + 2*gc[1])  *(1.0 - gc[2]);
3132   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])* (1 - gc[0] + gc[1] + 2*gc[2]);
3133
3134   devFunValue[6] = 0.125*(-1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3135   devFunValue[7] = 0.125*(1.0 + gc[0])* (-1 + gc[0] - gc[2] + 2*gc[1]) *(1.0 - gc[2]);
3136   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[0] - gc[1] + 2*gc[2]);
3137
3138   devFunValue[9] = 0.125*(1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3139   devFunValue[10] = 0.125*(1.0 - gc[0])*(-1 - gc[0] - gc[2] + 2*gc[1])*(1.0 - gc[2]);
3140   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[0] - gc[1] + 2*gc[2]);
3141
3142   devFunValue[12] = 0.125*(1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3143   devFunValue[13] = 0.125*(1.0 - gc[0])*(1 + gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
3144   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1] + 2*gc[2]);
3145
3146   devFunValue[15] = 0.125*(-1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3147   devFunValue[16] = 0.125*(1.0 + gc[0])*(1 - gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
3148   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1] + 2*gc[2]);
3149
3150   devFunValue[18] = 0.125*(-1.0 + gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3151   devFunValue[19] = 0.125*(1.0 + gc[0])*(-1 + gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
3152   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1] + 2*gc[2]);
3153
3154   devFunValue[21] = 0.125*(1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3155   devFunValue[22] = 0.125*(1.0 - gc[0])*(-1 - gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
3156   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1] + 2*gc[2]);
3157
3158   devFunValue[24] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3159   devFunValue[25] = 0.25*(1.0 - gc[0]*gc[0])*(-1)*(1.0 - gc[2]);
3160   devFunValue[26] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(-1.0);
3161
3162   devFunValue[27] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[2]);
3163   devFunValue[28] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3164   devFunValue[29] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(-1.0);
3165
3166   devFunValue[30] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3167   devFunValue[31] = 0.25*(1.0 - gc[0]*gc[0])*(1.0)*(1.0 - gc[2]);
3168   devFunValue[32] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(-1.0);
3169
3170   devFunValue[33] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 - gc[2]);
3171   devFunValue[34] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3172   devFunValue[35] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(-1.0);
3173
3174   devFunValue[36] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 - gc[1]);
3175   devFunValue[37] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(-1.0);
3176   devFunValue[38] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3177
3178   devFunValue[39] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[1]);
3179   devFunValue[40] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(-1.0);
3180   devFunValue[41] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3181
3182   devFunValue[42] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[1]);
3183   devFunValue[43] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0]);
3184   devFunValue[44] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3185
3186   devFunValue[45] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 + gc[1]);
3187   devFunValue[46] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0]);
3188   devFunValue[47] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3189
3190   devFunValue[48] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3191   devFunValue[49] = 0.25*(1.0 - gc[0]*gc[0])*(-1.0)*(1.0 + gc[2]);
3192   devFunValue[50] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
3193
3194   devFunValue[51] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[2]);
3195   devFunValue[52] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3196   devFunValue[53] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
3197
3198   devFunValue[54] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3199   devFunValue[55] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[2]);
3200   devFunValue[56] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
3201
3202   devFunValue[57] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 + gc[2]);
3203   devFunValue[58] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3204   devFunValue[59] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
3205
3206   SHAPE_FUN_MACRO_END;
3207 }
3208
3209 /*!
3210  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
3211  * Case B.
3212  */
3213 void GaussInfo::hexa20bInit()
3214 {
3215   LOCAL_COORD_MACRO_BEGIN;
3216   case 0:
3217     coords[0] = HEXA20B_REF[0];
3218     coords[1] = HEXA20B_REF[1];
3219     coords[2] = HEXA20B_REF[2];
3220     break;
3221   case 1:
3222     coords[0] = HEXA20B_REF[3];
3223     coords[1] = HEXA20B_REF[4];
3224     coords[2] = HEXA20B_REF[5];
3225     break;
3226   case 2:
3227     coords[0] = HEXA20B_REF[6];
3228     coords[1] = HEXA20B_REF[7];
3229     coords[2] = HEXA20B_REF[8];
3230     break;
3231   case 3:
3232     coords[0] = HEXA20B_REF[9];
3233     coords[1] = HEXA20B_REF[10];
3234     coords[2] = HEXA20B_REF[11];
3235     break;
3236   case 4:
3237     coords[0] = HEXA20B_REF[12];
3238     coords[1] = HEXA20B_REF[13];
3239     coords[2] = HEXA20B_REF[14];
3240     break;
3241   case 5:
3242     coords[0] = HEXA20B_REF[15];
3243     coords[1] = HEXA20B_REF[16];
3244     coords[2] = HEXA20B_REF[17];
3245     break;
3246   case 6:
3247     coords[0] = HEXA20B_REF[18];
3248     coords[1] = HEXA20B_REF[19];
3249     coords[2] = HEXA20B_REF[20];
3250     break;
3251   case 7:
3252     coords[0] = HEXA20B_REF[21];
3253     coords[1] = HEXA20B_REF[22];
3254     coords[2] = HEXA20B_REF[23];
3255     break;
3256   case 8:
3257     coords[0] = HEXA20B_REF[24];
3258     coords[1] = HEXA20B_REF[25];
3259     coords[2] = HEXA20B_REF[26];
3260     break;
3261   case 9:
3262     coords[0] = HEXA20B_REF[27];
3263     coords[1] = HEXA20B_REF[28];
3264     coords[2] = HEXA20B_REF[29];
3265     break;
3266   case 10:
3267     coords[0] = HEXA20B_REF[30];
3268     coords[1] = HEXA20B_REF[31];
3269     coords[2] = HEXA20B_REF[32];
3270     break;
3271   case 11:
3272     coords[0] = HEXA20B_REF[33];
3273     coords[1] = HEXA20B_REF[34];
3274     coords[2] = HEXA20B_REF[35];
3275     break;
3276   case 12:
3277     coords[0] = HEXA20B_REF[36];
3278     coords[1] = HEXA20B_REF[37];
3279     coords[2] = HEXA20B_REF[38];
3280     break;
3281   case 13:
3282     coords[0] = HEXA20B_REF[39];
3283     coords[1] = HEXA20B_REF[40];
3284     coords[2] = HEXA20B_REF[41];
3285     break;
3286   case 14:
3287     coords[0] = HEXA20B_REF[42];
3288     coords[1] = HEXA20B_REF[43];
3289     coords[2] = HEXA20B_REF[44];
3290     break;
3291   case 15:
3292     coords[0] = HEXA20B_REF[45];
3293     coords[1] = HEXA20B_REF[46];
3294     coords[2] = HEXA20B_REF[47];
3295     break;
3296   case 16:
3297     coords[0] = HEXA20B_REF[48];
3298     coords[1] = HEXA20B_REF[49];
3299     coords[2] = HEXA20B_REF[50];
3300     break;
3301   case 17:
3302     coords[0] = HEXA20B_REF[51];
3303     coords[1] = HEXA20B_REF[52];
3304     coords[2] = HEXA20B_REF[53];
3305     break;
3306   case 18:
3307     coords[0] = HEXA20B_REF[54];
3308     coords[1] = HEXA20B_REF[55];
3309     coords[2] = HEXA20B_REF[56];
3310     break;
3311   case 19:
3312     coords[0] = HEXA20B_REF[57];
3313     coords[1] = HEXA20B_REF[58];
3314     coords[2] = HEXA20B_REF[59];
3315     break;
3316   LOCAL_COORD_MACRO_END;
3317
3318   SHAPE_FUN_MACRO_BEGIN;
3319
3320   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3321     (-2.0 - gc[0] - gc[1] - gc[2]);
3322   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3323     (-2.0 + gc[0] - gc[1] - gc[2]);
3324   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3325     (-2.0 + gc[0] + gc[1] - gc[2]);
3326   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3327     (-2.0 - gc[0] + gc[1] - gc[2]);
3328   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3329     (-2.0 - gc[0] - gc[1] + gc[2]);
3330   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3331     (-2.0 + gc[0] - gc[1] + gc[2]);
3332   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3333     (-2.0 + gc[0] + gc[1] + gc[2]);
3334   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3335     (-2.0 - gc[0] + gc[1] + gc[2]);
3336
3337   funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3338   funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3339   funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3340   funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3341   funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3342   funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3343   funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3344   funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3345   funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3346   funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3347   funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3348   funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3349   SHAPE_FUN_MACRO_END;
3350 }
3351
3352 void GaussInfo::hexa27aInit()
3353 {
3354   LOCAL_COORD_MACRO_BEGIN;
3355   case 0:
3356     coords[0] = HEXA27A_REF[0];
3357     coords[1] = HEXA27A_REF[1];
3358     coords[2] = HEXA27A_REF[2];
3359     break;
3360   case 1:
3361     coords[0] = HEXA27A_REF[3];
3362     coords[1] = HEXA27A_REF[4];
3363     coords[2] = HEXA27A_REF[5];
3364     break;
3365   case 2:
3366     coords[0] = HEXA27A_REF[6];
3367     coords[1] = HEXA27A_REF[7];
3368     coords[2] = HEXA27A_REF[8];
3369     break;
3370   case 3:
3371     coords[0] = HEXA27A_REF[9];
3372     coords[1] = HEXA27A_REF[10];
3373     coords[2] = HEXA27A_REF[11];
3374     break;
3375   case 4:
3376     coords[0] = HEXA27A_REF[12];
3377     coords[1] = HEXA27A_REF[13];
3378     coords[2] = HEXA27A_REF[14];
3379     break;
3380   case 5:
3381     coords[0] = HEXA27A_REF[15];
3382     coords[1] = HEXA27A_REF[16];
3383     coords[2] = HEXA27A_REF[17];
3384     break;
3385   case 6:
3386     coords[0] = HEXA27A_REF[18];
3387     coords[1] = HEXA27A_REF[19];
3388     coords[2] = HEXA27A_REF[20];
3389     break;
3390   case 7:
3391     coords[0] = HEXA27A_REF[21];
3392     coords[1] = HEXA27A_REF[22];
3393     coords[2] = HEXA27A_REF[23];
3394     break;
3395   case 8:
3396     coords[0] = HEXA27A_REF[24];
3397     coords[1] = HEXA27A_REF[25];
3398     coords[2] = HEXA27A_REF[26];
3399     break;
3400   case 9:
3401     coords[0] = HEXA27A_REF[27];
3402     coords[1] = HEXA27A_REF[28];
3403     coords[2] = HEXA27A_REF[29];
3404     break;
3405   case 10:
3406     coords[0] = HEXA27A_REF[30];
3407     coords[1] = HEXA27A_REF[31];
3408     coords[2] = HEXA27A_REF[32];
3409     break;
3410   case 11:
3411     coords[0] = HEXA27A_REF[33];
3412     coords[1] = HEXA27A_REF[34];
3413     coords[2] = HEXA27A_REF[35];
3414     break;
3415   case 12:
3416     coords[0] = HEXA27A_REF[36];
3417     coords[1] = HEXA27A_REF[37];
3418     coords[2] = HEXA27A_REF[38];
3419     break;
3420   case 13:
3421     coords[0] = HEXA27A_REF[39];
3422     coords[1] = HEXA27A_REF[40];
3423     coords[2] = HEXA27A_REF[41];
3424     break;
3425   case 14:
3426     coords[0] = HEXA27A_REF[42];
3427     coords[1] = HEXA27A_REF[43];
3428     coords[2] = HEXA27A_REF[44];
3429     break;
3430   case 15:
3431     coords[0] = HEXA27A_REF[45];
3432     coords[1] = HEXA27A_REF[46];
3433     coords[2] = HEXA27A_REF[47];
3434     break;
3435   case 16:
3436     coords[0] = HEXA27A_REF[48];
3437     coords[1] = HEXA27A_REF[49];
3438     coords[2] = HEXA27A_REF[50];
3439     break;
3440   case 17:
3441     coords[0] = HEXA27A_REF[51];
3442     coords[1] = HEXA27A_REF[52];
3443     coords[2] = HEXA27A_REF[53];
3444     break;
3445   case 18:
3446     coords[0] = HEXA27A_REF[54];
3447     coords[1] = HEXA27A_REF[55];
3448     coords[2] = HEXA27A_REF[56];
3449     break;
3450   case 19:
3451     coords[0] = HEXA27A_REF[57];
3452     coords[1] = HEXA27A_REF[58];
3453     coords[2] = HEXA27A_REF[59];
3454     break;
3455   case 20:
3456     coords[0] = HEXA27A_REF[60];
3457     coords[1] = HEXA27A_REF[61];
3458     coords[2] = HEXA27A_REF[62];
3459     break;
3460   case 21:
3461     coords[0] = HEXA27A_REF[63];
3462     coords[1] = HEXA27A_REF[64];
3463     coords[2] = HEXA27A_REF[65];
3464     break;
3465   case 22:
3466     coords[0] = HEXA27A_REF[66];
3467     coords[1] = HEXA27A_REF[67];
3468     coords[2] = HEXA27A_REF[68];
3469     break;
3470   case 23:
3471     coords[0] = HEXA27A_REF[69];
3472     coords[1] = HEXA27A_REF[70];
3473     coords[2] = HEXA27A_REF[71];
3474     break;
3475   case 24:
3476     coords[0] = HEXA27A_REF[72];
3477     coords[1] = HEXA27A_REF[73];
3478     coords[2] = HEXA27A_REF[74];
3479     break;
3480   case 25:
3481     coords[0] = HEXA27A_REF[75];
3482     coords[1] = HEXA27A_REF[76];
3483     coords[2] = HEXA27A_REF[77];
3484     break;
3485   case 26:
3486     coords[0] = HEXA27A_REF[78];
3487     coords[1] = HEXA27A_REF[79];
3488     coords[2] = HEXA27A_REF[80];
3489     break;
3490   LOCAL_COORD_MACRO_END;
3491
3492   SHAPE_FUN_MACRO_BEGIN;
3493
3494   funValue[0] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3495   funValue[1] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3496   funValue[2] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3497   funValue[3] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3498   funValue[4] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3499   funValue[5] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3500   funValue[6] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3501   funValue[7] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3502   funValue[8] =0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3503   funValue[9] =0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3504   funValue[10]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3505   funValue[11]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3506   funValue[12]=0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3507   funValue[13]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3508   funValue[14]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3509   funValue[15]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3510   funValue[16]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3511   funValue[17]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3512   funValue[18]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3513   funValue[19]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3514   funValue[20]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3515   funValue[21]=0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3516   funValue[22]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3517   funValue[23]=0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3518   funValue[24]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3519   funValue[25]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3520   funValue[26]=(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3521
3522   SHAPE_FUN_MACRO_END;
3523
3524   DEV_SHAPE_FUN_MACRO_BEGIN;
3525
3526   devFunValue[0] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3527   devFunValue[1] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3528   devFunValue[2] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3529
3530   devFunValue[3] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3531   devFunValue[4] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3532   devFunValue[5] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3533
3534   devFunValue[6] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3535   devFunValue[7] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3536   devFunValue[8] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3537
3538   devFunValue[9]  = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3539   devFunValue[10] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3540   devFunValue[11] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3541
3542   devFunValue[12] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3543   devFunValue[13] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3544   devFunValue[14] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3545
3546   devFunValue[15] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3547   devFunValue[16] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3548   devFunValue[17] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3549
3550   devFunValue[18] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3551   devFunValue[19] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3552   devFunValue[20] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3553
3554   devFunValue[21] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3555   devFunValue[22] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3556   devFunValue[23] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3557
3558   devFunValue[24] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3559   devFunValue[25] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
3560   devFunValue[26] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3561
3562   devFunValue[27] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3563   devFunValue[28] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3564   devFunValue[29] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3565
3566   devFunValue[30] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3567   devFunValue[31] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
3568   devFunValue[32] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3569
3570   devFunValue[33] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3571   devFunValue[34] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3572   devFunValue[35] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3573
3574   devFunValue[36] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3575   devFunValue[37] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
3576   devFunValue[38] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3577
3578   devFunValue[39] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3579   devFunValue[40] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3580   devFunValue[41] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3581
3582   devFunValue[42] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3583   devFunValue[43] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
3584   devFunValue[44] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3585
3586   devFunValue[45] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3587   devFunValue[46] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3588   devFunValue[47] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3589
3590   devFunValue[48] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3591   devFunValue[49] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3592   devFunValue[50] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3593
3594   devFunValue[51] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3595   devFunValue[52] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3596   devFunValue[53] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3597
3598   devFunValue[54] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3599   devFunValue[55] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3600   devFunValue[56] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3601
3602   devFunValue[57] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3603   devFunValue[58] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3604   devFunValue[59] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3605
3606   devFunValue[60] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3607   devFunValue[61] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]-1.);
3608   devFunValue[62] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3609
3610   devFunValue[63] = 0.5*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3611   devFunValue[64] = 0.5*gc[0]*(gc[0]-1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
3612   devFunValue[65] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
3613
3614   devFunValue[66] = 0.5*(-2*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3615   devFunValue[67] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3616   devFunValue[68] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3617
3618   devFunValue[69] = 0.5*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3619   devFunValue[70] = 0.5*gc[0]*(gc[0]+1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
3620   devFunValue[71] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
3621
3622   devFunValue[72] = 0.5*(-2*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3623   devFunValue[73] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3624   devFunValue[74] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3625
3626   devFunValue[75] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3627   devFunValue[76] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]+1.);
3628   devFunValue[77] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3629
3630   devFunValue[78] = (-2*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3631   devFunValue[79] = (1.-gc[0]*gc[0])*(-2*gc[1])*(1.-gc[2]*gc[2]);
3632   devFunValue[80] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(-2*gc[2]);
3633
3634   DEV_SHAPE_FUN_MACRO_END;
3635 }
3636
3637 ////////////////////////////////////////////////////////////////////////////////////////////////
3638 //                                GAUSS COORD CLASS                                           //
3639 ////////////////////////////////////////////////////////////////////////////////////////////////
3640 /*!
3641  * Constructor
3642  */
3643 GaussCoords::GaussCoords()
3644 {
3645 }
3646
3647 /*!
3648  * Destructor
3649  */
3650 GaussCoords::~GaussCoords()
3651 {
3652   GaussInfoVector::iterator it = _my_gauss_info.begin();
3653   for( ; it != _my_gauss_info.end(); it++ )
3654     {
3655       if((*it) != NULL)
3656         delete (*it);
3657     }
3658 }
3659
3660 /*!
3661  * Add Gauss localization info
3662  */
3663 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
3664                                 int coordDim,
3665                                 const double* theGaussCoord,
3666                                 mcIdType theNbGauss,
3667                                 const double* theReferenceCoord,
3668                                 mcIdType theNbRef)
3669 {
3670   GaussInfoVector::iterator it = _my_gauss_info.begin();
3671   for( ; it != _my_gauss_info.end(); it++ )
3672     {
3673       if( (*it)->getCellType() == theGeometry )
3674         {
3675           break;
3676         }
3677     }
3678
3679   DataVector aGaussCoord;
3680   for(int i = 0 ; i < theNbGauss*coordDim; i++ )
3681     aGaussCoord.push_back(theGaussCoord[i]);
3682
3683   DataVector aReferenceCoord;
3684   for(int i = 0 ; i < theNbRef*coordDim; i++ )
3685     aReferenceCoord.push_back(theReferenceCoord[i]);
3686
3687
3688   GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, FromIdType<int>(theNbGauss), aReferenceCoord, FromIdType<int>(theNbRef));
3689   info->initLocalInfo();
3690
3691   //If info with cell type doesn't exist add it
3692   if( it == _my_gauss_info.end() )
3693     {
3694       _my_gauss_info.push_back(info);
3695
3696       // If information exists, update it
3697     }
3698   else
3699     {
3700       std::size_t index = std::distance(_my_gauss_info.begin(),it);
3701       delete (*it);
3702       _my_gauss_info[index] = info;
3703     }
3704 }
3705
3706
3707 /*!
3708  * Calculate gauss points coordinates
3709  */
3710 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
3711                                       const double *theNodeCoords,
3712                                       const int theSpaceDim,
3713                                       const mcIdType *theIndex)
3714 {
3715   const GaussInfo *info = getInfoGivenCellType(theGeometry);
3716   int nbCoords = theSpaceDim * info->getNbGauss();
3717   double *aCoords = new double[nbCoords];
3718   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
3719   return aCoords;
3720 }
3721
3722
3723 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
3724 {
3725   const GaussInfo *info = getInfoGivenCellType(theGeometry);
3726   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
3727 }
3728
3729 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
3730 {
3731   int aConn = info->getNbRef();
3732
3733   int nbCoords = theSpaceDim * info->getNbGauss();
3734   std::fill(result,result+nbCoords,0.);
3735
3736   for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
3737     {
3738       double *coord=result+gaussId*theSpaceDim;
3739       const double *function=info->getFunctionValues(gaussId);
3740       for ( int connId = 0; connId < aConn ; connId++ )
3741         {
3742           const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
3743           for( int dimId = 0; dimId < theSpaceDim; dimId++ )
3744             coord[dimId] += nodeCoord[dimId]*function[connId];
3745         }
3746     }
3747 }
3748
3749 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
3750 {
3751   GaussInfoVector::const_iterator it = _my_gauss_info.begin();
3752   //Try to find gauss localization info
3753   for( ; it != _my_gauss_info.end() ; it++ )
3754     if( (*it)->getCellType()==cellType)
3755       return (*it);
3756   throw INTERP_KERNEL::Exception("Can't find gauss localization information !");
3757 }