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