src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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