src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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