src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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