src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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