src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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