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   DEV_SHAPE_FUN_MACRO_BEGIN;
 921
 922   devFunValue[0] = 0.0 ;
 923   devFunValue[1] = 0.5 ;
 924
 925   devFunValue[2] = -0.5;
 926   devFunValue[3] = -0.5;
 927
 928   devFunValue[4] = 0.5;
 929   devFunValue[5] = 0.0;
 930
 931   DEV_SHAPE_FUN_MACRO_END;
 932 }
 933
 934 /*!
 935  * Init Triangle Reference coordinates ans Shape function.
 936  * Case B.
 937  */
 938 void GaussInfo::tria3bInit()
 939 {
 940   LOCAL_COORD_MACRO_BEGIN;
 941   case 0:
 942     coords[0] = TRIA3B_REF[0];
 943     coords[1] = TRIA3B_REF[1];
 944     break;
 945   case 1:
 946     coords[0] = TRIA3B_REF[2];
 947     coords[1] = TRIA3B_REF[3];
 948     break;
 949   case 2:
 950     coords[0] = TRIA3B_REF[4];
 951     coords[1] = TRIA3B_REF[5];
 952     break;
 953   LOCAL_COORD_MACRO_END;
 954
 955   SHAPE_FUN_MACRO_BEGIN;
 956   funValue[0] = 1.0 - gc[0] - gc[1];
 957   funValue[1] = gc[0];
 958   funValue[2] = gc[1];
 959   SHAPE_FUN_MACRO_END;
 960
 961   DEV_SHAPE_FUN_MACRO_BEGIN;
 962
 963   devFunValue[0] = -1.0 ;
 964   devFunValue[1] = -1.0 ;
 965
 966   devFunValue[2] = 1.0 ;
 967   devFunValue[3] = 0.0 ;
 968
 969   devFunValue[4] = 0.0 ;
 970   devFunValue[5] = 1.0 ;
 971
 972   DEV_SHAPE_FUN_MACRO_END;
 973 }
 974
 975 /*!
 976  * Init Quadratic Triangle Reference coordinates ans Shape function.
 977  * Case A.
 978  */
 979 void GaussInfo::tria6aInit()
 980 {
 981   LOCAL_COORD_MACRO_BEGIN;
 982   case 0:
 983     coords[0] = TRIA6A_REF[0];
 984     coords[1] = TRIA6A_REF[1];
 985     break;
 986   case 1:
 987     coords[0] = TRIA6A_REF[2];
 988     coords[1] = TRIA6A_REF[3];
 989     break;
 990   case 2:
 991     coords[0] = TRIA6A_REF[4];
 992     coords[1] = TRIA6A_REF[5];
 993     break;
 994   case 3:
 995     coords[0] = TRIA6A_REF[6];
 996     coords[1] = TRIA6A_REF[7];
 997     break;
 998   case 4:
 999     coords[0] = TRIA6A_REF[8];
1000     coords[1] = TRIA6A_REF[9];
1001     break;
1002   case 5:
1003     coords[0] = TRIA6A_REF[10];
1004     coords[1] = TRIA6A_REF[11];
1005     break;
1006   LOCAL_COORD_MACRO_END;
1007
1008   SHAPE_FUN_MACRO_BEGIN;
1009   funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
1010   funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
1011   funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
1012   funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
1013   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
1014   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
1015   SHAPE_FUN_MACRO_END;
1016
1017   DEV_SHAPE_FUN_MACRO_BEGIN;
1018
1019   devFunValue[0] = 0.0;
1020   devFunValue[1] = 0.5*( 2*gc[1] + 1.0 );
1021
1022   devFunValue[2] = 0.5*( 2*gc[0] + 2.0*gc[1] + 1.0);
1023   devFunValue[3] = 0.5*( 2*gc[1] + 2.0*gc[0] + 1.0);
1024
1025   devFunValue[4] = gc[0] + 0.5;
1026   devFunValue[5] = 0.0;
1027
1028   devFunValue[6] = -1.0*(1.0 + gc[1]);
1029   devFunValue[7] = -1.0*(2*gc[1]+gc[0]+1.0);
1030
1031   devFunValue[8] = -1.0*(2*gc[0]+gc[1]+1.0);
1032   devFunValue[9] = -1.0*(1.0 + gc[0]);
1033
1034   devFunValue[10] = 0.0;
1035   devFunValue[11] = (2*gc[1]+2.0);
1036
1037   DEV_SHAPE_FUN_MACRO_END;
1038 }
1039
1040 /*!
1041  * Init Quadratic Triangle Reference coordinates ans Shape function.
1042  * Case B.
1043  */
1044 void GaussInfo::tria6bInit()
1045 {
1046   LOCAL_COORD_MACRO_BEGIN;
1047   case 0:
1048     coords[0] = TRIA6B_REF[0];
1049     coords[1] = TRIA6B_REF[1];
1050     break;
1051   case 1:
1052     coords[0] = TRIA6B_REF[2];
1053     coords[1] = TRIA6B_REF[3];
1054     break;
1055   case 2:
1056     coords[0] = TRIA6B_REF[4];
1057     coords[1] = TRIA6B_REF[5];
1058     break;
1059   case 3:
1060     coords[0] = TRIA6B_REF[6];
1061     coords[1] = TRIA6B_REF[7];
1062     break;
1063   case 4:
1064     coords[0] = TRIA6B_REF[8];
1065     coords[1] = TRIA6B_REF[9];
1066     break;
1067   case 5:
1068     coords[0] = TRIA6B_REF[10];
1069     coords[1] = TRIA6B_REF[11];
1070     break;
1071   LOCAL_COORD_MACRO_END;
1072
1073   SHAPE_FUN_MACRO_BEGIN;
1074   funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
1075   funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
1076   funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
1077   funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
1078   funValue[4] = 4.0*gc[0]*gc[1];
1079   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
1080   SHAPE_FUN_MACRO_END;
1081
1082   DEV_SHAPE_FUN_MACRO_BEGIN;
1083
1084   devFunValue[0] = 4*gc[0] + 4*gc[1] - 3.0 ;
1085   devFunValue[1] = 4*gc[1] + 4*gc[0] - 3.0 ;
1086
1087   devFunValue[2] = 4*gc[0] - 1.0 ;
1088   devFunValue[3] = 0.0 ;
1089
1090   devFunValue[4] = 0.0 ;
1091   devFunValue[5] = 4*gc[1] - 1.0 ;
1092
1093   devFunValue[6] = -8.0*gc[0] - 4.0 * gc[1] + 4.0 ;
1094   devFunValue[7] = -4.0*gc[0] ;
1095
1096   devFunValue[8] = 4.0*gc[1];
1097   devFunValue[9] = 4.0*gc[0];
1098
1099   devFunValue[10] = -4.0*gc[1] ;
1100   devFunValue[11] = -8.0*gc[1] - 4.0*gc[0] + 4.0 ;
1101
1102   DEV_SHAPE_FUN_MACRO_END;
1103 }
1104
1105 void GaussInfo::tria7aInit()
1106 {
1107   LOCAL_COORD_MACRO_BEGIN;
1108   case 0:
1109     coords[0] = TRIA7A_REF[0];
1110     coords[1] = TRIA7A_REF[1];
1111     break;
1112   case 1:
1113     coords[0] = TRIA7A_REF[2];
1114     coords[1] = TRIA7A_REF[3];
1115     break;
1116   case 2:
1117     coords[0] = TRIA7A_REF[4];
1118     coords[1] = TRIA7A_REF[5];
1119     break;
1120   case 3:
1121     coords[0] = TRIA7A_REF[6];
1122     coords[1] = TRIA7A_REF[7];
1123     break;
1124   case 4:
1125     coords[0] = TRIA7A_REF[8];
1126     coords[1] = TRIA7A_REF[9];
1127     break;
1128   case 5:
1129     coords[0] = TRIA7A_REF[10];
1130     coords[1] = TRIA7A_REF[11];
1131     break;
1132   case 6:
1133     coords[0] = TRIA7A_REF[12];
1134     coords[1] = TRIA7A_REF[13];
1135     break;
1136   LOCAL_COORD_MACRO_END;
1137
1138   SHAPE_FUN_MACRO_BEGIN;
1139   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]);
1140   funValue[1]=gc[0]*(-1+2*gc[0]+3*gc[1]-3*gc[1]*(gc[0]+gc[1]));
1141   funValue[2]=gc[1]*(-1.+3.*gc[0]+2.*gc[1]-3.*gc[0]*(gc[0]+gc[1]));
1142   funValue[3]=4*gc[0]*(1-gc[0]-4*gc[1]+3*gc[1]*(gc[0]+gc[1]));
1143   funValue[4]=4*gc[0]*gc[1]*(-2+3*(gc[0]+gc[1]));
1144   funValue[5]=4*gc[1]*(1-4*gc[0]-gc[1]+3*gc[0]*(gc[0]+gc[1]));
1145   funValue[6]=27*gc[0]*gc[1]*(1-gc[0]-gc[1]);
1146   SHAPE_FUN_MACRO_END;
1147 }
1148
1149 /*!
1150  * Init Quadrangle Reference coordinates ans Shape function.
1151  * Case A.
1152  */
1153 void GaussInfo::quad4aInit()
1154 {
1155   LOCAL_COORD_MACRO_BEGIN;
1156   case 0:
1157     coords[0] = QUAD4A_REF[0];
1158     coords[1] = QUAD4A_REF[1];
1159     break;
1160   case 1:
1161     coords[0] = QUAD4A_REF[2];
1162     coords[1] = QUAD4A_REF[3];
1163     break;
1164   case 2:
1165     coords[0] = QUAD4A_REF[4];
1166     coords[1] = QUAD4A_REF[5];
1167     break;
1168   case 3:
1169     coords[0] = QUAD4A_REF[6];
1170     coords[1] = QUAD4A_REF[7];
1171     break;
1172   LOCAL_COORD_MACRO_END;
1173
1174   SHAPE_FUN_MACRO_BEGIN;
1175   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
1176   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
1177   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
1178   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0]);
1179   SHAPE_FUN_MACRO_END;
1180
1181   DEV_SHAPE_FUN_MACRO_BEGIN;
1182
1183   devFunValue[0] = -0.25*(1.0 + gc[1]);
1184   devFunValue[1] = 0.25*(1.0 - gc[0]);
1185
1186   devFunValue[2] = -0.25*(1.0 - gc[1]);
1187   devFunValue[3] = -0.25*(1.0 - gc[0]);
1188
1189   devFunValue[4] = 0.25*(1.0 - gc[1]);
1190   devFunValue[5] = -0.25*(1.0 + gc[0]);
1191
1192   devFunValue[6] = 0.25*(1.0 + gc[1]);
1193   devFunValue[7] = 0.25*(1.0 + gc[0]);
1194
1195   DEV_SHAPE_FUN_MACRO_END;
1196 }
1197
1198 /*!
1199  * Init Quadrangle Reference coordinates ans Shape function.
1200  * Case B.
1201  */
1202 void GaussInfo::quad4bInit()
1203 {
1204   LOCAL_COORD_MACRO_BEGIN;
1205   case 0:
1206     coords[0] = QUAD4B_REF[0];
1207     coords[1] = QUAD4B_REF[1];
1208     break;
1209   case 1:
1210     coords[0] = QUAD4B_REF[2];
1211     coords[1] = QUAD4B_REF[3];
1212     break;
1213   case 2:
1214     coords[0] = QUAD4B_REF[4];
1215     coords[1] = QUAD4B_REF[5];
1216     break;
1217   case 3:
1218     coords[0] = QUAD4B_REF[6];
1219     coords[1] = QUAD4B_REF[7];
1220     break;
1221   LOCAL_COORD_MACRO_END;
1222
1223   SHAPE_FUN_MACRO_BEGIN;
1224   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1225   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1226   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1227   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1228   SHAPE_FUN_MACRO_END;
1229
1230   DEV_SHAPE_FUN_MACRO_BEGIN;
1231
1232   devFunValue[0] = -0.25*(1.0 - gc[1]);
1233   devFunValue[1] = -0.25*(1.0 - gc[0]);
1234
1235   devFunValue[2] = 0.25*(1.0 - gc[1]);
1236   devFunValue[3] = -0.25*(1.0 + gc[0]);
1237
1238   devFunValue[4] = 0.25*(1.0 + gc[1]);
1239   devFunValue[5] = 0.25*(1.0 + gc[0]);
1240
1241   devFunValue[6] = -0.25*(1.0 + gc[1]);
1242   devFunValue[7] = 0.25*(1.0 - gc[0]);
1243
1244   DEV_SHAPE_FUN_MACRO_END;
1245 }
1246
1247 void GaussInfo::quad4cInit()
1248 {
1249   LOCAL_COORD_MACRO_BEGIN;
1250  case  0:
1251    coords[0] = -1.0;
1252    coords[1] = -1.0;
1253    break;
1254  case  1:
1255    coords[0] = -1.0;
1256    coords[1] =  1.0;
1257    break;
1258  case  2:
1259    coords[0] =  1.0;
1260    coords[1] =  1.0;
1261    break;
1262  case  3:
1263    coords[0] =  1.0;
1264    coords[1] = -1.0;
1265    break;
1266
1267    LOCAL_COORD_MACRO_END;
1268
1269    SHAPE_FUN_MACRO_BEGIN;
1270    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1271    funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1272    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1273    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1274    SHAPE_FUN_MACRO_END;
1275
1276   DEV_SHAPE_FUN_MACRO_BEGIN;
1277
1278   devFunValue[0] = -0.25*(1.0 - gc[1]);
1279   devFunValue[1] = -0.25*(1.0 - gc[0]);
1280
1281   devFunValue[2] = -0.25*(1.0 + gc[1]);
1282   devFunValue[3] = 0.25*(1.0 - gc[0]);
1283
1284   devFunValue[4] = 0.25*(1.0 + gc[0]);
1285   devFunValue[5] = 0.25*(1.0 + gc[1]);
1286
1287   devFunValue[6] = 0.25*(1.0 - gc[1]);
1288   devFunValue[7] = -0.25*(1.0 + gc[0]);
1289
1290   DEV_SHAPE_FUN_MACRO_END;
1291 }
1292
1293 /*!
1294  * This shapefunc map is same as degenerated seg2aInit
1295  */
1296 void GaussInfo::quad4DegSeg2Init()
1297 {
1298   LOCAL_COORD_MACRO_BEGIN;
1299  case  0:
1300    coords[0] = -1.0;
1301    coords[1] =  0.0;
1302    break;
1303  case  1:
1304    coords[0] =  1.0;
1305    coords[1] =  0.0;
1306    break;
1307  case  2:
1308    coords[0] =  0.0;
1309    coords[1] =  0.0;
1310    break;
1311  case  3:
1312    coords[0] =  0.0;
1313    coords[1] =  0.0;
1314    break;
1315    LOCAL_COORD_MACRO_END;
1316
1317    SHAPE_FUN_MACRO_BEGIN;
1318    funValue[0] = 0.5*(1.0 - gc[0]);
1319    funValue[1] = 0.5*(1.0 + gc[0]);
1320    funValue[2] = 0.;
1321    funValue[3] = 0.;
1322    SHAPE_FUN_MACRO_END;
1323 }
1324
1325 /*!
1326  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1327  * Case A.
1328  */
1329 void GaussInfo::quad8aInit()
1330 {
1331   LOCAL_COORD_MACRO_BEGIN;
1332   case 0:
1333     coords[0] = QUAD8A_REF[0];
1334     coords[1] = QUAD8A_REF[1];
1335     break;
1336   case 1:
1337     coords[0] = QUAD8A_REF[2];
1338     coords[1] = QUAD8A_REF[3];
1339     break;
1340   case 2:
1341     coords[0] = QUAD8A_REF[4];
1342     coords[1] = QUAD8A_REF[5];
1343     break;
1344   case 3:
1345     coords[0] = QUAD8A_REF[6];
1346     coords[1] = QUAD8A_REF[7];
1347     break;
1348   case 4:
1349     coords[0] = QUAD8A_REF[8];
1350     coords[1] = QUAD8A_REF[9];
1351     break;
1352   case 5:
1353     coords[0] = QUAD8A_REF[10];
1354     coords[1] = QUAD8A_REF[11];
1355     break;
1356   case 6:
1357     coords[0] = QUAD8A_REF[12];
1358     coords[1] = QUAD8A_REF[13];
1359     break;
1360   case 7:
1361     coords[0] = QUAD8A_REF[14];
1362     coords[1] = QUAD8A_REF[15];
1363     break;
1364   LOCAL_COORD_MACRO_END;
1365
1366   SHAPE_FUN_MACRO_BEGIN;
1367   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
1368   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
1369   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
1370   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
1371   funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1372   funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1373   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1374   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1375   SHAPE_FUN_MACRO_END;
1376
1377   DEV_SHAPE_FUN_MACRO_BEGIN;
1378
1379   devFunValue[0] = 0.25*(1.0 + gc[1])*(2*gc[0]-gc[1]);
1380   devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1]-gc[0]);
1381
1382   devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0]+gc[1]);
1383   devFunValue[3] = 0.25*(1.0 - gc[0])*(2*gc[1]+gc[0]);
1384
1385   devFunValue[4] = 0.25*(1.0 - gc[1])*(2*gc[0]-gc[1]);
1386   devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1]-gc[0]);
1387
1388   devFunValue[6] = 0.25*(1.0 + gc[1])*(2*gc[0]+gc[1]);
1389   devFunValue[7] = 0.25*(1.0 + gc[0])*(2*gc[1]+gc[0]);
1390
1391   devFunValue[8] = -0.5*(1.0 - gc[1])*(1.0 + gc[1]);
1392   devFunValue[9] = 0.5*(1.0 - gc[0])*(-2*gc[1]);
1393
1394   devFunValue[10] = 0.5*(1.0 - gc[1])*(-2*gc[0]);
1395   devFunValue[11] = -0.5*(1.0 - gc[0])*(1.0 + gc[0]);
1396
1397   devFunValue[12] = 0.5*(1.0 - gc[1])*(1.0 + gc[1]);
1398   devFunValue[13] = 0.5*(1.0 + gc[0])*(-2*gc[1]);
1399
1400   devFunValue[14] = 0.5*(1.0 + gc[1])*(-2*gc[0]);
1401   devFunValue[15] = 0.5*(1.0 - gc[0])*(1.0 + gc[0]);
1402
1403   DEV_SHAPE_FUN_MACRO_END;
1404 }
1405
1406 /*!
1407  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1408  * Case B.
1409  */
1410 void GaussInfo::quad8bInit()
1411 {
1412   LOCAL_COORD_MACRO_BEGIN;
1413   case 0:
1414     coords[0] = QUAD8B_REF[0];
1415     coords[1] = QUAD8B_REF[1];
1416     break;
1417   case 1:
1418     coords[0] = QUAD8B_REF[2];
1419     coords[1] = QUAD8B_REF[3];
1420     break;
1421   case 2:
1422     coords[0] = QUAD8B_REF[4];
1423     coords[1] = QUAD8B_REF[5];
1424     break;
1425   case 3:
1426     coords[0] = QUAD8B_REF[6];
1427     coords[1] = QUAD8B_REF[7];
1428     break;
1429   case 4:
1430     coords[0] = QUAD8B_REF[8];
1431     coords[1] = QUAD8B_REF[9];
1432     break;
1433   case 5:
1434     coords[0] = QUAD8B_REF[10];
1435     coords[1] = QUAD8B_REF[11];
1436     break;
1437   case 6:
1438     coords[0] = QUAD8B_REF[12];
1439     coords[1] = QUAD8B_REF[13];
1440     break;
1441   case 7:
1442     coords[0] = QUAD8B_REF[14];
1443     coords[1] = QUAD8B_REF[15];
1444     break;
1445   LOCAL_COORD_MACRO_END;
1446
1447   SHAPE_FUN_MACRO_BEGIN;
1448   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
1449   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
1450   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
1451   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
1452   funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
1453   funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
1454   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
1455   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
1456   SHAPE_FUN_MACRO_END;
1457
1458   DEV_SHAPE_FUN_MACRO_BEGIN;
1459
1460   devFunValue[0] = 0.25*(1.0 - gc[1])*(2*gc[0] + gc[1] );
1461   devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1] + gc[0] );
1462
1463   devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0] - gc[1] );
1464   devFunValue[3] = 0.25*(1.0 + gc[0])*(2*gc[1] - gc[0] );
1465
1466   devFunValue[4] = 0.25*(1.0 + gc[1])*(2*gc[0] + gc[1]);
1467   devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1] + gc[0]);
1468
1469   devFunValue[6] = 0.25*( 1.0 + gc[1])*(2*gc[0] - gc[1] );
1470   devFunValue[7] = 0.25*(1.0 - gc[0])*(2*gc[1] - gc[0]);
1471
1472   devFunValue[8] = -(1.0 - gc[1])*gc[0];
1473   devFunValue[9] = -0.5*(1.0 - gc[0]*gc[0]);
1474
1475   devFunValue[10] = 0.5*(1.0 - gc[1]*gc[1]);
1476   devFunValue[11] = -(1.0 + gc[0])*gc[1];
1477
1478   devFunValue[12] = -(1.0 + gc[1])*gc[0];
1479   devFunValue[13] = 0.5*(1.0 - gc[0]*gc[0]);
1480
1481   devFunValue[14] = -0.5*(1.0 - gc[1]*gc[1]);
1482   devFunValue[15] = -(1.0 - gc[0])*gc[1];
1483
1484   DEV_SHAPE_FUN_MACRO_END;
1485 }
1486
1487 void GaussInfo::quad9aInit()
1488 {
1489   LOCAL_COORD_MACRO_BEGIN;
1490   case 0:
1491     coords[0] = QUAD9A_REF[0];
1492     coords[1] = QUAD9A_REF[1];
1493     break;
1494   case 1:
1495     coords[0] = QUAD9A_REF[2];
1496     coords[1] = QUAD9A_REF[3];
1497     break;
1498   case 2:
1499     coords[0] = QUAD9A_REF[4];
1500     coords[1] = QUAD9A_REF[5];
1501     break;
1502   case 3:
1503     coords[0] = QUAD9A_REF[6];
1504     coords[1] = QUAD9A_REF[7];
1505     break;
1506   case 4:
1507     coords[0] = QUAD9A_REF[8];
1508     coords[1] = QUAD9A_REF[9];
1509     break;
1510   case 5:
1511     coords[0] = QUAD9A_REF[10];
1512     coords[1] = QUAD9A_REF[11];
1513     break;
1514   case 6:
1515     coords[0] = QUAD9A_REF[12];
1516     coords[1] = QUAD9A_REF[13];
1517     break;
1518   case 7:
1519     coords[0] = QUAD9A_REF[14];
1520     coords[1] = QUAD9A_REF[15];
1521     break;
1522   case 8:
1523     coords[0] = QUAD9A_REF[16];
1524     coords[1] = QUAD9A_REF[17];
1525     break;
1526   LOCAL_COORD_MACRO_END;
1527
1528   SHAPE_FUN_MACRO_BEGIN;
1529   funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
1530   funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
1531   funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
1532   funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
1533   funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
1534   funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
1535   funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
1536   funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
1537   funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
1538   SHAPE_FUN_MACRO_END;
1539 }
1540
1541 /*!
1542  * Init Tetrahedron Reference coordinates ans Shape function.
1543  * Case A.
1544  */
1545 void GaussInfo::tetra4aInit()
1546 {
1547   LOCAL_COORD_MACRO_BEGIN;
1548   case 0:
1549     coords[0] = TETRA4A_REF[0];
1550     coords[1] = TETRA4A_REF[1];
1551     coords[2] = TETRA4A_REF[2];
1552     break;
1553   case 1:
1554     coords[0] = TETRA4A_REF[3];
1555     coords[1] = TETRA4A_REF[4];
1556     coords[2] = TETRA4A_REF[5];
1557     break;
1558   case 2:
1559     coords[0] = TETRA4A_REF[6];
1560     coords[1] = TETRA4A_REF[7];
1561     coords[2] = TETRA4A_REF[8];
1562     break;
1563   case 3:
1564     coords[0] = TETRA4A_REF[9];
1565     coords[1] = TETRA4A_REF[10];
1566     coords[2] = TETRA4A_REF[11];
1567     break;
1568   LOCAL_COORD_MACRO_END;
1569
1570   SHAPE_FUN_MACRO_BEGIN;
1571   funValue[0] = gc[1];
1572   funValue[1] = gc[2];
1573   funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
1574   funValue[3] = gc[0];
1575   SHAPE_FUN_MACRO_END;
1576 }
1577
1578 /*!
1579  * Init Tetrahedron Reference coordinates ans Shape function.
1580  * Case B.
1581  */
1582 void GaussInfo::tetra4bInit()
1583 {
1584   LOCAL_COORD_MACRO_BEGIN;
1585   case 0:
1586     coords[0] = TETRA4B_REF[0];
1587     coords[1] = TETRA4B_REF[1];
1588     coords[2] = TETRA4B_REF[2];
1589     break;
1590   case 1:
1591     coords[0] = TETRA4B_REF[3];
1592     coords[1] = TETRA4B_REF[4];
1593     coords[2] = TETRA4B_REF[5];
1594     break;
1595   case 2:
1596     coords[0] = TETRA4B_REF[6];
1597     coords[1] = TETRA4B_REF[7];
1598     coords[2] = TETRA4B_REF[8];
1599     break;
1600   case 3:
1601     coords[0] = TETRA4B_REF[9];
1602     coords[1] = TETRA4B_REF[10];
1603     coords[2] = TETRA4B_REF[11];
1604     break;
1605   LOCAL_COORD_MACRO_END;
1606
1607   SHAPE_FUN_MACRO_BEGIN;
1608   funValue[0] = gc[1];
1609   funValue[2] = gc[2];
1610   funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
1611   funValue[3] = gc[0];
1612   SHAPE_FUN_MACRO_END;
1613 }
1614
1615 /*!
1616  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1617  * Case A.
1618  */
1619 void GaussInfo::tetra10aInit()
1620 {
1621   LOCAL_COORD_MACRO_BEGIN;
1622   case 0:
1623     coords[0] = TETRA10A_REF[0];
1624     coords[1] = TETRA10A_REF[1];
1625     coords[2] = TETRA10A_REF[2];
1626     break;
1627   case 1:
1628     coords[0] = TETRA10A_REF[3];
1629     coords[1] = TETRA10A_REF[4];
1630     coords[2] = TETRA10A_REF[5];
1631     break;
1632   case 2:
1633     coords[0] = TETRA10A_REF[6];
1634     coords[1] = TETRA10A_REF[7];
1635     coords[2] = TETRA10A_REF[8];
1636     break;
1637   case 3:
1638     coords[0] = TETRA10A_REF[9];
1639     coords[1] = TETRA10A_REF[10];
1640     coords[2] = TETRA10A_REF[11];
1641     break;
1642   case 4:
1643     coords[0] = TETRA10A_REF[12];
1644     coords[1] = TETRA10A_REF[13];
1645     coords[2] = TETRA10A_REF[14];
1646     break;
1647   case 5:
1648     coords[0] = TETRA10A_REF[15];
1649     coords[1] = TETRA10A_REF[16];
1650     coords[2] = TETRA10A_REF[17];
1651     break;
1652   case 6:
1653     coords[0] = TETRA10A_REF[18];
1654     coords[1] = TETRA10A_REF[19];
1655     coords[2] = TETRA10A_REF[20];
1656     break;
1657   case 7:
1658     coords[0] = TETRA10A_REF[21];
1659     coords[1] = TETRA10A_REF[22];
1660     coords[2] = TETRA10A_REF[23];
1661     break;
1662   case 8:
1663     coords[0] = TETRA10A_REF[24];
1664     coords[1] = TETRA10A_REF[25];
1665     coords[2] = TETRA10A_REF[26];
1666     break;
1667   case 9:
1668     coords[0] = TETRA10A_REF[27];
1669     coords[1] = TETRA10A_REF[28];
1670     coords[2] = TETRA10A_REF[29];
1671     break;
1672   LOCAL_COORD_MACRO_END;
1673
1674   SHAPE_FUN_MACRO_BEGIN;
1675   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1676   funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
1677   funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1678   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1679   funValue[4] = 4.0*gc[1]*gc[2];
1680   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1681   funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1682   funValue[7] = 4.0*gc[0]*gc[1];
1683   funValue[8] = 4.0*gc[0]*gc[2];
1684   funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1685   SHAPE_FUN_MACRO_END;
1686 }
1687
1688 /*!
1689  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1690  * Case B.
1691  */
1692 void GaussInfo::tetra10bInit()
1693 {
1694   LOCAL_COORD_MACRO_BEGIN;
1695   case 0:
1696     coords[0] = TETRA10B_REF[0];
1697     coords[1] = TETRA10B_REF[1];
1698     coords[2] = TETRA10B_REF[2];
1699     break;
1700   case 1:
1701     coords[0] = TETRA10B_REF[3];
1702     coords[1] = TETRA10B_REF[4];
1703     coords[2] = TETRA10B_REF[5];
1704     break;
1705   case 2:
1706     coords[0] = TETRA10B_REF[6];
1707     coords[1] = TETRA10B_REF[7];
1708     coords[2] = TETRA10B_REF[8];
1709     break;
1710   case 3:
1711     coords[0] = TETRA10B_REF[9];
1712     coords[1] = TETRA10B_REF[10];
1713     coords[2] = TETRA10B_REF[11];
1714     break;
1715   case 4:
1716     coords[0] = TETRA10B_REF[12];
1717     coords[1] = TETRA10B_REF[13];
1718     coords[2] = TETRA10B_REF[14];
1719     break;
1720   case 5:
1721     coords[0] = TETRA10B_REF[15];
1722     coords[1] = TETRA10B_REF[16];
1723     coords[2] = TETRA10B_REF[17];
1724     break;
1725   case 6:
1726     coords[0] = TETRA10B_REF[18];
1727     coords[1] = TETRA10B_REF[19];
1728     coords[2] = TETRA10B_REF[20];
1729     break;
1730   case 7:
1731     coords[0] = TETRA10B_REF[21];
1732     coords[1] = TETRA10B_REF[22];
1733     coords[2] = TETRA10B_REF[23];
1734     break;
1735   case 8:
1736     coords[0] = TETRA10B_REF[24];
1737     coords[1] = TETRA10B_REF[25];
1738     coords[2] = TETRA10B_REF[26];
1739     break;
1740   case 9:
1741     coords[0] = TETRA10B_REF[27];
1742     coords[1] = TETRA10B_REF[28];
1743     coords[2] = TETRA10B_REF[29];
1744     break;
1745   LOCAL_COORD_MACRO_END;
1746   SHAPE_FUN_MACRO_BEGIN;
1747   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1748   funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1749   funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1750   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1751   funValue[6] = 4.0*gc[1]*gc[2];
1752   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1753   funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1754   funValue[7] = 4.0*gc[0]*gc[1];
1755   funValue[9] = 4.0*gc[0]*gc[2];
1756   funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1757   SHAPE_FUN_MACRO_END;
1758 }
1759
1760 /*!
1761  * Init Pyramid Reference coordinates ans Shape function.
1762  * Case A.
1763  */
1764 void GaussInfo::pyra5aInit()
1765 {
1766   LOCAL_COORD_MACRO_BEGIN;
1767   case 0:
1768     coords[0] = PYRA5A_REF[0];
1769     coords[1] = PYRA5A_REF[1];
1770     coords[2] = PYRA5A_REF[2];
1771     break;
1772   case 1:
1773     coords[0] = PYRA5A_REF[3];
1774     coords[1] = PYRA5A_REF[4];
1775     coords[2] = PYRA5A_REF[5];
1776     break;
1777   case 2:
1778     coords[0] = PYRA5A_REF[6];
1779     coords[1] = PYRA5A_REF[7];
1780     coords[2] = PYRA5A_REF[8];
1781     break;
1782   case 3:
1783     coords[0] = PYRA5A_REF[9];
1784     coords[1] = PYRA5A_REF[10];
1785     coords[2] = PYRA5A_REF[11];
1786     break;
1787   case 4:
1788     coords[0] = PYRA5A_REF[12];
1789     coords[1] = PYRA5A_REF[13];
1790     coords[2] = PYRA5A_REF[14];
1791     break;
1792   LOCAL_COORD_MACRO_END;
1793
1794   SHAPE_FUN_MACRO_BEGIN;
1795   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1796   funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1797   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1798   funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1799   funValue[4] = gc[2];
1800   SHAPE_FUN_MACRO_END;
1801
1802   DEV_SHAPE_FUN_MACRO_BEGIN;
1803
1804   devFunValue[0] = 0.5*(gc[0]+1)*(1.0 - gc[2]);
1805   devFunValue[1] = -0.5*gc[1]*(1.0 - gc[2]);
1806   devFunValue[2] = -0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0);
1807
1808   devFunValue[3] = -0.5*gc[0]*(1.0 - gc[2]);
1809   devFunValue[4] = 0.5*(gc[1]+1)*(1.0 - gc[2]);
1810   devFunValue[5] = -0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0);
1811
1812   devFunValue[6] = 0.5*(gc[0]-1)*(1.0 - gc[2]);
1813   devFunValue[7] = -0.5*gc[1]*(1.0 - gc[2]);
1814   devFunValue[8] = -0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0);
1815
1816   devFunValue[9] = -0.5*gc[0]*(1.0 - gc[2]);
1817   devFunValue[10] = 0.5*(1*gc[1]-1)*(1.0 - gc[2]);
1818   devFunValue[11] = -0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0);
1819
1820   devFunValue[12] = 0.0;
1821   devFunValue[13] = 0.0;
1822   devFunValue[14] = 1.0;
1823
1824   DEV_SHAPE_FUN_MACRO_END;
1825 }
1826 /*!
1827  * Init Pyramid Reference coordinates ans Shape function.
1828  * Case B.
1829  */
1830 void GaussInfo::pyra5bInit()
1831 {
1832   LOCAL_COORD_MACRO_BEGIN;
1833   case 0:
1834     coords[0] = PYRA5B_REF[0];
1835     coords[1] = PYRA5B_REF[1];
1836     coords[2] = PYRA5B_REF[2];
1837     break;
1838   case 1:
1839     coords[0] = PYRA5B_REF[3];
1840     coords[1] = PYRA5B_REF[4];
1841     coords[2] = PYRA5B_REF[5];
1842     break;
1843   case 2:
1844     coords[0] = PYRA5B_REF[6];
1845     coords[1] = PYRA5B_REF[7];
1846     coords[2] = PYRA5B_REF[8];
1847     break;
1848   case 3:
1849     coords[0] = PYRA5B_REF[9];
1850     coords[1] = PYRA5B_REF[10];
1851     coords[2] = PYRA5B_REF[11];
1852     break;
1853   case 4:
1854     coords[0] = PYRA5B_REF[12];
1855     coords[1] = PYRA5B_REF[13];
1856     coords[2] = PYRA5B_REF[14];
1857     break;
1858   LOCAL_COORD_MACRO_END;
1859
1860   SHAPE_FUN_MACRO_BEGIN;
1861   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1862   funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1863   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1864   funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1865   funValue[4] = gc[2];
1866   SHAPE_FUN_MACRO_END;
1867 }
1868
1869 /*!
1870  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1871  * Case A.
1872  */
1873 void GaussInfo::pyra13aInit()
1874 {
1875   LOCAL_COORD_MACRO_BEGIN;
1876   case 0:
1877     coords[0] = PYRA13A_REF[0];
1878     coords[1] = PYRA13A_REF[1];
1879     coords[2] = PYRA13A_REF[2];
1880     break;
1881   case 1:
1882     coords[0] = PYRA13A_REF[3];
1883     coords[1] = PYRA13A_REF[4];
1884     coords[2] = PYRA13A_REF[5];
1885     break;
1886   case 2:
1887     coords[0] = PYRA13A_REF[6];
1888     coords[1] = PYRA13A_REF[7];
1889     coords[2] = PYRA13A_REF[8];
1890     break;
1891   case 3:
1892     coords[0] = PYRA13A_REF[9];
1893     coords[1] = PYRA13A_REF[10];
1894     coords[2] = PYRA13A_REF[11];
1895     break;
1896   case 4:
1897     coords[0] = PYRA13A_REF[12];
1898     coords[1] = PYRA13A_REF[13];
1899     coords[2] = PYRA13A_REF[14];
1900     break;
1901   case 5:
1902     coords[0] = PYRA13A_REF[15];
1903     coords[1] = PYRA13A_REF[16];
1904     coords[2] = PYRA13A_REF[17];
1905     break;
1906   case 6:
1907     coords[0] = PYRA13A_REF[18];
1908     coords[1] = PYRA13A_REF[19];
1909     coords[2] = PYRA13A_REF[20];
1910     break;
1911   case 7:
1912     coords[0] = PYRA13A_REF[21];
1913     coords[1] = PYRA13A_REF[22];
1914     coords[2] = PYRA13A_REF[23];
1915     break;
1916   case 8:
1917     coords[0] = PYRA13A_REF[24];
1918     coords[1] = PYRA13A_REF[25];
1919     coords[2] = PYRA13A_REF[26];
1920     break;
1921   case 9:
1922     coords[0] = PYRA13A_REF[27];
1923     coords[1] = PYRA13A_REF[28];
1924     coords[2] = PYRA13A_REF[29];
1925     break;
1926   case 10:
1927     coords[0] = PYRA13A_REF[30];
1928     coords[1] = PYRA13A_REF[31];
1929     coords[2] = PYRA13A_REF[32];
1930     break;
1931   case 11:
1932     coords[0] = PYRA13A_REF[33];
1933     coords[1] = PYRA13A_REF[34];
1934     coords[2] = PYRA13A_REF[35];
1935     break;
1936   case 12:
1937     coords[0] = PYRA13A_REF[36];
1938     coords[1] = PYRA13A_REF[37];
1939     coords[2] = PYRA13A_REF[38];
1940     break;
1941   LOCAL_COORD_MACRO_END;
1942
1943   SHAPE_FUN_MACRO_BEGIN;
1944   funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1945     (gc[0] - 0.5)/(1.0 - gc[2]);
1946   funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1947     (gc[1] - 0.5)/(1.0 - gc[2]);
1948   funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1949     (-gc[0] - 0.5)/(1.0 - gc[2]);
1950   funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1951     (-gc[1] - 0.5)/(1.0 - gc[2]);
1952
1953   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1954
1955   funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1956     (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1957   funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1958     (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1959   funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1960     (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1961   funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1962     (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1963
1964   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1965     (1.0 - gc[2]);
1966   funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1967     (1.0 - gc[2]);
1968   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1969     (1.0 - gc[2]);
1970   funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1971     (1.0 - gc[2]);
1972   SHAPE_FUN_MACRO_END;
1973 }
1974
1975 /*!
1976  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1977  * Case B.
1978  */
1979 void GaussInfo::pyra13bInit()
1980 {
1981   LOCAL_COORD_MACRO_BEGIN;
1982   case 0:
1983     coords[0] = PYRA13B_REF[0];
1984     coords[1] = PYRA13B_REF[1];
1985     coords[2] = PYRA13B_REF[2];
1986     break;
1987   case 1:
1988     coords[0] = PYRA13B_REF[3];
1989     coords[1] = PYRA13B_REF[4];
1990     coords[2] = PYRA13B_REF[5];
1991     break;
1992   case 2:
1993     coords[0] = PYRA13B_REF[6];
1994     coords[1] = PYRA13B_REF[7];
1995     coords[2] = PYRA13B_REF[8];
1996     break;
1997   case 3:
1998     coords[0] = PYRA13B_REF[9];
1999     coords[1] = PYRA13B_REF[10];
2000     coords[2] = PYRA13B_REF[11];
2001     break;
2002   case 4:
2003     coords[0] = PYRA13B_REF[12];
2004     coords[1] = PYRA13B_REF[13];
2005     coords[2] = PYRA13B_REF[14];
2006     break;
2007   case 5:
2008     coords[0] = PYRA13B_REF[15];
2009     coords[1] = PYRA13B_REF[16];
2010     coords[2] = PYRA13B_REF[17];
2011     break;
2012   case 6:
2013     coords[0] = PYRA13B_REF[18];
2014     coords[1] = PYRA13B_REF[19];
2015     coords[2] = PYRA13B_REF[20];
2016     break;
2017   case 7:
2018     coords[0] = PYRA13B_REF[21];
2019     coords[1] = PYRA13B_REF[22];
2020     coords[2] = PYRA13B_REF[23];
2021     break;
2022   case 8:
2023     coords[0] = PYRA13B_REF[24];
2024     coords[1] = PYRA13B_REF[25];
2025     coords[2] = PYRA13B_REF[26];
2026     break;
2027   case 9:
2028     coords[0] = PYRA13B_REF[27];
2029     coords[1] = PYRA13B_REF[28];
2030     coords[2] = PYRA13B_REF[29];
2031     break;
2032   case 10:
2033     coords[0] = PYRA13B_REF[30];
2034     coords[1] = PYRA13B_REF[31];
2035     coords[2] = PYRA13B_REF[32];
2036     break;
2037   case 11:
2038     coords[0] = PYRA13B_REF[33];
2039     coords[1] = PYRA13B_REF[34];
2040     coords[2] = PYRA13B_REF[35];
2041     break;
2042   case 12:
2043     coords[0] = PYRA13B_REF[36];
2044     coords[1] = PYRA13B_REF[37];
2045     coords[2] = PYRA13B_REF[38];
2046     break;
2047   LOCAL_COORD_MACRO_END;
2048
2049   SHAPE_FUN_MACRO_BEGIN;
2050   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]);
2051   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]);
2052   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]);
2053   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]);
2054
2055   funValue[4] =2.*gc[2]*(gc[2]-0.5);
2056
2057   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]);
2058   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]);
2059   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]);
2060   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]);
2061
2062   funValue[9] =gc[2]*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2063   funValue[10]=gc[2]*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2064   funValue[11]=gc[2]*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2065   funValue[12]=gc[2]*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2066
2067   SHAPE_FUN_MACRO_END;
2068 }
2069
2070
2071 /*!
2072  * Init Pentahedron Reference coordinates and Shape function.
2073  * Case A.
2074  */
2075 void GaussInfo::penta6aInit()
2076 {
2077   LOCAL_COORD_MACRO_BEGIN;
2078   case 0:
2079     coords[0] = PENTA6A_REF[0];
2080     coords[1] = PENTA6A_REF[1];
2081     coords[2] = PENTA6A_REF[2];
2082     break;
2083   case 1:
2084     coords[0] = PENTA6A_REF[3];
2085     coords[1] = PENTA6A_REF[4];
2086     coords[2] = PENTA6A_REF[5];
2087     break;
2088   case 2:
2089     coords[0] = PENTA6A_REF[6];
2090     coords[1] = PENTA6A_REF[7];
2091     coords[2] = PENTA6A_REF[8];
2092     break;
2093   case 3:
2094     coords[0] = PENTA6A_REF[9];
2095     coords[1] = PENTA6A_REF[10];
2096     coords[2] = PENTA6A_REF[11];
2097     break;
2098   case 4:
2099     coords[0] = PENTA6A_REF[12];
2100     coords[1] = PENTA6A_REF[13];
2101     coords[2] = PENTA6A_REF[14];
2102     break;
2103   case 5:
2104     coords[0] = PENTA6A_REF[15];
2105     coords[1] = PENTA6A_REF[16];
2106     coords[2] = PENTA6A_REF[17];
2107     break;
2108   LOCAL_COORD_MACRO_END;
2109
2110   SHAPE_FUN_MACRO_BEGIN;
2111   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
2112   funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
2113   funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2114
2115   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
2116   funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
2117   funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2118   SHAPE_FUN_MACRO_END;
2119
2120   DEV_SHAPE_FUN_MACRO_BEGIN;
2121
2122   devFunValue[0] = 0.5*gc[1]*(-1.0);
2123   devFunValue[1] = 0.5*(1.0 - gc[0]);
2124   devFunValue[2] = 0.0;
2125
2126   devFunValue[3] = 0.5*gc[2]*(-1.0);
2127   devFunValue[4] = 0.0;
2128   devFunValue[5] = 0.5*(1.0 - gc[0]);
2129
2130   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(-1.0);
2131   devFunValue[7] = 0.5*(-1.0)*(1.0 - gc[0]);
2132   devFunValue[8] = 0.5*(-1.0)*(1.0 - gc[0]);
2133
2134   devFunValue[9] = 0.5*gc[1];
2135   devFunValue[10] = 0.5*(gc[0] + 1.0);
2136   devFunValue[11] = 0.0;
2137
2138   devFunValue[12] = 0.5*gc[2];
2139   devFunValue[13] = 0.0;
2140   devFunValue[14] = 0.5*(gc[0] + 1.0);
2141
2142   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2]);
2143   devFunValue[16] = 0.5*(-1.0)*(1.0 + gc[0]);
2144   devFunValue[17] = 0.5*(-1.0)*(1.0 + gc[0]);
2145
2146   DEV_SHAPE_FUN_MACRO_END;
2147 }
2148
2149 /*!
2150  * Init Pentahedron Reference coordinates and Shape function.
2151  * Case B.
2152  */
2153 void GaussInfo::penta6bInit()
2154 {
2155   LOCAL_COORD_MACRO_BEGIN;
2156   case 0:
2157     coords[0] = PENTA6B_REF[0];
2158     coords[1] = PENTA6B_REF[1];
2159     coords[2] = PENTA6B_REF[2];
2160     break;
2161   case 1:
2162     coords[0] = PENTA6B_REF[3];
2163     coords[1] = PENTA6B_REF[4];
2164     coords[2] = PENTA6B_REF[5];
2165     break;
2166   case 2:
2167     coords[0] = PENTA6B_REF[6];
2168     coords[1] = PENTA6B_REF[7];
2169     coords[2] = PENTA6B_REF[8];
2170     break;
2171   case 3:
2172     coords[0] = PENTA6B_REF[9];
2173     coords[1] = PENTA6B_REF[10];
2174     coords[2] = PENTA6B_REF[11];
2175     break;
2176   case 4:
2177     coords[0] = PENTA6B_REF[12];
2178     coords[1] = PENTA6B_REF[13];
2179     coords[2] = PENTA6B_REF[14];
2180     break;
2181   case 5:
2182     coords[0] = PENTA6B_REF[15];
2183     coords[1] = PENTA6B_REF[16];
2184     coords[2] = PENTA6B_REF[17];
2185     break;
2186   LOCAL_COORD_MACRO_END;
2187
2188   SHAPE_FUN_MACRO_BEGIN;
2189   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
2190   funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
2191   funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2192   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
2193   funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
2194   funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2195   SHAPE_FUN_MACRO_END;
2196 }
2197
2198 /*!
2199  * This shapefunc map is same as degenerated tria3aInit
2200  */
2201 void GaussInfo::penta6DegTria3aInit()
2202 {
2203   LOCAL_COORD_MACRO_BEGIN;
2204  case  0:
2205    coords[0] = -1.0;
2206    coords[1] =  1.0;
2207    coords[2] =  0.0;
2208    break;
2209  case  1:
2210    coords[0] = -1.0;
2211    coords[1] = -1.0;
2212    coords[2] =  0.0;
2213    break;
2214  case  2:
2215    coords[0] =  1.0;
2216    coords[1] = -1.0;
2217    coords[2] =  0.0;
2218    break;
2219  case  3:
2220    coords[0] =  0.0;
2221    coords[1] =  0.0;
2222    coords[2] =  0.0;
2223    break;
2224  case  4:
2225    coords[0] =  0.0;
2226    coords[1] =  0.0;
2227    coords[2] =  0.0;
2228    break;
2229  case  5:
2230    coords[0] =  0.0;
2231    coords[1] =  0.0;
2232    coords[2] =  0.0;
2233    break;
2234    LOCAL_COORD_MACRO_END;
2235
2236    SHAPE_FUN_MACRO_BEGIN;
2237    funValue[0] = 0.5*(1.0 + gc[1]);
2238    funValue[1] = -0.5*(gc[0] + gc[1]);
2239    funValue[2] = 0.5*(1.0 + gc[0]);
2240    funValue[3] = 0.;
2241    funValue[4] = 0.;
2242    funValue[5] = 0.;
2243    SHAPE_FUN_MACRO_END;
2244 }
2245
2246 /*!
2247  * This shapefunc map is same as degenerated tria3bInit
2248  */
2249 void GaussInfo::penta6DegTria3bInit()
2250 {
2251   LOCAL_COORD_MACRO_BEGIN;
2252  case  0:
2253    coords[0] =  0.0;
2254    coords[1] =  0.0;
2255    coords[2] =  0.0;
2256    break;
2257  case  1:
2258    coords[0] =  1.0;
2259    coords[1] =  0.0;
2260    coords[2] =  0.0;
2261    break;
2262  case  2:
2263    coords[0] =  0.0;
2264    coords[1] =  1.0;
2265    coords[2] =  0.0;
2266    break;
2267  case  3:
2268    coords[0] =  0.0;
2269    coords[1] =  0.0;
2270    coords[2] =  0.0;
2271    break;
2272  case  4:
2273    coords[0] =  0.0;
2274    coords[1] =  0.0;
2275    coords[2] =  0.0;
2276    break;
2277  case  5:
2278    coords[0] =  0.0;
2279    coords[1] =  0.0;
2280    coords[2] =  0.0;
2281    break;
2282    LOCAL_COORD_MACRO_END;
2283
2284    SHAPE_FUN_MACRO_BEGIN;
2285    funValue[0] = 1.0 - gc[0] - gc[1];
2286    funValue[1] = gc[0];
2287    funValue[2] = gc[1];
2288    funValue[3] = 0.;
2289    funValue[4] = 0.;
2290    funValue[5] = 0.;
2291    SHAPE_FUN_MACRO_END;
2292 }
2293
2294 /*!
2295  * Init Pentahedron Reference coordinates and Shape function.
2296  * Case A.
2297  */
2298 void GaussInfo::penta15aInit()
2299 {
2300   LOCAL_COORD_MACRO_BEGIN;
2301   case 0:
2302     coords[0] = PENTA15A_REF[0];
2303     coords[1] = PENTA15A_REF[1];
2304     coords[2] = PENTA15A_REF[2];
2305     break;
2306   case 1:
2307     coords[0] = PENTA15A_REF[3];
2308     coords[1] = PENTA15A_REF[4];
2309     coords[2] = PENTA15A_REF[5];
2310     break;
2311   case 2:
2312     coords[0] = PENTA15A_REF[6];
2313     coords[1] = PENTA15A_REF[7];
2314     coords[2] = PENTA15A_REF[8];
2315     break;
2316   case 3:
2317     coords[0] = PENTA15A_REF[9];
2318     coords[1] = PENTA15A_REF[10];
2319     coords[2] = PENTA15A_REF[11];
2320     break;
2321   case 4:
2322     coords[0] = PENTA15A_REF[12];
2323     coords[1] = PENTA15A_REF[13];
2324     coords[2] = PENTA15A_REF[14];
2325     break;
2326   case 5:
2327     coords[0] = PENTA15A_REF[15];
2328     coords[1] = PENTA15A_REF[16];
2329     coords[2] = PENTA15A_REF[17];
2330     break;
2331   case 6:
2332     coords[0] = PENTA15A_REF[18];
2333     coords[1] = PENTA15A_REF[19];
2334     coords[2] = PENTA15A_REF[20];
2335     break;
2336   case 7:
2337     coords[0] = PENTA15A_REF[21];
2338     coords[1] = PENTA15A_REF[22];
2339     coords[2] = PENTA15A_REF[23];
2340     break;
2341   case 8:
2342     coords[0] = PENTA15A_REF[24];
2343     coords[1] = PENTA15A_REF[25];
2344     coords[2] = PENTA15A_REF[26];
2345     break;
2346   case 9:
2347     coords[0] = PENTA15A_REF[27];
2348     coords[1] = PENTA15A_REF[28];
2349     coords[2] = PENTA15A_REF[29];
2350     break;
2351   case 10:
2352     coords[0] = PENTA15A_REF[30];
2353     coords[1] = PENTA15A_REF[31];
2354     coords[2] = PENTA15A_REF[32];
2355     break;
2356   case 11:
2357     coords[0] = PENTA15A_REF[33];
2358     coords[1] = PENTA15A_REF[34];
2359     coords[2] = PENTA15A_REF[35];
2360     break;
2361   case 12:
2362     coords[0] = PENTA15A_REF[36];
2363     coords[1] = PENTA15A_REF[37];
2364     coords[2] = PENTA15A_REF[38];
2365     break;
2366   case 13:
2367     coords[0] = PENTA15A_REF[39];
2368     coords[1] = PENTA15A_REF[40];
2369     coords[2] = PENTA15A_REF[41];
2370     break;
2371   case 14:
2372     coords[0] = PENTA15A_REF[42];
2373     coords[1] = PENTA15A_REF[43];
2374     coords[2] = PENTA15A_REF[44];
2375     break;
2376   LOCAL_COORD_MACRO_END;
2377
2378   SHAPE_FUN_MACRO_BEGIN;
2379   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2380   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2381   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2382
2383   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2384   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2385   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2386
2387   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2388   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2389   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2390
2391   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
2392   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
2393   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2394
2395   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2396   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2397   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2398   SHAPE_FUN_MACRO_END;
2399
2400   DEV_SHAPE_FUN_MACRO_BEGIN;
2401
2402   devFunValue[0] = 0.5*gc[1]*(2 * gc[0] - 2 * gc[1] + 1.0 );
2403   devFunValue[1] = (1.0 - gc[0])*(2.0*gc[1] -1 - 0.5*gc[0]);
2404   devFunValue[2] = 0.0;
2405
2406   devFunValue[3] = 0.5*gc[2]*(2 * gc[0] - 2 * gc[2] + 1.0 );
2407   devFunValue[4] = 0.0;
2408   devFunValue[5] = (1.0 - gc[0])*(2.0*gc[2] -1 - 0.5*gc[0]);
2409
2410   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(2*gc[0] -1.0 + 2*gc[1] + 2*gc[2]);
2411   devFunValue[7] = 0.5*(gc[0] - 1.0)* ( -4*gc[1] -gc[0] -4*gc[2] + 2.0);
2412   devFunValue[8] = 0.5*(gc[0] - 1.0)* ( -4*gc[2] -4*gc[1] -gc[0] + 2.0);
2413
2414   devFunValue[9] = 0.5*gc[1]*( 2*gc[0] + 2*gc[1] -1.0 );
2415   devFunValue[10] = (1.0 + gc[0])*(2.0*gc[1] - 1.0 + 0.5*gc[0]);
2416   devFunValue[11] = 0.0;
2417
2418   devFunValue[12] = 0.5*gc[2]*( 2*gc[0] + 2*gc[2] -1.0 );
2419   devFunValue[13] = 0.0;
2420   devFunValue[14] = (1.0 + gc[0])*(2.0*gc[2] - 1.0 + 0.5*gc[0]);
2421
2422   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2])*(2 * gc[0] - 2.0 * gc[1] - 2.0 * gc[2] + 1.0 );
2423   devFunValue[16] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[1] + gc[0] - 4.0 * gc[2] + 2.0) ;
2424   devFunValue[17] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[2] - 4*gc[1] + gc[0] + 2.0 );
2425
2426   devFunValue[18] = - 2.0*gc[1]*gc[2];
2427   devFunValue[19] = 2.0*gc[2]*(1.0 - gc[0]);
2428   devFunValue[20] = 2.0*gc[1]*(1.0 - gc[0]);
2429
2430   devFunValue[21] = -2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2431   devFunValue[22] = -2.0*gc[2]*(1.0 - gc[0]);
2432   devFunValue[23] = (-4*gc[2]-2*gc[1]+2.0)*(1.0 - gc[0]);
2433
2434   devFunValue[24] = -2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2435   devFunValue[25] = (-4*gc[1]-2*gc[2]+2)*(1.0 - gc[0]);
2436   devFunValue[26] = -2.0*gc[1]*(1.0 - gc[0]);
2437
2438   devFunValue[27] = gc[1]*(- 2.0 * gc[0]);
2439   devFunValue[28] = (1.0 - gc[0]*gc[0]);
2440   devFunValue[29] = 0.0;
2441
2442   devFunValue[30] = -2.0 * gc[2] *gc[0];
2443   devFunValue[31] = 0.0;
2444   devFunValue[32] = (1.0 - gc[0]*gc[0]);
2445
2446   devFunValue[33] = -2.0 * (1.0 - gc[1] - gc[2]) * gc[0];
2447   devFunValue[34] = -1.0*(1.0 - gc[0]*gc[0]);
2448   devFunValue[35] = -1.0*(1.0 - gc[0]*gc[0]);
2449
2450   devFunValue[36] = 2.0*gc[1]*gc[2];
2451   devFunValue[37] = 2.0*gc[2]*(1.0 + gc[0]);
2452   devFunValue[38] = 2.0*gc[1]*(1.0 + gc[0]);
2453
2454   devFunValue[39] = 2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2455   devFunValue[40] = -2.0*gc[2]*(1.0 + gc[0]);
2456   devFunValue[41] = (2.0 - 2.0 * gc[1] - 4.0 * gc[2])*(1.0 + gc[0]);
2457
2458   devFunValue[42] = 2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2459   devFunValue[43] = (2.0 - 4*gc[1] - 2*gc[2])*(1.0 + gc[0]);
2460   devFunValue[44] = -2.0*gc[1]*(1.0 + gc[0]);
2461
2462   DEV_SHAPE_FUN_MACRO_END;
2463 }
2464
2465 /*!
2466  * Init Qaudratic Pentahedron Reference coordinates and Shape function.
2467  * Case B.
2468  */
2469 void GaussInfo::penta15bInit()
2470 {
2471   LOCAL_COORD_MACRO_BEGIN;
2472   case 0:
2473     coords[0] = PENTA15B_REF[0];
2474     coords[1] = PENTA15B_REF[1];
2475     coords[2] = PENTA15B_REF[2];
2476     break;
2477   case 1:
2478     coords[0] = PENTA15B_REF[3];
2479     coords[1] = PENTA15B_REF[4];
2480     coords[2] = PENTA15B_REF[5];
2481     break;
2482   case 2:
2483     coords[0] = PENTA15B_REF[6];
2484     coords[1] = PENTA15B_REF[7];
2485     coords[2] = PENTA15B_REF[8];
2486     break;
2487   case 3:
2488     coords[0] = PENTA15B_REF[9];
2489     coords[1] = PENTA15B_REF[10];
2490     coords[2] = PENTA15B_REF[11];
2491     break;
2492   case 4:
2493     coords[0] = PENTA15B_REF[12];
2494     coords[1] = PENTA15B_REF[13];
2495     coords[2] = PENTA15B_REF[14];
2496     break;
2497   case 5:
2498     coords[0] = PENTA15B_REF[15];
2499     coords[1] = PENTA15B_REF[16];
2500     coords[2] = PENTA15B_REF[17];
2501     break;
2502   case 6:
2503     coords[0] = PENTA15B_REF[18];
2504     coords[1] = PENTA15B_REF[19];
2505     coords[2] = PENTA15B_REF[20];
2506     break;
2507   case 7:
2508     coords[0] = PENTA15B_REF[21];
2509     coords[1] = PENTA15B_REF[22];
2510     coords[2] = PENTA15B_REF[23];
2511     break;
2512   case 8:
2513     coords[0] = PENTA15B_REF[24];
2514     coords[1] = PENTA15B_REF[25];
2515     coords[2] = PENTA15B_REF[26];
2516     break;
2517   case 9:
2518     coords[0] = PENTA15B_REF[27];
2519     coords[1] = PENTA15B_REF[28];
2520     coords[2] = PENTA15B_REF[29];
2521     break;
2522   case 10:
2523     coords[0] = PENTA15B_REF[30];
2524     coords[1] = PENTA15B_REF[31];
2525     coords[2] = PENTA15B_REF[32];
2526     break;
2527   case 11:
2528     coords[0] = PENTA15B_REF[33];
2529     coords[1] = PENTA15B_REF[34];
2530     coords[2] = PENTA15B_REF[35];
2531     break;
2532   case 12:
2533     coords[0] = PENTA15B_REF[36];
2534     coords[1] = PENTA15B_REF[37];
2535     coords[2] = PENTA15B_REF[38];
2536     break;
2537   case 13:
2538     coords[0] = PENTA15B_REF[39];
2539     coords[1] = PENTA15B_REF[40];
2540     coords[2] = PENTA15B_REF[41];
2541     break;
2542   case 14:
2543     coords[0] = PENTA15B_REF[42];
2544     coords[1] = PENTA15B_REF[43];
2545     coords[2] = PENTA15B_REF[44];
2546     break;
2547   LOCAL_COORD_MACRO_END;
2548
2549   SHAPE_FUN_MACRO_BEGIN;
2550   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2551   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2552   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2553
2554   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2555   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2556   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2557
2558   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2559   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2560   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2561
2562   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
2563   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
2564   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2565
2566   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2567   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2568   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2569   SHAPE_FUN_MACRO_END;
2570 }
2571
2572 void GaussInfo::penta18aInit()
2573 {
2574   LOCAL_COORD_MACRO_BEGIN;
2575   case 0:
2576     coords[0] = PENTA18A_REF[0];
2577     coords[1] = PENTA18A_REF[1];
2578     coords[2] = PENTA18A_REF[2];
2579     break;
2580   case 1:
2581     coords[0] = PENTA18A_REF[3];
2582     coords[1] = PENTA18A_REF[4];
2583     coords[2] = PENTA18A_REF[5];
2584     break;
2585   case 2:
2586     coords[0] = PENTA18A_REF[6];
2587     coords[1] = PENTA18A_REF[7];
2588     coords[2] = PENTA18A_REF[8];
2589     break;
2590   case 3:
2591     coords[0] = PENTA18A_REF[9];
2592     coords[1] = PENTA18A_REF[10];
2593     coords[2] = PENTA18A_REF[11];
2594     break;
2595   case 4:
2596     coords[0] = PENTA18A_REF[12];
2597     coords[1] = PENTA18A_REF[13];
2598     coords[2] = PENTA18A_REF[14];
2599     break;
2600   case 5:
2601     coords[0] = PENTA18A_REF[15];
2602     coords[1] = PENTA18A_REF[16];
2603     coords[2] = PENTA18A_REF[17];
2604     break;
2605   case 6:
2606     coords[0] = PENTA18A_REF[18];
2607     coords[1] = PENTA18A_REF[19];
2608     coords[2] = PENTA18A_REF[20];
2609     break;
2610   case 7:
2611     coords[0] = PENTA18A_REF[21];
2612     coords[1] = PENTA18A_REF[22];
2613     coords[2] = PENTA18A_REF[23];
2614     break;
2615   case 8:
2616     coords[0] = PENTA18A_REF[24];
2617     coords[1] = PENTA18A_REF[25];
2618     coords[2] = PENTA18A_REF[26];
2619     break;
2620   case 9:
2621     coords[0] = PENTA18A_REF[27];
2622     coords[1] = PENTA18A_REF[28];
2623     coords[2] = PENTA18A_REF[29];
2624     break;
2625   case 10:
2626     coords[0] = PENTA18A_REF[30];
2627     coords[1] = PENTA18A_REF[31];
2628     coords[2] = PENTA18A_REF[32];
2629     break;
2630   case 11:
2631     coords[0] = PENTA18A_REF[33];
2632     coords[1] = PENTA18A_REF[34];
2633     coords[2] = PENTA18A_REF[35];
2634     break;
2635   case 12:
2636     coords[0] = PENTA18A_REF[36];
2637     coords[1] = PENTA18A_REF[37];
2638     coords[2] = PENTA18A_REF[38];
2639     break;
2640   case 13:
2641     coords[0] = PENTA18A_REF[39];
2642     coords[1] = PENTA18A_REF[40];
2643     coords[2] = PENTA18A_REF[41];
2644     break;
2645   case 14:
2646     coords[0] = PENTA18A_REF[42];
2647     coords[1] = PENTA18A_REF[43];
2648     coords[2] = PENTA18A_REF[44];
2649     break;
2650   case 15:
2651     coords[0] = PENTA18A_REF[45];
2652     coords[1] = PENTA18A_REF[46];
2653     coords[2] = PENTA18A_REF[47];
2654     break;
2655   case 16:
2656     coords[0] = PENTA18A_REF[48];
2657     coords[1] = PENTA18A_REF[49];
2658     coords[2] = PENTA18A_REF[50];
2659     break;
2660   case 17:
2661     coords[0] = PENTA18A_REF[51];
2662     coords[1] = PENTA18A_REF[52];
2663     coords[2] = PENTA18A_REF[53];
2664     break;
2665   LOCAL_COORD_MACRO_END;
2666
2667   SHAPE_FUN_MACRO_BEGIN;
2668   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2669   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2670   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2671
2672   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2673   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2674   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2675
2676   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2677   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2678   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2679
2680   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
2681   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
2682   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2683
2684   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2685   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2686   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2687
2688   funValue[15] = 4.0*gc[1]*gc[2]*(1.0 - gc[0] * gc[0]);
2689   funValue[16] = 4.0*gc[2]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2690   funValue[17] = 4.0*gc[1]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2691   SHAPE_FUN_MACRO_END;
2692 }
2693
2694 void GaussInfo::penta18bInit()
2695 {
2696   LOCAL_COORD_MACRO_BEGIN;
2697   case 0:
2698     coords[0] = PENTA18B_REF[0];
2699     coords[1] = PENTA18B_REF[1];
2700     coords[2] = PENTA18B_REF[2];
2701     break;
2702   case 1:
2703     coords[0] = PENTA18B_REF[3];
2704     coords[1] = PENTA18B_REF[4];
2705     coords[2] = PENTA18B_REF[5];
2706     break;
2707   case 2:
2708     coords[0] = PENTA18B_REF[6];
2709     coords[1] = PENTA18B_REF[7];
2710     coords[2] = PENTA18B_REF[8];
2711     break;
2712   case 3:
2713     coords[0] = PENTA18B_REF[9];
2714     coords[1] = PENTA18B_REF[10];
2715     coords[2] = PENTA18B_REF[11];
2716     break;
2717   case 4:
2718     coords[0] = PENTA18B_REF[12];
2719     coords[1] = PENTA18B_REF[13];
2720     coords[2] = PENTA18B_REF[14];
2721     break;
2722   case 5:
2723     coords[0] = PENTA18B_REF[15];
2724     coords[1] = PENTA18B_REF[16];
2725     coords[2] = PENTA18B_REF[17];
2726     break;
2727   case 6:
2728     coords[0] = PENTA18B_REF[18];
2729     coords[1] = PENTA18B_REF[19];
2730     coords[2] = PENTA18B_REF[20];
2731     break;
2732   case 7:
2733     coords[0] = PENTA18B_REF[21];
2734     coords[1] = PENTA18B_REF[22];
2735     coords[2] = PENTA18B_REF[23];
2736     break;
2737   case 8:
2738     coords[0] = PENTA18B_REF[24];
2739     coords[1] = PENTA18B_REF[25];
2740     coords[2] = PENTA18B_REF[26];
2741     break;
2742   case 9:
2743     coords[0] = PENTA18B_REF[27];
2744     coords[1] = PENTA18B_REF[28];
2745     coords[2] = PENTA18B_REF[29];
2746     break;
2747   case 10:
2748     coords[0] = PENTA18B_REF[30];
2749     coords[1] = PENTA18B_REF[31];
2750     coords[2] = PENTA18B_REF[32];
2751     break;
2752   case 11:
2753     coords[0] = PENTA18B_REF[33];
2754     coords[1] = PENTA18B_REF[34];
2755     coords[2] = PENTA18B_REF[35];
2756     break;
2757   case 12:
2758     coords[0] = PENTA18B_REF[36];
2759     coords[1] = PENTA18B_REF[37];
2760     coords[2] = PENTA18B_REF[38];
2761     break;
2762   case 13:
2763     coords[0] = PENTA18B_REF[39];
2764     coords[1] = PENTA18B_REF[40];
2765     coords[2] = PENTA18B_REF[41];
2766     break;
2767   case 14:
2768     coords[0] = PENTA18B_REF[42];
2769     coords[1] = PENTA18B_REF[43];
2770     coords[2] = PENTA18B_REF[44];
2771     break;
2772   case 15:
2773     coords[0] = PENTA18B_REF[45];
2774     coords[1] = PENTA18B_REF[46];
2775     coords[2] = PENTA18B_REF[47];
2776     break;
2777   case 16:
2778     coords[0] = PENTA18B_REF[48];
2779     coords[1] = PENTA18B_REF[49];
2780     coords[2] = PENTA18B_REF[50];
2781     break;
2782   case 17:
2783     coords[0] = PENTA18B_REF[51];
2784     coords[1] = PENTA18B_REF[52];
2785     coords[2] = PENTA18B_REF[53];
2786     break;
2787   LOCAL_COORD_MACRO_END;
2788
2789   SHAPE_FUN_MACRO_BEGIN;
2790   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2791   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2792   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2793
2794   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2795   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2796   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2797
2798   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2799   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2800   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2801
2802   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
2803   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
2804   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2805
2806   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2807   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2808   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2809
2810   funValue[17] = 4.0*gc[1]*gc[2]*(1.0 - gc[0] * gc[0]);
2811   funValue[16] = 4.0*gc[2]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2812   funValue[15] = 4.0*gc[1]*(gc[0] * gc[0] - 1) * ( gc[2] + gc[1] - 1);
2813   SHAPE_FUN_MACRO_END;
2814 }
2815
2816 /*!
2817  * Init Hehahedron Reference coordinates and Shape function.
2818  * Case A.
2819  */
2820 void GaussInfo::hexa8aInit()
2821 {
2822   LOCAL_COORD_MACRO_BEGIN;
2823   case 0:
2824     coords[0] = HEXA8A_REF[0];
2825     coords[1] = HEXA8A_REF[1];
2826     coords[2] = HEXA8A_REF[2];
2827     break;
2828   case 1:
2829     coords[0] = HEXA8A_REF[3];
2830     coords[1] = HEXA8A_REF[4];
2831     coords[2] = HEXA8A_REF[5];
2832     break;
2833   case 2:
2834     coords[0] = HEXA8A_REF[6];
2835     coords[1] = HEXA8A_REF[7];
2836     coords[2] = HEXA8A_REF[8];
2837     break;
2838   case 3:
2839     coords[0] = HEXA8A_REF[9];
2840     coords[1] = HEXA8A_REF[10];
2841     coords[2] = HEXA8A_REF[11];
2842     break;
2843   case 4:
2844     coords[0] = HEXA8A_REF[12];
2845     coords[1] = HEXA8A_REF[13];
2846     coords[2] = HEXA8A_REF[14];
2847     break;
2848   case 5:
2849     coords[0] = HEXA8A_REF[15];
2850     coords[1] = HEXA8A_REF[16];
2851     coords[2] = HEXA8A_REF[17];
2852     break;
2853   case 6:
2854     coords[0] = HEXA8A_REF[18];
2855     coords[1] = HEXA8A_REF[19];
2856     coords[2] = HEXA8A_REF[20];
2857     break;
2858   case 7:
2859     coords[0] = HEXA8A_REF[21];
2860     coords[1] = HEXA8A_REF[22];
2861     coords[2] = HEXA8A_REF[23];
2862     break;
2863   LOCAL_COORD_MACRO_END;
2864
2865   SHAPE_FUN_MACRO_BEGIN;
2866   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2867   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2868   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2869   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2870
2871   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2872   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2873   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2874   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2875   SHAPE_FUN_MACRO_END;
2876
2877   DEV_SHAPE_FUN_MACRO_BEGIN;
2878
2879   devFunValue[0] = 0.125 * (-1.0) * (1.0 - gc[1])*(1.0 - gc[2]);
2880   devFunValue[1] = 0.125 * (1.0 - gc[0]) * (-1.0) * (1.0 - gc[2]);
2881   devFunValue[2] = 0.125 * (1.0 - gc[0]) * (1.0 - gc[1]) *(-1.0);
2882
2883   devFunValue[3] = 0.125*(1.0 - gc[1])*(1.0 - gc[2]);
2884   devFunValue[4] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 - gc[2]);
2885   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0);
2886
2887   devFunValue[6] = 0.125*(1.0 + gc[1])*(1.0 - gc[2]);
2888   devFunValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[2]);
2889   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0);
2890
2891   devFunValue[9 ] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 - gc[2]);
2892   devFunValue[10] = 0.125*(1.0 - gc[0])*(1.0 - gc[2]);
2893   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0);
2894
2895   devFunValue[12] = 0.125*(-1.0)*(1.0 - gc[1])*(1.0 + gc[2]);
2896   devFunValue[13] = 0.125*(1.0 - gc[0])*(-1.0)*(1.0 + gc[2]);
2897   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1]);
2898
2899   devFunValue[15] = 0.125*(1.0 - gc[1])*(1.0 + gc[2]);
2900   devFunValue[16] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 + gc[2]);
2901   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1]);
2902
2903   devFunValue[18] = 0.125*(1.0 + gc[1])*(1.0 + gc[2]);
2904   devFunValue[19] = 0.125*(1.0 + gc[0])*(1.0 + gc[2]);
2905   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1]);
2906
2907   devFunValue[21] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 + gc[2]);
2908   devFunValue[22] = 0.125*(1.0 - gc[0])*(1.0 + gc[2]);
2909   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1]);
2910
2911   DEV_SHAPE_FUN_MACRO_END;
2912 }
2913
2914 /*!
2915  * Init Hehahedron Reference coordinates and Shape function.
2916  * Case B.
2917  */
2918 void GaussInfo::hexa8bInit()
2919 {
2920   LOCAL_COORD_MACRO_BEGIN;
2921   case 0:
2922     coords[0] = HEXA8B_REF[0];
2923     coords[1] = HEXA8B_REF[1];
2924     coords[2] = HEXA8B_REF[2];
2925     break;
2926   case 1:
2927     coords[0] = HEXA8B_REF[3];
2928     coords[1] = HEXA8B_REF[4];
2929     coords[2] = HEXA8B_REF[5];
2930     break;
2931   case 2:
2932     coords[0] = HEXA8B_REF[6];
2933     coords[1] = HEXA8B_REF[7];
2934     coords[2] = HEXA8B_REF[8];
2935     break;
2936   case 3:
2937     coords[0] = HEXA8B_REF[9];
2938     coords[1] = HEXA8B_REF[10];
2939     coords[2] = HEXA8B_REF[11];
2940     break;
2941   case 4:
2942     coords[0] = HEXA8B_REF[12];
2943     coords[1] = HEXA8B_REF[13];
2944     coords[2] = HEXA8B_REF[14];
2945     break;
2946   case 5:
2947     coords[0] = HEXA8B_REF[15];
2948     coords[1] = HEXA8B_REF[16];
2949     coords[2] = HEXA8B_REF[17];
2950     break;
2951   case 6:
2952     coords[0] = HEXA8B_REF[18];
2953     coords[1] = HEXA8B_REF[19];
2954     coords[2] = HEXA8B_REF[20];
2955     break;
2956   case 7:
2957     coords[0] = HEXA8B_REF[21];
2958     coords[1] = HEXA8B_REF[22];
2959     coords[2] = HEXA8B_REF[23];
2960     break;
2961   LOCAL_COORD_MACRO_END;
2962
2963   SHAPE_FUN_MACRO_BEGIN;
2964   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2965   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2966   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2967   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2968
2969   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2970   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2971   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2972   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2973   SHAPE_FUN_MACRO_END;
2974 }
2975
2976 /*!
2977  * This shapefunc map is same as degenerated quad4bInit
2978  */
2979 void GaussInfo::hexa8DegQuad4aInit()
2980 {
2981   LOCAL_COORD_MACRO_BEGIN;
2982  case  0:
2983    coords[0] = -1.0;
2984    coords[1] =  1.0;
2985    coords[2] =  0.0;
2986    break;
2987  case  1:
2988    coords[0] = -1.0;
2989    coords[1] = -1.0;
2990    coords[2] =  0.0;
2991    break;
2992  case  2:
2993    coords[0] =  1.0;
2994    coords[1] = -1.0;
2995    coords[2] =  0.0;
2996    break;
2997  case  3:
2998    coords[0] =  1.0;
2999    coords[1] =  1.0;
3000    coords[2] =  0.0;
3001    break;
3002  case  4:
3003    coords[0] =  0.0;
3004    coords[1] =  0.0;
3005    coords[2] =  0.0;
3006    break;
3007  case  5:
3008    coords[0] =  0.0;
3009    coords[1] =  0.0;
3010    coords[2] =  0.0;
3011    break;
3012  case  6:
3013    coords[0] =  0.0;
3014    coords[1] =  0.0;
3015    coords[2] =  0.0;
3016    break;
3017  case  7:
3018    coords[0] =  0.0;
3019    coords[1] =  0.0;
3020    coords[2] =  0.0;
3021    break;
3022   LOCAL_COORD_MACRO_END;
3023
3024   SHAPE_FUN_MACRO_BEGIN;
3025   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
3026   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
3027   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
3028   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3029   funValue[4] = 0.;
3030   funValue[5] = 0.;
3031   funValue[6] = 0.;
3032   funValue[7] = 0.;
3033   SHAPE_FUN_MACRO_END;
3034 }
3035
3036 /*!
3037  * This shapefunc map is same as degenerated quad4bInit
3038  */
3039 void GaussInfo::hexa8DegQuad4bInit()
3040 {
3041   LOCAL_COORD_MACRO_BEGIN;
3042  case  0:
3043    coords[0] = -1.0;
3044    coords[1] = -1.0;
3045    coords[2] =  0.0;
3046    break;
3047  case  1:
3048    coords[0] =  1.0;
3049    coords[1] = -1.0;
3050    coords[2] =  0.0;
3051    break;
3052  case  2:
3053    coords[0] =  1.0;
3054    coords[1] =  1.0;
3055    coords[2] =  0.0;
3056    break;
3057  case  3:
3058    coords[0] = -1.0;
3059    coords[1] =  1.0;
3060    coords[2] =  0.0;
3061    break;
3062  case  4:
3063    coords[0] =  0.0;
3064    coords[1] =  0.0;
3065    coords[2] =  0.0;
3066    break;
3067  case  5:
3068    coords[0] =  0.0;
3069    coords[1] =  0.0;
3070    coords[2] =  0.0;
3071    break;
3072  case  6:
3073    coords[0] =  0.0;
3074    coords[1] =  0.0;
3075    coords[2] =  0.0;
3076    break;
3077  case  7:
3078    coords[0] =  0.0;
3079    coords[1] =  0.0;
3080    coords[2] =  0.0;
3081    break;
3082    LOCAL_COORD_MACRO_END;
3083
3084    SHAPE_FUN_MACRO_BEGIN;
3085    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
3086    funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
3087    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3088    funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
3089    funValue[4] = 0.;
3090    funValue[5] = 0.;
3091    funValue[6] = 0.;
3092    funValue[7] = 0.;
3093    SHAPE_FUN_MACRO_END;
3094 }
3095
3096 /*!
3097  * This shapefunc map is same as degenerated quad4cInit
3098  */
3099 void GaussInfo::hexa8DegQuad4cInit()
3100 {
3101   LOCAL_COORD_MACRO_BEGIN;
3102  case  0:
3103    coords[0] = -1.0;
3104    coords[1] = -1.0;
3105    coords[2] =  0.0;
3106    break;
3107  case  1:
3108    coords[0] = -1.0;
3109    coords[1] =  1.0;
3110    coords[2] =  0.0;
3111    break;
3112  case  2:
3113    coords[0] =  1.0;
3114    coords[1] =  1.0;
3115    coords[2] =  0.0;
3116    break;
3117  case  3:
3118    coords[0] =  1.0;
3119    coords[1] = -1.0;
3120    coords[2] =  0.0;
3121    break;
3122  case  4:
3123    coords[0] =  0.0;
3124    coords[1] =  0.0;
3125    coords[2] =  0.0;
3126    break;
3127  case  5:
3128    coords[0] =  0.0;
3129    coords[1] =  0.0;
3130    coords[2] =  0.0;
3131    break;
3132  case  6:
3133    coords[0] =  0.0;
3134    coords[1] =  0.0;
3135    coords[2] =  0.0;
3136    break;
3137  case  7:
3138    coords[0] =  0.0;
3139    coords[1] =  0.0;
3140    coords[2] =  0.0;
3141    break;
3142
3143    LOCAL_COORD_MACRO_END;
3144
3145    SHAPE_FUN_MACRO_BEGIN;
3146    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
3147    funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
3148    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3149    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
3150    funValue[4] = 0. ;
3151    funValue[5] = 0. ;
3152    funValue[6] = 0. ;
3153    funValue[7] = 0. ;
3154    SHAPE_FUN_MACRO_END;
3155 }
3156
3157 /*!
3158  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
3159  * Case A.
3160  */
3161 void GaussInfo::hexa20aInit()
3162 {
3163   LOCAL_COORD_MACRO_BEGIN;
3164   case 0:
3165     coords[0] = HEXA20A_REF[0];
3166     coords[1] = HEXA20A_REF[1];
3167     coords[2] = HEXA20A_REF[2];
3168     break;
3169   case 1:
3170     coords[0] = HEXA20A_REF[3];
3171     coords[1] = HEXA20A_REF[4];
3172     coords[2] = HEXA20A_REF[5];
3173     break;
3174   case 2:
3175     coords[0] = HEXA20A_REF[6];
3176     coords[1] = HEXA20A_REF[7];
3177     coords[2] = HEXA20A_REF[8];
3178     break;
3179   case 3:
3180     coords[0] = HEXA20A_REF[9];
3181     coords[1] = HEXA20A_REF[10];
3182     coords[2] = HEXA20A_REF[11];
3183     break;
3184   case 4:
3185     coords[0] = HEXA20A_REF[12];
3186     coords[1] = HEXA20A_REF[13];
3187     coords[2] = HEXA20A_REF[14];
3188     break;
3189   case 5:
3190     coords[0] = HEXA20A_REF[15];
3191     coords[1] = HEXA20A_REF[16];
3192     coords[2] = HEXA20A_REF[17];
3193     break;
3194   case 6:
3195     coords[0] = HEXA20A_REF[18];
3196     coords[1] = HEXA20A_REF[19];
3197     coords[2] = HEXA20A_REF[20];
3198     break;
3199   case 7:
3200     coords[0] = HEXA20A_REF[21];
3201     coords[1] = HEXA20A_REF[22];
3202     coords[2] = HEXA20A_REF[23];
3203     break;
3204   case 8:
3205     coords[0] = HEXA20A_REF[24];
3206     coords[1] = HEXA20A_REF[25];
3207     coords[2] = HEXA20A_REF[26];
3208     break;
3209   case 9:
3210     coords[0] = HEXA20A_REF[27];
3211     coords[1] = HEXA20A_REF[28];
3212     coords[2] = HEXA20A_REF[29];
3213     break;
3214   case 10:
3215     coords[0] = HEXA20A_REF[30];
3216     coords[1] = HEXA20A_REF[31];
3217     coords[2] = HEXA20A_REF[32];
3218     break;
3219   case 11:
3220     coords[0] = HEXA20A_REF[33];
3221     coords[1] = HEXA20A_REF[34];
3222     coords[2] = HEXA20A_REF[35];
3223     break;
3224   case 12:
3225     coords[0] = HEXA20A_REF[36];
3226     coords[1] = HEXA20A_REF[37];
3227     coords[2] = HEXA20A_REF[38];
3228     break;
3229   case 13:
3230     coords[0] = HEXA20A_REF[39];
3231     coords[1] = HEXA20A_REF[40];
3232     coords[2] = HEXA20A_REF[41];
3233     break;
3234   case 14:
3235     coords[0] = HEXA20A_REF[42];
3236     coords[1] = HEXA20A_REF[43];
3237     coords[2] = HEXA20A_REF[44];
3238     break;
3239   case 15:
3240     coords[0] = HEXA20A_REF[45];
3241     coords[1] = HEXA20A_REF[46];
3242     coords[2] = HEXA20A_REF[47];
3243     break;
3244   case 16:
3245     coords[0] = HEXA20A_REF[48];
3246     coords[1] = HEXA20A_REF[49];
3247     coords[2] = HEXA20A_REF[50];
3248     break;
3249   case 17:
3250     coords[0] = HEXA20A_REF[51];
3251     coords[1] = HEXA20A_REF[52];
3252     coords[2] = HEXA20A_REF[53];
3253     break;
3254   case 18:
3255     coords[0] = HEXA20A_REF[54];
3256     coords[1] = HEXA20A_REF[55];
3257     coords[2] = HEXA20A_REF[56];
3258     break;
3259   case 19:
3260     coords[0] = HEXA20A_REF[57];
3261     coords[1] = HEXA20A_REF[58];
3262     coords[2] = HEXA20A_REF[59];
3263     break;
3264   LOCAL_COORD_MACRO_END;
3265
3266   SHAPE_FUN_MACRO_BEGIN;
3267   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3268     (-2.0 - gc[0] - gc[1] - gc[2]);
3269   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3270     (-2.0 + gc[0] - gc[1] - gc[2]);
3271   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3272     (-2.0 + gc[0] + gc[1] - gc[2]);
3273   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3274     (-2.0 - gc[0] + gc[1] - gc[2]);
3275   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3276     (-2.0 - gc[0] - gc[1] + gc[2]);
3277   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3278     (-2.0 + gc[0] - gc[1] + gc[2]);
3279   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3280     (-2.0 + gc[0] + gc[1] + gc[2]);
3281   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3282     (-2.0 - gc[0] + gc[1] + gc[2]);
3283
3284   funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3285   funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3286   funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3287   funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3288   funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3289   funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3290   funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3291   funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3292   funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3293   funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3294   funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3295   funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3296   SHAPE_FUN_MACRO_END;
3297
3298   DEV_SHAPE_FUN_MACRO_BEGIN;
3299
3300   devFunValue[0] = 0.125*(1. + gc[1] + gc[2] + 2* gc[0]) * (1.0 - gc[1])*(1.0 - gc[2]);
3301   devFunValue[1] = 0.125*(1.0 - gc[0])*(1. + gc[0] + gc[2] + 2 * gc[1])*(1.0 - gc[2]);
3302   devFunValue[2] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[0] + gc[1] + 2 *gc[2]);
3303
3304   devFunValue[3] = 0.125*(-1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3305   devFunValue[4] = 0.125*(1.0 + gc[0])*  (1 - gc[0] + gc[2] + 2*gc[1])  *(1.0 - gc[2]);
3306   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])* (1 - gc[0] + gc[1] + 2*gc[2]);
3307
3308   devFunValue[6] = 0.125*(-1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3309   devFunValue[7] = 0.125*(1.0 + gc[0])* (-1 + gc[0] - gc[2] + 2*gc[1]) *(1.0 - gc[2]);
3310   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[0] - gc[1] + 2*gc[2]);
3311
3312   devFunValue[9] = 0.125*(1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3313   devFunValue[10] = 0.125*(1.0 - gc[0])*(-1 - gc[0] - gc[2] + 2*gc[1])*(1.0 - gc[2]);
3314   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[0] - gc[1] + 2*gc[2]);
3315
3316   devFunValue[12] = 0.125*(1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3317   devFunValue[13] = 0.125*(1.0 - gc[0])*(1 + gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
3318   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1] + 2*gc[2]);
3319
3320   devFunValue[15] = 0.125*(-1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3321   devFunValue[16] = 0.125*(1.0 + gc[0])*(1 - gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
3322   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1] + 2*gc[2]);
3323
3324   devFunValue[18] = 0.125*(-1.0 + gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3325   devFunValue[19] = 0.125*(1.0 + gc[0])*(-1 + gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
3326   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1] + 2*gc[2]);
3327
3328   devFunValue[21] = 0.125*(1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3329   devFunValue[22] = 0.125*(1.0 - gc[0])*(-1 - gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
3330   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1] + 2*gc[2]);
3331
3332   devFunValue[24] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3333   devFunValue[25] = 0.25*(1.0 - gc[0]*gc[0])*(-1)*(1.0 - gc[2]);
3334   devFunValue[26] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(-1.0);
3335
3336   devFunValue[27] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[2]);
3337   devFunValue[28] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3338   devFunValue[29] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(-1.0);
3339
3340   devFunValue[30] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3341   devFunValue[31] = 0.25*(1.0 - gc[0]*gc[0])*(1.0)*(1.0 - gc[2]);
3342   devFunValue[32] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(-1.0);
3343
3344   devFunValue[33] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 - gc[2]);
3345   devFunValue[34] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3346   devFunValue[35] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(-1.0);
3347
3348   devFunValue[36] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 - gc[1]);
3349   devFunValue[37] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(-1.0);
3350   devFunValue[38] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3351
3352   devFunValue[39] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[1]);
3353   devFunValue[40] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(-1.0);
3354   devFunValue[41] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3355
3356   devFunValue[42] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[1]);
3357   devFunValue[43] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0]);
3358   devFunValue[44] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3359
3360   devFunValue[45] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 + gc[1]);
3361   devFunValue[46] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0]);
3362   devFunValue[47] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3363
3364   devFunValue[48] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3365   devFunValue[49] = 0.25*(1.0 - gc[0]*gc[0])*(-1.0)*(1.0 + gc[2]);
3366   devFunValue[50] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
3367
3368   devFunValue[51] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[2]);
3369   devFunValue[52] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3370   devFunValue[53] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
3371
3372   devFunValue[54] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3373   devFunValue[55] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[2]);
3374   devFunValue[56] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
3375
3376   devFunValue[57] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 + gc[2]);
3377   devFunValue[58] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3378   devFunValue[59] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
3379
3380   SHAPE_FUN_MACRO_END;
3381 }
3382
3383 /*!
3384  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
3385  * Case B.
3386  */
3387 void GaussInfo::hexa20bInit()
3388 {
3389   LOCAL_COORD_MACRO_BEGIN;
3390   case 0:
3391     coords[0] = HEXA20B_REF[0];
3392     coords[1] = HEXA20B_REF[1];
3393     coords[2] = HEXA20B_REF[2];
3394     break;
3395   case 1:
3396     coords[0] = HEXA20B_REF[3];
3397     coords[1] = HEXA20B_REF[4];
3398     coords[2] = HEXA20B_REF[5];
3399     break;
3400   case 2:
3401     coords[0] = HEXA20B_REF[6];
3402     coords[1] = HEXA20B_REF[7];
3403     coords[2] = HEXA20B_REF[8];
3404     break;
3405   case 3:
3406     coords[0] = HEXA20B_REF[9];
3407     coords[1] = HEXA20B_REF[10];
3408     coords[2] = HEXA20B_REF[11];
3409     break;
3410   case 4:
3411     coords[0] = HEXA20B_REF[12];
3412     coords[1] = HEXA20B_REF[13];
3413     coords[2] = HEXA20B_REF[14];
3414     break;
3415   case 5:
3416     coords[0] = HEXA20B_REF[15];
3417     coords[1] = HEXA20B_REF[16];
3418     coords[2] = HEXA20B_REF[17];
3419     break;
3420   case 6:
3421     coords[0] = HEXA20B_REF[18];
3422     coords[1] = HEXA20B_REF[19];
3423     coords[2] = HEXA20B_REF[20];
3424     break;
3425   case 7:
3426     coords[0] = HEXA20B_REF[21];
3427     coords[1] = HEXA20B_REF[22];
3428     coords[2] = HEXA20B_REF[23];
3429     break;
3430   case 8:
3431     coords[0] = HEXA20B_REF[24];
3432     coords[1] = HEXA20B_REF[25];
3433     coords[2] = HEXA20B_REF[26];
3434     break;
3435   case 9:
3436     coords[0] = HEXA20B_REF[27];
3437     coords[1] = HEXA20B_REF[28];
3438     coords[2] = HEXA20B_REF[29];
3439     break;
3440   case 10:
3441     coords[0] = HEXA20B_REF[30];
3442     coords[1] = HEXA20B_REF[31];
3443     coords[2] = HEXA20B_REF[32];
3444     break;
3445   case 11:
3446     coords[0] = HEXA20B_REF[33];
3447     coords[1] = HEXA20B_REF[34];
3448     coords[2] = HEXA20B_REF[35];
3449     break;
3450   case 12:
3451     coords[0] = HEXA20B_REF[36];
3452     coords[1] = HEXA20B_REF[37];
3453     coords[2] = HEXA20B_REF[38];
3454     break;
3455   case 13:
3456     coords[0] = HEXA20B_REF[39];
3457     coords[1] = HEXA20B_REF[40];
3458     coords[2] = HEXA20B_REF[41];
3459     break;
3460   case 14:
3461     coords[0] = HEXA20B_REF[42];
3462     coords[1] = HEXA20B_REF[43];
3463     coords[2] = HEXA20B_REF[44];
3464     break;
3465   case 15:
3466     coords[0] = HEXA20B_REF[45];
3467     coords[1] = HEXA20B_REF[46];
3468     coords[2] = HEXA20B_REF[47];
3469     break;
3470   case 16:
3471     coords[0] = HEXA20B_REF[48];
3472     coords[1] = HEXA20B_REF[49];
3473     coords[2] = HEXA20B_REF[50];
3474     break;
3475   case 17:
3476     coords[0] = HEXA20B_REF[51];
3477     coords[1] = HEXA20B_REF[52];
3478     coords[2] = HEXA20B_REF[53];
3479     break;
3480   case 18:
3481     coords[0] = HEXA20B_REF[54];
3482     coords[1] = HEXA20B_REF[55];
3483     coords[2] = HEXA20B_REF[56];
3484     break;
3485   case 19:
3486     coords[0] = HEXA20B_REF[57];
3487     coords[1] = HEXA20B_REF[58];
3488     coords[2] = HEXA20B_REF[59];
3489     break;
3490   LOCAL_COORD_MACRO_END;
3491
3492   SHAPE_FUN_MACRO_BEGIN;
3493
3494   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3495     (-2.0 - gc[0] - gc[1] - gc[2]);
3496   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
3497     (-2.0 + gc[0] - gc[1] - gc[2]);
3498   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3499     (-2.0 + gc[0] + gc[1] - gc[2]);
3500   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
3501     (-2.0 - gc[0] + gc[1] - gc[2]);
3502   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3503     (-2.0 - gc[0] - gc[1] + gc[2]);
3504   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
3505     (-2.0 + gc[0] - gc[1] + gc[2]);
3506   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3507     (-2.0 + gc[0] + gc[1] + gc[2]);
3508   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
3509     (-2.0 - gc[0] + gc[1] + gc[2]);
3510
3511   funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3512   funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
3513   funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3514   funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
3515   funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
3516   funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
3517   funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
3518   funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
3519   funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3520   funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
3521   funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3522   funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
3523   SHAPE_FUN_MACRO_END;
3524 }
3525
3526 void GaussInfo::hexa27aInit()
3527 {
3528   LOCAL_COORD_MACRO_BEGIN;
3529   case 0:
3530     coords[0] = HEXA27A_REF[0];
3531     coords[1] = HEXA27A_REF[1];
3532     coords[2] = HEXA27A_REF[2];
3533     break;
3534   case 1:
3535     coords[0] = HEXA27A_REF[3];
3536     coords[1] = HEXA27A_REF[4];
3537     coords[2] = HEXA27A_REF[5];
3538     break;
3539   case 2:
3540     coords[0] = HEXA27A_REF[6];
3541     coords[1] = HEXA27A_REF[7];
3542     coords[2] = HEXA27A_REF[8];
3543     break;
3544   case 3:
3545     coords[0] = HEXA27A_REF[9];
3546     coords[1] = HEXA27A_REF[10];
3547     coords[2] = HEXA27A_REF[11];
3548     break;
3549   case 4:
3550     coords[0] = HEXA27A_REF[12];
3551     coords[1] = HEXA27A_REF[13];
3552     coords[2] = HEXA27A_REF[14];
3553     break;
3554   case 5:
3555     coords[0] = HEXA27A_REF[15];
3556     coords[1] = HEXA27A_REF[16];
3557     coords[2] = HEXA27A_REF[17];
3558     break;
3559   case 6:
3560     coords[0] = HEXA27A_REF[18];
3561     coords[1] = HEXA27A_REF[19];
3562     coords[2] = HEXA27A_REF[20];
3563     break;
3564   case 7:
3565     coords[0] = HEXA27A_REF[21];
3566     coords[1] = HEXA27A_REF[22];
3567     coords[2] = HEXA27A_REF[23];
3568     break;
3569   case 8:
3570     coords[0] = HEXA27A_REF[24];
3571     coords[1] = HEXA27A_REF[25];
3572     coords[2] = HEXA27A_REF[26];
3573     break;
3574   case 9:
3575     coords[0] = HEXA27A_REF[27];
3576     coords[1] = HEXA27A_REF[28];
3577     coords[2] = HEXA27A_REF[29];
3578     break;
3579   case 10:
3580     coords[0] = HEXA27A_REF[30];
3581     coords[1] = HEXA27A_REF[31];
3582     coords[2] = HEXA27A_REF[32];
3583     break;
3584   case 11:
3585     coords[0] = HEXA27A_REF[33];
3586     coords[1] = HEXA27A_REF[34];
3587     coords[2] = HEXA27A_REF[35];
3588     break;
3589   case 12:
3590     coords[0] = HEXA27A_REF[36];
3591     coords[1] = HEXA27A_REF[37];
3592     coords[2] = HEXA27A_REF[38];
3593     break;
3594   case 13:
3595     coords[0] = HEXA27A_REF[39];
3596     coords[1] = HEXA27A_REF[40];
3597     coords[2] = HEXA27A_REF[41];
3598     break;
3599   case 14:
3600     coords[0] = HEXA27A_REF[42];
3601     coords[1] = HEXA27A_REF[43];
3602     coords[2] = HEXA27A_REF[44];
3603     break;
3604   case 15:
3605     coords[0] = HEXA27A_REF[45];
3606     coords[1] = HEXA27A_REF[46];
3607     coords[2] = HEXA27A_REF[47];
3608     break;
3609   case 16:
3610     coords[0] = HEXA27A_REF[48];
3611     coords[1] = HEXA27A_REF[49];
3612     coords[2] = HEXA27A_REF[50];
3613     break;
3614   case 17:
3615     coords[0] = HEXA27A_REF[51];
3616     coords[1] = HEXA27A_REF[52];
3617     coords[2] = HEXA27A_REF[53];
3618     break;
3619   case 18:
3620     coords[0] = HEXA27A_REF[54];
3621     coords[1] = HEXA27A_REF[55];
3622     coords[2] = HEXA27A_REF[56];
3623     break;
3624   case 19:
3625     coords[0] = HEXA27A_REF[57];
3626     coords[1] = HEXA27A_REF[58];
3627     coords[2] = HEXA27A_REF[59];
3628     break;
3629   case 20:
3630     coords[0] = HEXA27A_REF[60];
3631     coords[1] = HEXA27A_REF[61];
3632     coords[2] = HEXA27A_REF[62];
3633     break;
3634   case 21:
3635     coords[0] = HEXA27A_REF[63];
3636     coords[1] = HEXA27A_REF[64];
3637     coords[2] = HEXA27A_REF[65];
3638     break;
3639   case 22:
3640     coords[0] = HEXA27A_REF[66];
3641     coords[1] = HEXA27A_REF[67];
3642     coords[2] = HEXA27A_REF[68];
3643     break;
3644   case 23:
3645     coords[0] = HEXA27A_REF[69];
3646     coords[1] = HEXA27A_REF[70];
3647     coords[2] = HEXA27A_REF[71];
3648     break;
3649   case 24:
3650     coords[0] = HEXA27A_REF[72];
3651     coords[1] = HEXA27A_REF[73];
3652     coords[2] = HEXA27A_REF[74];
3653     break;
3654   case 25:
3655     coords[0] = HEXA27A_REF[75];
3656     coords[1] = HEXA27A_REF[76];
3657     coords[2] = HEXA27A_REF[77];
3658     break;
3659   case 26:
3660     coords[0] = HEXA27A_REF[78];
3661     coords[1] = HEXA27A_REF[79];
3662     coords[2] = HEXA27A_REF[80];
3663     break;
3664   LOCAL_COORD_MACRO_END;
3665
3666   SHAPE_FUN_MACRO_BEGIN;
3667
3668   funValue[0] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3669   funValue[1] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3670   funValue[2] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3671   funValue[3] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3672   funValue[4] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3673   funValue[5] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3674   funValue[6] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3675   funValue[7] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3676   funValue[8] =0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3677   funValue[9] =0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3678   funValue[10]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3679   funValue[11]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3680   funValue[12]=0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3681   funValue[13]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3682   funValue[14]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3683   funValue[15]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3684   funValue[16]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3685   funValue[17]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3686   funValue[18]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3687   funValue[19]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3688   funValue[20]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3689   funValue[21]=0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3690   funValue[22]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3691   funValue[23]=0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3692   funValue[24]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3693   funValue[25]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3694   funValue[26]=(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3695
3696   SHAPE_FUN_MACRO_END;
3697
3698   DEV_SHAPE_FUN_MACRO_BEGIN;
3699
3700   devFunValue[0] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3701   devFunValue[1] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3702   devFunValue[2] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3703
3704   devFunValue[3] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3705   devFunValue[4] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3706   devFunValue[5] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3707
3708   devFunValue[6] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3709   devFunValue[7] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3710   devFunValue[8] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3711
3712   devFunValue[9]  = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3713   devFunValue[10] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3714   devFunValue[11] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3715
3716   devFunValue[12] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3717   devFunValue[13] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3718   devFunValue[14] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3719
3720   devFunValue[15] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3721   devFunValue[16] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3722   devFunValue[17] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3723
3724   devFunValue[18] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3725   devFunValue[19] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3726   devFunValue[20] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3727
3728   devFunValue[21] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3729   devFunValue[22] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3730   devFunValue[23] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3731
3732   devFunValue[24] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3733   devFunValue[25] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
3734   devFunValue[26] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3735
3736   devFunValue[27] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
3737   devFunValue[28] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
3738   devFunValue[29] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
3739
3740   devFunValue[30] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3741   devFunValue[31] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
3742   devFunValue[32] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3743
3744   devFunValue[33] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
3745   devFunValue[34] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
3746   devFunValue[35] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
3747
3748   devFunValue[36] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3749   devFunValue[37] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
3750   devFunValue[38] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3751
3752   devFunValue[39] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
3753   devFunValue[40] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
3754   devFunValue[41] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
3755
3756   devFunValue[42] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3757   devFunValue[43] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
3758   devFunValue[44] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3759
3760   devFunValue[45] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
3761   devFunValue[46] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
3762   devFunValue[47] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
3763
3764   devFunValue[48] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3765   devFunValue[49] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3766   devFunValue[50] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3767
3768   devFunValue[51] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3769   devFunValue[52] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3770   devFunValue[53] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3771
3772   devFunValue[54] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3773   devFunValue[55] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3774   devFunValue[56] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3775
3776   devFunValue[57] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3777   devFunValue[58] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3778   devFunValue[59] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3779
3780   devFunValue[60] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
3781   devFunValue[61] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]-1.);
3782   devFunValue[62] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
3783
3784   devFunValue[63] = 0.5*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3785   devFunValue[64] = 0.5*gc[0]*(gc[0]-1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
3786   devFunValue[65] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
3787
3788   devFunValue[66] = 0.5*(-2*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
3789   devFunValue[67] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
3790   devFunValue[68] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(-2*gc[2]);
3791
3792   devFunValue[69] = 0.5*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3793   devFunValue[70] = 0.5*gc[0]*(gc[0]+1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
3794   devFunValue[71] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
3795
3796   devFunValue[72] = 0.5*(-2*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
3797   devFunValue[73] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
3798   devFunValue[74] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(-2*gc[2]);
3799
3800   devFunValue[75] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
3801   devFunValue[76] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]+1.);
3802   devFunValue[77] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
3803
3804   devFunValue[78] = (-2*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
3805   devFunValue[79] = (1.-gc[0]*gc[0])*(-2*gc[1])*(1.-gc[2]*gc[2]);
3806   devFunValue[80] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(-2*gc[2]);
3807
3808   DEV_SHAPE_FUN_MACRO_END;
3809 }
3810
3811 ////////////////////////////////////////////////////////////////////////////////////////////////
3812 //                                GAUSS COORD CLASS                                           //
3813 ////////////////////////////////////////////////////////////////////////////////////////////////
3814 /*!
3815  * Constructor
3816  */
3817 GaussCoords::GaussCoords()
3818 {
3819 }
3820
3821 /*!
3822  * Destructor
3823  */
3824 GaussCoords::~GaussCoords()
3825 {
3826   GaussInfoVector::iterator it = _my_gauss_info.begin();
3827   for( ; it != _my_gauss_info.end(); it++ )
3828     {
3829       if((*it) != NULL)
3830         delete (*it);
3831     }
3832 }
3833
3834 /*!
3835  * Add Gauss localization info
3836  */
3837 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
3838                                 int coordDim,
3839                                 const double* theGaussCoord,
3840                                 mcIdType theNbGauss,
3841                                 const double* theReferenceCoord,
3842                                 mcIdType theNbRef)
3843 {
3844   GaussInfoVector::iterator it = _my_gauss_info.begin();
3845   for( ; it != _my_gauss_info.end(); it++ )
3846     {
3847       if( (*it)->getCellType() == theGeometry )
3848         {
3849           break;
3850         }
3851     }
3852
3853   DataVector aGaussCoord;
3854   for(int i = 0 ; i < theNbGauss*coordDim; i++ )
3855     aGaussCoord.push_back(theGaussCoord[i]);
3856
3857   DataVector aReferenceCoord;
3858   for(int i = 0 ; i < theNbRef*coordDim; i++ )
3859     aReferenceCoord.push_back(theReferenceCoord[i]);
3860
3861
3862   GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, FromIdType<int>(theNbGauss), aReferenceCoord, FromIdType<int>(theNbRef));
3863   info->initLocalInfo();
3864
3865   //If info with cell type doesn't exist add it
3866   if( it == _my_gauss_info.end() )
3867     {
3868       _my_gauss_info.push_back(info);
3869
3870       // If information exists, update it
3871     }
3872   else
3873     {
3874       std::size_t index = std::distance(_my_gauss_info.begin(),it);
3875       delete (*it);
3876       _my_gauss_info[index] = info;
3877     }
3878 }
3879
3880
3881 /*!
3882  * Calculate gauss points coordinates
3883  */
3884 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
3885                                       const double *theNodeCoords,
3886                                       const int theSpaceDim,
3887                                       const mcIdType *theIndex)
3888 {
3889   const GaussInfo *info = getInfoGivenCellType(theGeometry);
3890   int nbCoords = theSpaceDim * info->getNbGauss();
3891   double *aCoords = new double[nbCoords];
3892   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
3893   return aCoords;
3894 }
3895
3896
3897 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
3898 {
3899   const GaussInfo *info = getInfoGivenCellType(theGeometry);
3900   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
3901 }
3902
3903 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
3904 {
3905   int aConn = info->getNbRef();
3906
3907   int nbCoords = theSpaceDim * info->getNbGauss();
3908   std::fill(result,result+nbCoords,0.);
3909
3910   for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
3911     {
3912       double *coord=result+gaussId*theSpaceDim;
3913       const double *function=info->getFunctionValues(gaussId);
3914       for ( int connId = 0; connId < aConn ; connId++ )
3915         {
3916           const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
3917           for( int dimId = 0; dimId < theSpaceDim; dimId++ )
3918             coord[dimId] += nodeCoord[dimId]*function[connId];
3919         }
3920     }
3921 }
3922
3923 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
3924 {
3925   GaussInfoVector::const_iterator it = _my_gauss_info.begin();
3926   //Try to find gauss localization info
3927   for( ; it != _my_gauss_info.end() ; it++ )
3928     if( (*it)->getCellType()==cellType)
3929       return (*it);
3930   throw INTERP_KERNEL::Exception("Can't find gauss localization information !");
3931 }