src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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