src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

   1 // Copyright (C) 2007-2013  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.
   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 /*!
 191  * Initialize the internal vectors
 192  */
 193 void GaussInfo::initLocalInfo()
 194 {
 195   bool aSatify = false;
 196   const CellModel& cellModel=CellModel::GetCellModel(_my_geometry);
 197   switch( _my_geometry )
 198     {
 199     case NORM_SEG2:
 200       _my_local_ref_dim = 1;
 201       _my_local_nb_ref  = 2;
 202       seg2Init();
 203       aSatify = isSatisfy();
 204       CHECK_MACRO;
 205       break;
 206
 207     case NORM_SEG3:
 208       _my_local_ref_dim = 1;
 209       _my_local_nb_ref  = 3;
 210       seg3Init();
 211       aSatify = isSatisfy();
 212       CHECK_MACRO;
 213       break;
 214
 215     case NORM_TRI3:
 216       _my_local_ref_dim = 2;
 217       _my_local_nb_ref  = 3;
 218       tria3aInit();
 219       aSatify = isSatisfy();
 220
 221       if(!aSatify)
 222         {
 223           tria3bInit();
 224           aSatify = isSatisfy();
 225           CHECK_MACRO;
 226         }
 227       break;
 228
 229     case NORM_TRI6:
 230       _my_local_ref_dim = 2;
 231       _my_local_nb_ref  = 6;
 232       tria6aInit();
 233       aSatify = isSatisfy();
 234       if(!aSatify)
 235         {
 236           tria6bInit();
 237           aSatify = isSatisfy();
 238           CHECK_MACRO;
 239         }
 240       break;
 241
 242     case NORM_QUAD4:
 243       _my_local_ref_dim = 2;
 244       _my_local_nb_ref  = 4;
 245       quad4aInit();
 246       aSatify = isSatisfy();
 247
 248       if(!aSatify)
 249         {
 250           quad4bInit();
 251           aSatify = isSatisfy();
 252           CHECK_MACRO;
 253         }
 254       break;
 255
 256     case NORM_QUAD8:
 257       _my_local_ref_dim = 2;
 258       _my_local_nb_ref  = 8;
 259       quad8aInit();
 260       aSatify = isSatisfy();
 261
 262       if(!aSatify)
 263         {
 264           quad8bInit();
 265           aSatify = isSatisfy();
 266           CHECK_MACRO;
 267         }
 268       break;
 269
 270     case NORM_QUAD9:
 271       _my_local_ref_dim = 2;
 272       _my_local_nb_ref  = 9;
 273       quad9aInit();
 274       aSatify = isSatisfy();
 275       CHECK_MACRO;
 276       break;
 277
 278     case NORM_TETRA4:
 279       _my_local_ref_dim = 3;
 280       _my_local_nb_ref  = 4;
 281       tetra4aInit();
 282       aSatify = isSatisfy();
 283
 284       if(!aSatify)
 285         {
 286           tetra4bInit();
 287           aSatify = isSatisfy();
 288           CHECK_MACRO;
 289         }
 290       break;
 291
 292     case NORM_TETRA10:
 293       _my_local_ref_dim = 3;
 294       _my_local_nb_ref  = 10;
 295       tetra10aInit();
 296       aSatify = isSatisfy();
 297
 298       if(!aSatify)
 299         {
 300           tetra10bInit();
 301           aSatify = isSatisfy();
 302           CHECK_MACRO;
 303         }
 304       break;
 305
 306     case NORM_PYRA5:
 307       _my_local_ref_dim = 3;
 308       _my_local_nb_ref  = 5;
 309       pyra5aInit();
 310       aSatify = isSatisfy();
 311
 312       if(!aSatify)
 313         {
 314           pyra5bInit();
 315           aSatify = isSatisfy();
 316           CHECK_MACRO;
 317         }
 318       break;
 319
 320     case NORM_PYRA13:
 321       _my_local_ref_dim = 3;
 322       _my_local_nb_ref  = 13;
 323       pyra13aInit();
 324       aSatify = isSatisfy();
 325
 326       if(!aSatify)
 327         {
 328           pyra13bInit();
 329           aSatify = isSatisfy();
 330           CHECK_MACRO;
 331         }
 332       break;
 333
 334     case NORM_PENTA6:
 335       _my_local_ref_dim = 3;
 336       _my_local_nb_ref  = 6;
 337       penta6aInit();
 338       aSatify = isSatisfy();
 339
 340       if(!aSatify)
 341         {
 342           penta6bInit();
 343           aSatify = isSatisfy();
 344           CHECK_MACRO;
 345         }
 346       break;
 347
 348     case NORM_PENTA15:
 349       _my_local_ref_dim = 3;
 350       _my_local_nb_ref  = 15;
 351       penta15aInit();
 352       aSatify = isSatisfy();
 353
 354       if(!aSatify)
 355         {
 356           penta15bInit();
 357           aSatify = isSatisfy();
 358           CHECK_MACRO;
 359         }
 360       break;
 361
 362     case NORM_HEXA8:
 363       _my_local_ref_dim = 3;
 364       _my_local_nb_ref  = 8;
 365       hexa8aInit();
 366       aSatify = isSatisfy();
 367
 368       if(!aSatify)
 369         {
 370           hexa8bInit();
 371           aSatify = isSatisfy();
 372           CHECK_MACRO;
 373         }
 374       break;
 375
 376     case NORM_HEXA20:
 377       _my_local_ref_dim = 3;
 378       _my_local_nb_ref  = 20;
 379       hexa20aInit();
 380       aSatify = isSatisfy();
 381
 382       if(!aSatify)
 383         {
 384           hexa20bInit();
 385           aSatify = isSatisfy();
 386           CHECK_MACRO;
 387         }
 388       break;
 389
 390     default:
 391       throw INTERP_KERNEL::Exception("Not manged cell type !");
 392       break;
 393     }
 394 }
 395
 396 /**
 397  * Return shape function value by node id
 398  */
 399 const double* GaussInfo::getFunctionValues( const int theGaussId ) const
 400 {
 401   return &_my_function_value[ _my_nb_ref*theGaussId ];
 402 }
 403
 404 /*!
 405  * Init Segment 2 Reference coordinates ans Shape function.
 406  */
 407 void GaussInfo::seg2Init()
 408 {
 409   LOCAL_COORD_MACRO_BEGIN;
 410   case  0:
 411     coords[0] = -1.0;
 412     break;
 413   case  1:
 414     coords[0] =  1.0;
 415     break;
 416    LOCAL_COORD_MACRO_END;
 417
 418    SHAPE_FUN_MACRO_BEGIN;
 419    funValue[0] = 0.5*(1.0 - gc[0]);
 420    funValue[1] = 0.5*(1.0 + gc[0]);
 421    SHAPE_FUN_MACRO_END;
 422 }
 423
 424 /*!
 425  * Init Segment 3 Reference coordinates ans Shape function.
 426  */
 427 void GaussInfo::seg3Init()
 428 {
 429   LOCAL_COORD_MACRO_BEGIN;
 430   case  0:
 431     coords[0] = -1.0;
 432     break;
 433   case  1:
 434     coords[0] =  1.0;
 435     break;
 436   case  2:
 437     coords[0] =  0.0;
 438     break;
 439    LOCAL_COORD_MACRO_END;
 440
 441    SHAPE_FUN_MACRO_BEGIN;
 442    funValue[0] = -0.5*(1.0 - gc[0])*gc[0];
 443    funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
 444    funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
 445    SHAPE_FUN_MACRO_END;
 446 }
 447
 448 /*!
 449  * Init Triangle Reference coordinates ans Shape function.
 450  * Case A.
 451  */
 452 void GaussInfo::tria3aInit()
 453 {
 454   LOCAL_COORD_MACRO_BEGIN;
 455  case  0:
 456    coords[0] = -1.0;
 457    coords[1] =  1.0;
 458    break;
 459  case  1:
 460    coords[0] = -1.0;
 461    coords[1] = -1.0;
 462    break;
 463  case  2:
 464    coords[0] =  1.0;
 465    coords[1] = -1.0;
 466    break;
 467    LOCAL_COORD_MACRO_END;
 468
 469    SHAPE_FUN_MACRO_BEGIN;
 470    funValue[0] = 0.5*(1.0 + gc[1]);
 471    funValue[1] = -0.5*(gc[0] + gc[1]);
 472    funValue[2] = 0.5*(1.0 + gc[0]);
 473    SHAPE_FUN_MACRO_END;
 474 }
 475
 476 /*!
 477  * Init Triangle Reference coordinates ans Shape function.
 478  * Case B.
 479  */
 480 void GaussInfo::tria3bInit()
 481 {
 482   LOCAL_COORD_MACRO_BEGIN;
 483  case  0:
 484    coords[0] =  0.0;
 485    coords[1] =  0.0;
 486    break;
 487  case  1:
 488    coords[0] =  1.0;
 489    coords[1] =  0.0;
 490    break;
 491  case  2:
 492    coords[0] =  0.0;
 493    coords[1] =  1.0;
 494    break;
 495    LOCAL_COORD_MACRO_END;
 496
 497    SHAPE_FUN_MACRO_BEGIN;
 498    funValue[0] = 1.0 - gc[0] - gc[1];
 499    funValue[1] = gc[0];
 500    funValue[2] = gc[1];
 501    SHAPE_FUN_MACRO_END;
 502 }
 503
 504 /*!
 505  * Init Quadratic Triangle Reference coordinates ans Shape function.
 506  * Case A.
 507  */
 508 void GaussInfo::tria6aInit()
 509 {
 510   LOCAL_COORD_MACRO_BEGIN;
 511  case  0:
 512    coords[0] = -1.0;
 513    coords[1] =  1.0;
 514    break;
 515  case  1:
 516    coords[0] = -1.0;
 517    coords[1] = -1.0;
 518    break;
 519  case  2:
 520    coords[0] =  1.0;
 521    coords[1] = -1.0;
 522    break;
 523  case  3:
 524    coords[0] = -1.0;
 525    coords[1] =  1.0;
 526    break;
 527  case  4:
 528    coords[0] =  0.0;
 529    coords[1] = -1.0;
 530    break;
 531  case  5:
 532    coords[0] =  0.0;
 533    coords[1] =  0.0;
 534    break;
 535    LOCAL_COORD_MACRO_END;
 536
 537    SHAPE_FUN_MACRO_BEGIN;
 538    funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
 539    funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
 540    funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
 541    funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
 542    funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
 543    funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
 544    SHAPE_FUN_MACRO_END;
 545 }
 546
 547 /*!
 548  * Init Quadratic Triangle Reference coordinates ans Shape function.
 549  * Case B.
 550  */
 551 void GaussInfo::tria6bInit()
 552 {
 553   LOCAL_COORD_MACRO_BEGIN;
 554  case  0:
 555    coords[0] =  0.0;
 556    coords[1] =  0.0;
 557    break;
 558
 559  case  1:
 560    coords[0] =  1.0;
 561    coords[1] =  0.0;
 562    break;
 563
 564  case  2:
 565    coords[0] =  0.0;
 566    coords[1] =  1.0;
 567    break;
 568
 569  case  3:
 570    coords[0] =  0.5;
 571    coords[1] =  0.0;
 572    break;
 573
 574  case  4:
 575    coords[0] =  0.5;
 576    coords[1] =  0.5;
 577    break;
 578
 579  case  5:
 580    coords[0] =  0.0;
 581    coords[1] =  0.5;
 582    break;
 583
 584    LOCAL_COORD_MACRO_END;
 585
 586    SHAPE_FUN_MACRO_BEGIN;
 587    funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
 588    funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
 589    funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
 590    funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
 591    funValue[4] = 4.0*gc[0]*gc[1];
 592    funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
 593    SHAPE_FUN_MACRO_END;
 594 }
 595
 596 /*!
 597  * Init Quadrangle Reference coordinates ans Shape function.
 598  * Case A.
 599  */
 600 void GaussInfo::quad4aInit()
 601 {
 602   LOCAL_COORD_MACRO_BEGIN;
 603  case  0:
 604    coords[0] = -1.0;
 605    coords[1] =  1.0;
 606    break;
 607  case  1:
 608    coords[0] = -1.0;
 609    coords[1] = -1.0;
 610    break;
 611  case  2:
 612    coords[0] =  1.0;
 613    coords[1] = -1.0;
 614    break;
 615  case  3:
 616    coords[0] =  1.0;
 617    coords[1] =  1.0;
 618    break;
 619
 620    LOCAL_COORD_MACRO_END;
 621
 622    SHAPE_FUN_MACRO_BEGIN;
 623    funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
 624    funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
 625    funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
 626    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
 627    SHAPE_FUN_MACRO_END;
 628 }
 629
 630 /*!
 631  * Init Quadrangle Reference coordinates ans Shape function.
 632  * Case B.
 633  */
 634 void GaussInfo::quad4bInit()
 635 {
 636   LOCAL_COORD_MACRO_BEGIN;
 637  case  0:
 638    coords[0] = -1.0;
 639    coords[1] = -1.0;
 640    break;
 641  case  1:
 642    coords[0] =  1.0;
 643    coords[1] = -1.0;
 644    break;
 645  case  2:
 646    coords[0] =  1.0;
 647    coords[1] =  1.0;
 648    break;
 649  case  3:
 650    coords[0] = -1.0;
 651    coords[1] =  1.0;
 652    break;
 653
 654    LOCAL_COORD_MACRO_END;
 655
 656    SHAPE_FUN_MACRO_BEGIN;
 657    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
 658    funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
 659    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
 660    funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
 661    SHAPE_FUN_MACRO_END;
 662 }
 663
 664
 665 /*!
 666  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
 667  * Case A.
 668  */
 669 void GaussInfo::quad8aInit()
 670 {
 671   LOCAL_COORD_MACRO_BEGIN;
 672  case  0:
 673    coords[0] = -1.0;
 674    coords[1] =  1.0;
 675    break;
 676  case  1:
 677    coords[0] = -1.0;
 678    coords[1] = -1.0;
 679    break;
 680  case  2:
 681    coords[0] =  1.0;
 682    coords[1] = -1.0;
 683    break;
 684  case  3:
 685    coords[0] =  1.0;
 686    coords[1] =  1.0;
 687    break;
 688  case  4:
 689    coords[0] = -1.0;
 690    coords[1] =  0.0;
 691    break;
 692  case  5:
 693    coords[0] =  0.0;
 694    coords[1] = -1.0;
 695    break;
 696  case  6:
 697    coords[0] =  1.0;
 698    coords[1] =  0.0;
 699    break;
 700  case  7:
 701    coords[0] =  0.0;
 702    coords[1] =  1.0;
 703    break;
 704    LOCAL_COORD_MACRO_END;
 705
 706    SHAPE_FUN_MACRO_BEGIN;
 707    funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
 708    funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
 709    funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
 710    funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
 711    funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
 712    funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
 713    funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
 714    funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
 715    SHAPE_FUN_MACRO_END;
 716 }
 717
 718 /*!
 719  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
 720  * Case B.
 721  */
 722 void GaussInfo::quad8bInit()
 723 {
 724   LOCAL_COORD_MACRO_BEGIN;
 725  case  0:
 726    coords[0] = -1.0;
 727    coords[1] = -1.0;
 728    break;
 729  case  1:
 730    coords[0] =  1.0;
 731    coords[1] = -1.0;
 732    break;
 733  case  2:
 734    coords[0] =  1.0;
 735    coords[1] =  1.0;
 736    break;
 737  case  3:
 738    coords[0] = -1.0;
 739    coords[1] =  1.0;
 740    break;
 741  case  4:
 742    coords[0] =  0.0;
 743    coords[1] = -1.0;
 744    break;
 745  case  5:
 746    coords[0] =  1.0;
 747    coords[1] =  0.0;
 748    break;
 749  case  6:
 750    coords[0] =  0.0;
 751    coords[1] =  1.0;
 752    break;
 753  case  7:
 754    coords[0] = -1.0;
 755    coords[1] =  0.0;
 756    break;
 757    LOCAL_COORD_MACRO_END;
 758
 759    SHAPE_FUN_MACRO_BEGIN;
 760    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
 761    funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
 762    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
 763    funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
 764    funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
 765    funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
 766    funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
 767    funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
 768    SHAPE_FUN_MACRO_END;
 769 }
 770
 771 void GaussInfo::quad9aInit()
 772 {
 773   LOCAL_COORD_MACRO_BEGIN;
 774  case  0:
 775    coords[0] = -1.0;
 776    coords[1] = -1.0;
 777    break;
 778  case  1:
 779    coords[0] =  1.0;
 780    coords[1] = -1.0;
 781    break;
 782  case  2:
 783    coords[0] =  1.0;
 784    coords[1] =  1.0;
 785    break;
 786  case  3:
 787    coords[0] = -1.0;
 788    coords[1] =  1.0;
 789    break;
 790  case  4:
 791    coords[0] =  0.0;
 792    coords[1] = -1.0;
 793    break;
 794  case  5:
 795    coords[0] =  1.0;
 796    coords[1] =  0.0;
 797    break;
 798  case  6:
 799    coords[0] =  0.0;
 800    coords[1] =  1.0;
 801    break;
 802  case  7:
 803    coords[0] = -1.0;
 804    coords[1] =  0.0;
 805    break;
 806  case  8:
 807    coords[0] =  0.0;
 808    coords[1] =  0.0;
 809    break;
 810    LOCAL_COORD_MACRO_END;
 811
 812    SHAPE_FUN_MACRO_BEGIN;
 813    funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
 814    funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
 815    funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
 816    funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
 817    funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
 818    funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
 819    funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
 820    funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
 821    funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
 822    SHAPE_FUN_MACRO_END;
 823 }
 824
 825 /*!
 826  * Init Tetrahedron Reference coordinates ans Shape function.
 827  * Case A.
 828  */
 829 void GaussInfo::tetra4aInit()
 830 {
 831   LOCAL_COORD_MACRO_BEGIN;
 832  case  0:
 833    coords[0] =  0.0;
 834    coords[1] =  1.0;
 835    coords[2] =  0.0;
 836    break;
 837  case  1:
 838    coords[0] =  0.0;
 839    coords[1] =  0.0;
 840    coords[2] =  1.0;
 841    break;
 842  case  2:
 843    coords[0] =  0.0;
 844    coords[1] =  0.0;
 845    coords[2] =  0.0;
 846    break;
 847  case  3:
 848    coords[0] =  1.0;
 849    coords[1] =  0.0;
 850    coords[2] =  0.0;
 851    break;
 852    LOCAL_COORD_MACRO_END;
 853
 854    SHAPE_FUN_MACRO_BEGIN;
 855    funValue[0] = gc[1];
 856    funValue[1] = gc[2];
 857    funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
 858    funValue[3] = gc[0];
 859    SHAPE_FUN_MACRO_END;
 860 }
 861
 862 /*!
 863  * Init Tetrahedron Reference coordinates ans Shape function.
 864  * Case B.
 865  */
 866 void GaussInfo::tetra4bInit()
 867 {
 868   LOCAL_COORD_MACRO_BEGIN;
 869  case  0:
 870    coords[0] =  0.0;
 871    coords[1] =  1.0;
 872    coords[2] =  0.0;
 873    break;
 874  case  2:
 875    coords[0] =  0.0;
 876    coords[1] =  0.0;
 877    coords[2] =  1.0;
 878    break;
 879  case  1:
 880    coords[0] =  0.0;
 881    coords[1] =  0.0;
 882    coords[2] =  0.0;
 883    break;
 884  case  3:
 885    coords[0] =  1.0;
 886    coords[1] =  0.0;
 887    coords[2] =  0.0;
 888    break;
 889    LOCAL_COORD_MACRO_END;
 890
 891    SHAPE_FUN_MACRO_BEGIN;
 892    funValue[0] = gc[1];
 893    funValue[2] = gc[2];
 894    funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
 895    funValue[3] = gc[0];
 896    SHAPE_FUN_MACRO_END;
 897
 898 }
 899
 900 /*!
 901  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
 902  * Case A.
 903  */
 904 void GaussInfo::tetra10aInit()
 905 {
 906   LOCAL_COORD_MACRO_BEGIN;
 907  case  0:
 908    coords[0] =  0.0;
 909    coords[1] =  1.0;
 910    coords[2] =  0.0;
 911    break;
 912  case  1:
 913    coords[0] =  0.0;
 914    coords[1] =  0.0;
 915    coords[2] =  1.0;
 916    break;
 917  case  2:
 918    coords[0] =  0.0;
 919    coords[1] =  0.0;
 920    coords[2] =  0.0;
 921    break;
 922  case  3:
 923    coords[0] =  1.0;
 924    coords[1] =  0.0;
 925    coords[2] =  0.0;
 926    break;
 927  case  4:
 928    coords[0] =  0.0;
 929    coords[1] =  0.5;
 930    coords[2] =  0.5;
 931    break;
 932  case  5:
 933    coords[0] =  0.0;
 934    coords[1] =  0.0;
 935    coords[2] =  0.5;
 936    break;
 937  case  6:
 938    coords[0] =  0.0;
 939    coords[1] =  0.5;
 940    coords[2] =  0.0;
 941    break;
 942  case  7:
 943    coords[0] =  0.5;
 944    coords[1] =  0.5;
 945    coords[2] =  0.0;
 946    break;
 947  case  8:
 948    coords[0] =  0.5;
 949    coords[1] =  0.0;
 950    coords[2] =  0.5;
 951    break;
 952  case  9:
 953    coords[0] =  0.5;
 954    coords[1] =  0.0;
 955    coords[2] =  0.0;
 956    break;
 957    LOCAL_COORD_MACRO_END;
 958
 959    SHAPE_FUN_MACRO_BEGIN;
 960    funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
 961    funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
 962    funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
 963    funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
 964    funValue[4] = 4.0*gc[1]*gc[2];
 965    funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
 966    funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
 967    funValue[7] = 4.0*gc[0]*gc[1];
 968    funValue[8] = 4.0*gc[0]*gc[2];
 969    funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
 970    SHAPE_FUN_MACRO_END;
 971 }
 972
 973 /*!
 974  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
 975  * Case B.
 976  */
 977 void GaussInfo::tetra10bInit()
 978 {
 979   LOCAL_COORD_MACRO_BEGIN;
 980  case  0:
 981    coords[0] =  0.0;
 982    coords[1] =  1.0;
 983    coords[2] =  0.0;
 984    break;
 985  case  2:
 986    coords[0] =  0.0;
 987    coords[1] =  0.0;
 988    coords[2] =  1.0;
 989    break;
 990  case  1:
 991    coords[0] =  0.0;
 992    coords[1] =  0.0;
 993    coords[2] =  0.0;
 994    break;
 995  case  3:
 996    coords[0] =  1.0;
 997    coords[1] =  0.0;
 998    coords[2] =  0.0;
 999    break;
1000  case  6:
1001    coords[0] =  0.0;
1002    coords[1] =  0.5;
1003    coords[2] =  0.5;
1004    break;
1005  case  5:
1006    coords[0] =  0.0;
1007    coords[1] =  0.0;
1008    coords[2] =  0.5;
1009    break;
1010  case  4:
1011    coords[0] =  0.0;
1012    coords[1] =  0.5;
1013    coords[2] =  0.0;
1014    break;
1015  case  7:
1016    coords[0] =  0.5;
1017    coords[1] =  0.5;
1018    coords[2] =  0.0;
1019    break;
1020  case  9:
1021    coords[0] =  0.5;
1022    coords[1] =  0.0;
1023    coords[2] =  0.5;
1024    break;
1025  case  8:
1026    coords[0] =  0.5;
1027    coords[1] =  0.0;
1028    coords[2] =  0.0;
1029    break;
1030    LOCAL_COORD_MACRO_END;
1031
1032    SHAPE_FUN_MACRO_BEGIN;
1033    funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1034    funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1035    funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1036    funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1037    funValue[6] = 4.0*gc[1]*gc[2];
1038    funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1039    funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1040    funValue[7] = 4.0*gc[0]*gc[1];
1041    funValue[9] = 4.0*gc[0]*gc[2];
1042    funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1043    SHAPE_FUN_MACRO_END;
1044 }
1045
1046 /*!
1047  * Init Pyramid Reference coordinates ans Shape function.
1048  * Case A.
1049  */
1050 void GaussInfo::pyra5aInit()
1051 {
1052   LOCAL_COORD_MACRO_BEGIN;
1053  case  0:
1054    coords[0] =  1.0;
1055    coords[1] =  0.0;
1056    coords[2] =  0.0;
1057    break;
1058  case  1:
1059    coords[0] =  0.0;
1060    coords[1] =  1.0;
1061    coords[2] =  0.0;
1062    break;
1063  case  2:
1064    coords[0] = -1.0;
1065    coords[1] =  0.0;
1066    coords[2] =  0.0;
1067    break;
1068  case  3:
1069    coords[0] =  0.0;
1070    coords[1] = -1.0;
1071    coords[2] =  0.0;
1072    break;
1073  case  4:
1074    coords[0] =  0.0;
1075    coords[1] =  0.0;
1076    coords[2] =  1.0;
1077    break;
1078    LOCAL_COORD_MACRO_END;
1079
1080    SHAPE_FUN_MACRO_BEGIN;
1081    funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1082    funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1083    funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1084    funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1085    funValue[4] = gc[2];
1086    SHAPE_FUN_MACRO_END;
1087 }
1088 /*!
1089  * Init Pyramid Reference coordinates ans Shape function.
1090  * Case B.
1091  */
1092 void GaussInfo::pyra5bInit()
1093 {
1094   LOCAL_COORD_MACRO_BEGIN;
1095  case  0:
1096    coords[0] =  1.0;
1097    coords[1] =  0.0;
1098    coords[2] =  0.0;
1099    break;
1100  case  3:
1101    coords[0] =  0.0;
1102    coords[1] =  1.0;
1103    coords[2] =  0.0;
1104    break;
1105  case  2:
1106    coords[0] = -1.0;
1107    coords[1] =  0.0;
1108    coords[2] =  0.0;
1109    break;
1110  case  1:
1111    coords[0] =  0.0;
1112    coords[1] = -1.0;
1113    coords[2] =  0.0;
1114    break;
1115  case  4:
1116    coords[0] =  0.0;
1117    coords[1] =  0.0;
1118    coords[2] =  1.0;
1119    break;
1120    LOCAL_COORD_MACRO_END;
1121
1122    SHAPE_FUN_MACRO_BEGIN;
1123    funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1124    funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1125    funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
1126    funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
1127    funValue[4] = gc[2];
1128    SHAPE_FUN_MACRO_END;
1129 }
1130
1131 /*!
1132  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1133  * Case A.
1134  */
1135 void GaussInfo::pyra13aInit()
1136 {
1137   LOCAL_COORD_MACRO_BEGIN;
1138  case  0:
1139    coords[0] =  1.0;
1140    coords[1] =  0.0;
1141    coords[2] =  0.0;
1142    break;
1143  case  1:
1144    coords[0] =  0.0;
1145    coords[1] =  1.0;
1146    coords[2] =  0.0;
1147    break;
1148  case  2:
1149    coords[0] = -1.0;
1150    coords[1] =  0.0;
1151    coords[2] =  0.0;
1152    break;
1153  case  3:
1154    coords[0] =  0.0;
1155    coords[1] = -1.0;
1156    coords[2] =  0.0;
1157    break;
1158  case  4:
1159    coords[0] =  0.0;
1160    coords[1] =  0.0;
1161    coords[2] =  1.0;
1162    break;
1163
1164  case  5:
1165    coords[0] =  0.5;
1166    coords[1] =  0.5;
1167    coords[2] =  0.0;
1168    break;
1169  case  6:
1170    coords[0] = -0.5;
1171    coords[1] =  0.5;
1172    coords[2] =  0.0;
1173    break;
1174  case  7:
1175    coords[0] = -0.5;
1176    coords[1] = -0.5;
1177    coords[2] =  0.0;
1178    break;
1179  case  8:
1180    coords[0] =  0.5;
1181    coords[1] = -0.5;
1182    coords[2] =  0.0;
1183    break;
1184  case  9:
1185    coords[0] =  0.5;
1186    coords[1] =  0.0;
1187    coords[2] =  0.5;
1188    break;
1189  case 10:
1190    coords[0] =  0.0;
1191    coords[1] =  0.5;
1192    coords[2] =  0.5;
1193    break;
1194  case 11:
1195    coords[0] = -0.5;
1196    coords[1] =  0.0;
1197    coords[2] =  0.5;
1198    break;
1199  case 12:
1200    coords[0] =  0.0;
1201    coords[1] = -0.5;
1202    coords[2] =  0.5;
1203    break;
1204    LOCAL_COORD_MACRO_END;
1205
1206    SHAPE_FUN_MACRO_BEGIN;
1207    funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1208      (gc[0] - 0.5)/(1.0 - gc[2]);
1209    funValue[1] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1210      (gc[1] - 0.5)/(1.0 - gc[2]);
1211    funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1212      (-gc[0] - 0.5)/(1.0 - gc[2]);
1213    funValue[3] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1214      (-gc[1] - 0.5)/(1.0 - gc[2]);
1215
1216    funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1217
1218    funValue[5] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1219      (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1220    funValue[6] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1221      (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1222    funValue[7] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1223      (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1224    funValue[8] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1225      (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1226
1227    funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1228      (1.0 - gc[2]);
1229    funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1230      (1.0 - gc[2]);
1231    funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1232      (1.0 - gc[2]);
1233    funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1234      (1.0 - gc[2]);
1235    SHAPE_FUN_MACRO_END;
1236 }
1237
1238 /*!
1239  * Init Quadratic Pyramid Reference coordinates ans Shape function.
1240  * Case B.
1241  */
1242 void GaussInfo::pyra13bInit()
1243 {
1244   LOCAL_COORD_MACRO_BEGIN;
1245  case  0:
1246    coords[0] =  1.0;
1247    coords[1] =  0.0;
1248    coords[2] =  0.0;
1249    break;
1250  case  3:
1251    coords[0] =  0.0;
1252    coords[1] =  1.0;
1253    coords[2] =  0.0;
1254    break;
1255  case  2:
1256    coords[0] = -1.0;
1257    coords[1] =  0.0;
1258    coords[2] =  0.0;
1259    break;
1260  case  1:
1261    coords[0] =  0.0;
1262    coords[1] = -1.0;
1263    coords[2] =  0.0;
1264    break;
1265  case  4:
1266    coords[0] =  0.0;
1267    coords[1] =  0.0;
1268    coords[2] =  1.0;
1269    break;
1270  case  8:
1271    coords[0] =  0.5;
1272    coords[1] =  0.5;
1273    coords[2] =  0.0;
1274    break;
1275  case  7:
1276    coords[0] = -0.5;
1277    coords[1] =  0.5;
1278    coords[2] =  0.0;
1279    break;
1280  case  6:
1281    coords[0] = -0.5;
1282    coords[1] = -0.5;
1283    coords[2] =  0.0;
1284    break;
1285  case  5:
1286    coords[0] =  0.5;
1287    coords[1] = -0.5;
1288    coords[2] =  0.0;
1289    break;
1290  case  9:
1291    coords[0] =  0.5;
1292    coords[1] =  0.0;
1293    coords[2] =  0.5;
1294    break;
1295  case 12:
1296    coords[0] =  0.0;
1297    coords[1] =  0.5;
1298    coords[2] =  0.5;
1299    break;
1300  case 11:
1301    coords[0] = -0.5;
1302    coords[1] =  0.0;
1303    coords[2] =  0.5;
1304    break;
1305  case 10:
1306    coords[0] =  0.0;
1307    coords[1] = -0.5;
1308    coords[2] =  0.5;
1309    break;
1310    LOCAL_COORD_MACRO_END;
1311
1312    SHAPE_FUN_MACRO_BEGIN;
1313    funValue[0] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1314      (gc[0] - 0.5)/(1.0 - gc[2]);
1315    funValue[3] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] - gc[1] + gc[2] - 1.0)*
1316      (gc[1] - 0.5)/(1.0 - gc[2]);
1317    funValue[2] = 0.5*(+gc[0] - gc[1] + gc[2] - 1.0)*(+gc[0] + gc[1] + gc[2] - 1.0)*
1318      (-gc[0] - 0.5)/(1.0 - gc[2]);
1319    funValue[1] = 0.5*(+gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1320      (-gc[1] - 0.5)/(1.0 - gc[2]);
1321
1322    funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
1323
1324    funValue[8] = 0.5*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)*
1325      (gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1326    funValue[7] = 0.5*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)*
1327      (gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1328    funValue[6] = 0.5*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)*
1329      (-gc[0] + gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1330    funValue[5] = 0.5*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)*
1331      (-gc[0] - gc[1] + gc[2] - 1.0)/(1.0 - gc[2]);
1332
1333    funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
1334      (1.0 - gc[2]);
1335    funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
1336      (1.0 - gc[2]);
1337    funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
1338      (1.0 - gc[2]);
1339    funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
1340      (1.0 - gc[2]);
1341    SHAPE_FUN_MACRO_END;
1342 }
1343
1344
1345 /*!
1346  * Init Pentahedron Reference coordinates and Shape function.
1347  * Case A.
1348  */
1349 void GaussInfo::penta6aInit()
1350 {
1351   LOCAL_COORD_MACRO_BEGIN;
1352  case  0:
1353    coords[0] = -1.0;
1354    coords[1] =  1.0;
1355    coords[2] =  0.0;
1356    break;
1357  case  1:
1358    coords[0] = -1.0;
1359    coords[1] = -0.0;
1360    coords[2] =  1.0;
1361    break;
1362  case  2:
1363    coords[0] = -1.0;
1364    coords[1] =  0.0;
1365    coords[2] =  0.0;
1366    break;
1367  case  3:
1368    coords[0] =  1.0;
1369    coords[1] =  1.0;
1370    coords[2] =  0.0;
1371    break;
1372  case  4:
1373    coords[0] =  1.0;
1374    coords[1] =  0.0;
1375    coords[2] =  1.0;
1376    break;
1377  case  5:
1378    coords[0] =  1.0;
1379    coords[1] =  0.0;
1380    coords[2] =  0.0;
1381    break;
1382    LOCAL_COORD_MACRO_END;
1383
1384    SHAPE_FUN_MACRO_BEGIN;
1385    funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1386    funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
1387    funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1388
1389    funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1390    funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
1391    funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1392    SHAPE_FUN_MACRO_END;
1393 }
1394
1395 /*!
1396  * Init Pentahedron Reference coordinates and Shape function.
1397  * Case B.
1398  */
1399 void GaussInfo::penta6bInit()
1400 {
1401   LOCAL_COORD_MACRO_BEGIN;
1402  case  0:
1403    coords[0] = -1.0;
1404    coords[1] =  1.0;
1405    coords[2] =  0.0;
1406    break;
1407  case  2:
1408    coords[0] = -1.0;
1409    coords[1] = -0.0;
1410    coords[2] =  1.0;
1411    break;
1412  case  1:
1413    coords[0] = -1.0;
1414    coords[1] =  0.0;
1415    coords[2] =  0.0;
1416    break;
1417  case  3:
1418    coords[0] =  1.0;
1419    coords[1] =  1.0;
1420    coords[2] =  0.0;
1421    break;
1422  case  5:
1423    coords[0] =  1.0;
1424    coords[1] =  0.0;
1425    coords[2] =  1.0;
1426    break;
1427  case  4:
1428    coords[0] =  1.0;
1429    coords[1] =  0.0;
1430    coords[2] =  0.0;
1431    break;
1432    LOCAL_COORD_MACRO_END;
1433
1434    SHAPE_FUN_MACRO_BEGIN;
1435    funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
1436    funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
1437    funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1438    funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
1439    funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
1440    funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1441    SHAPE_FUN_MACRO_END;
1442 }
1443 /*!
1444  * Init Pentahedron Reference coordinates and Shape function.
1445  * Case A.
1446  */
1447 void GaussInfo::penta15aInit()
1448 {
1449   LOCAL_COORD_MACRO_BEGIN;
1450  case  0:
1451    coords[0] = -1.0;
1452    coords[1] =  1.0;
1453    coords[2] =  0.0;
1454    break;
1455  case  1:
1456    coords[0] = -1.0;
1457    coords[1] = -0.0;
1458    coords[2] =  1.0;
1459    break;
1460  case  2:
1461    coords[0] = -1.0;
1462    coords[1] =  0.0;
1463    coords[2] =  0.0;
1464    break;
1465  case  3:
1466    coords[0] =  1.0;
1467    coords[1] =  1.0;
1468    coords[2] =  0.0;
1469    break;
1470  case  4:
1471    coords[0] =  1.0;
1472    coords[1] =  0.0;
1473    coords[2] =  1.0;
1474    break;
1475  case  5:
1476    coords[0] =  1.0;
1477    coords[1] =  0.0;
1478    coords[2] =  0.0;
1479    break;
1480
1481  case  6:
1482    coords[0] = -1.0;
1483    coords[1] =  0.5;
1484    coords[2] =  0.5;
1485    break;
1486  case  7:
1487    coords[0] = -1.0;
1488    coords[1] =  0.0;
1489    coords[2] =  0.5;
1490    break;
1491  case  8:
1492    coords[0] = -1.0;
1493    coords[1] =  0.5;
1494    coords[2] =  0.0;
1495    break;
1496  case  9:
1497    coords[0] =  0.0;
1498    coords[1] =  1.0;
1499    coords[2] =  0.0;
1500    break;
1501  case 10:
1502    coords[0] =  0.0;
1503    coords[1] =  0.0;
1504    coords[2] =  1.0;
1505    break;
1506  case 11:
1507    coords[0] =  0.0;
1508    coords[1] =  0.0;
1509    coords[2] =  0.0;
1510    break;
1511  case 12:
1512    coords[0] =  1.0;
1513    coords[1] =  0.5;
1514    coords[2] =  0.5;
1515    break;
1516  case 13:
1517    coords[0] =  1.0;
1518    coords[1] =  0.0;
1519    coords[2] =  0.5;
1520    break;
1521  case 14:
1522    coords[0] =  1.0;
1523    coords[1] =  0.5;
1524    coords[2] =  0.0;
1525    break;
1526    LOCAL_COORD_MACRO_END;
1527
1528    SHAPE_FUN_MACRO_BEGIN;
1529    funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1530    funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1531    funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1532
1533    funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1534    funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1535    funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1536
1537    funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1538    funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1539    funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1540
1541    funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
1542    funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
1543    funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1544
1545    funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1546    funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1547    funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1548    SHAPE_FUN_MACRO_END;
1549 }
1550
1551 /*!
1552  * Init Qaudratic Pentahedron Reference coordinates and Shape function.
1553  * Case B.
1554  */
1555 void GaussInfo::penta15bInit()
1556 {
1557   LOCAL_COORD_MACRO_BEGIN;
1558  case  0:
1559    coords[0] = -1.0;
1560    coords[1] =  1.0;
1561    coords[2] =  0.0;
1562    break;
1563  case  2:
1564    coords[0] = -1.0;
1565    coords[1] = -0.0;
1566    coords[2] =  1.0;
1567    break;
1568  case  1:
1569    coords[0] = -1.0;
1570    coords[1] =  0.0;
1571    coords[2] =  0.0;
1572    break;
1573  case  3:
1574    coords[0] =  1.0;
1575    coords[1] =  1.0;
1576    coords[2] =  0.0;
1577    break;
1578  case  5:
1579    coords[0] =  1.0;
1580    coords[1] =  0.0;
1581    coords[2] =  1.0;
1582    break;
1583  case  4:
1584    coords[0] =  1.0;
1585    coords[1] =  0.0;
1586    coords[2] =  0.0;
1587    break;
1588
1589  case  8:
1590    coords[0] = -1.0;
1591    coords[1] =  0.5;
1592    coords[2] =  0.5;
1593    break;
1594  case  7:
1595    coords[0] = -1.0;
1596    coords[1] =  0.0;
1597    coords[2] =  0.5;
1598    break;
1599  case  6:
1600    coords[0] = -1.0;
1601    coords[1] =  0.5;
1602    coords[2] =  0.0;
1603    break;
1604  case 12:
1605    coords[0] =  0.0;
1606    coords[1] =  1.0;
1607    coords[2] =  0.0;
1608    break;
1609  case 14:
1610    coords[0] =  0.0;
1611    coords[1] =  0.0;
1612    coords[2] =  1.0;
1613    break;
1614  case 13:
1615    coords[0] =  0.0;
1616    coords[1] =  0.0;
1617    coords[2] =  0.0;
1618    break;
1619  case 11:
1620    coords[0] =  1.0;
1621    coords[1] =  0.5;
1622    coords[2] =  0.5;
1623    break;
1624  case 10:
1625    coords[0] =  1.0;
1626    coords[1] =  0.0;
1627    coords[2] =  0.5;
1628    break;
1629  case  9:
1630    coords[0] =  1.0;
1631    coords[1] =  0.5;
1632    coords[2] =  0.0;
1633    break;
1634    LOCAL_COORD_MACRO_END;
1635
1636    SHAPE_FUN_MACRO_BEGIN;
1637    funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
1638    funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
1639    funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1640
1641    funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
1642    funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
1643    funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
1644
1645    funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
1646    funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1647    funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
1648
1649    funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
1650    funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
1651    funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
1652
1653    funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
1654    funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1655    funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
1656    SHAPE_FUN_MACRO_END;
1657 }
1658
1659 /*!
1660  * Init Hehahedron Reference coordinates and Shape function.
1661  * Case A.
1662  */
1663 void GaussInfo::hexa8aInit()
1664 {
1665   LOCAL_COORD_MACRO_BEGIN;
1666  case  0:
1667    coords[0] = -1.0;
1668    coords[1] = -1.0;
1669    coords[2] = -1.0;
1670    break;
1671  case  1:
1672    coords[0] =  1.0;
1673    coords[1] = -1.0;
1674    coords[2] = -1.0;
1675    break;
1676  case  2:
1677    coords[0] =  1.0;
1678    coords[1] =  1.0;
1679    coords[2] = -1.0;
1680    break;
1681  case  3:
1682    coords[0] = -1.0;
1683    coords[1] =  1.0;
1684    coords[2] = -1.0;
1685    break;
1686  case  4:
1687    coords[0] = -1.0;
1688    coords[1] = -1.0;
1689    coords[2] =  1.0;
1690    break;
1691  case  5:
1692    coords[0] =  1.0;
1693    coords[1] = -1.0;
1694    coords[2] =  1.0;
1695    break;
1696  case  6:
1697    coords[0] =  1.0;
1698    coords[1] =  1.0;
1699    coords[2] =  1.0;
1700    break;
1701  case  7:
1702    coords[0] = -1.0;
1703    coords[1] =  1.0;
1704    coords[2] =  1.0;
1705    break;
1706    LOCAL_COORD_MACRO_END;
1707
1708    SHAPE_FUN_MACRO_BEGIN;
1709    funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1710    funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1711    funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1712    funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1713
1714    funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1715    funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1716    funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1717    funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1718    SHAPE_FUN_MACRO_END;
1719 }
1720
1721 /*!
1722  * Init Hehahedron Reference coordinates and Shape function.
1723  * Case B.
1724  */
1725 void GaussInfo::hexa8bInit()
1726 {
1727   LOCAL_COORD_MACRO_BEGIN;
1728  case  0:
1729    coords[0] = -1.0;
1730    coords[1] = -1.0;
1731    coords[2] = -1.0;
1732    break;
1733  case  3:
1734    coords[0] =  1.0;
1735    coords[1] = -1.0;
1736    coords[2] = -1.0;
1737    break;
1738  case  2:
1739    coords[0] =  1.0;
1740    coords[1] =  1.0;
1741    coords[2] = -1.0;
1742    break;
1743  case  1:
1744    coords[0] = -1.0;
1745    coords[1] =  1.0;
1746    coords[2] = -1.0;
1747    break;
1748  case  4:
1749    coords[0] = -1.0;
1750    coords[1] = -1.0;
1751    coords[2] =  1.0;
1752    break;
1753  case  7:
1754    coords[0] =  1.0;
1755    coords[1] = -1.0;
1756    coords[2] =  1.0;
1757    break;
1758  case  6:
1759    coords[0] =  1.0;
1760    coords[1] =  1.0;
1761    coords[2] =  1.0;
1762    break;
1763  case  5:
1764    coords[0] = -1.0;
1765    coords[1] =  1.0;
1766    coords[2] =  1.0;
1767    break;
1768    LOCAL_COORD_MACRO_END;
1769
1770    SHAPE_FUN_MACRO_BEGIN;
1771    funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1772    funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1773    funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1774    funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1775
1776    funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1777    funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1778    funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1779    funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1780    SHAPE_FUN_MACRO_END;
1781 }
1782
1783 /*!
1784  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1785  * Case A.
1786  */
1787 void GaussInfo::hexa20aInit()
1788 {
1789   LOCAL_COORD_MACRO_BEGIN;
1790  case  0:
1791    coords[0] = -1.0;
1792    coords[1] = -1.0;
1793    coords[2] = -1.0;
1794    break;
1795  case  1:
1796    coords[0] =  1.0;
1797    coords[1] = -1.0;
1798    coords[2] = -1.0;
1799    break;
1800  case  2:
1801    coords[0] =  1.0;
1802    coords[1] =  1.0;
1803    coords[2] = -1.0;
1804    break;
1805  case  3:
1806    coords[0] = -1.0;
1807    coords[1] =  1.0;
1808    coords[2] = -1.0;
1809    break;
1810  case  4:
1811    coords[0] = -1.0;
1812    coords[1] = -1.0;
1813    coords[2] =  1.0;
1814    break;
1815  case  5:
1816    coords[0] =  1.0;
1817    coords[1] = -1.0;
1818    coords[2] =  1.0;
1819    break;
1820  case  6:
1821    coords[0] =  1.0;
1822    coords[1] =  1.0;
1823    coords[2] =  1.0;
1824    break;
1825  case  7:
1826    coords[0] = -1.0;
1827    coords[1] =  1.0;
1828    coords[2] =  1.0;
1829    break;
1830
1831  case  8:
1832    coords[0] =  0.0;
1833    coords[1] = -1.0;
1834    coords[2] = -1.0;
1835    break;
1836  case  9:
1837    coords[0] =  1.0;
1838    coords[1] =  0.0;
1839    coords[2] = -1.0;
1840    break;
1841  case 10:
1842    coords[0] =  0.0;
1843    coords[1] =  1.0;
1844    coords[2] = -1.0;
1845    break;
1846  case 11:
1847    coords[0] = -1.0;
1848    coords[1] =  0.0;
1849    coords[2] = -1.0;
1850    break;
1851  case 12:
1852    coords[0] = -1.0;
1853    coords[1] = -1.0;
1854    coords[2] =  0.0;
1855    break;
1856  case 13:
1857    coords[0] =  1.0;
1858    coords[1] = -1.0;
1859    coords[2] =  0.0;
1860    break;
1861  case 14:
1862    coords[0] =  1.0;
1863    coords[1] =  1.0;
1864    coords[2] =  0.0;
1865    break;
1866  case 15:
1867    coords[0] = -1.0;
1868    coords[1] =  1.0;
1869    coords[2] =  0.0;
1870    break;
1871  case 16:
1872    coords[0] =  0.0;
1873    coords[1] = -1.0;
1874    coords[2] =  1.0;
1875    break;
1876  case 17:
1877    coords[0] =  1.0;
1878    coords[1] =  0.0;
1879    coords[2] =  1.0;
1880    break;
1881  case 18:
1882    coords[0] =  0.0;
1883    coords[1] =  1.0;
1884    coords[2] =  1.0;
1885    break;
1886  case 19:
1887    coords[0] = -1.0;
1888    coords[1] =  0.0;
1889    coords[2] =  1.0;
1890    break;
1891    LOCAL_COORD_MACRO_END;
1892
1893    SHAPE_FUN_MACRO_BEGIN;
1894    funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1895      (-2.0 - gc[0] - gc[1] - gc[2]);
1896    funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
1897      (-2.0 + gc[0] - gc[1] - gc[2]);
1898    funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1899      (-2.0 + gc[0] + gc[1] - gc[2]);
1900    funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
1901      (-2.0 - gc[0] + gc[1] - gc[2]);
1902    funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1903      (-2.0 - gc[0] - gc[1] + gc[2]);
1904    funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
1905      (-2.0 + gc[0] - gc[1] + gc[2]);
1906    funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1907      (-2.0 + gc[0] + gc[1] + gc[2]);
1908    funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
1909      (-2.0 - gc[0] + gc[1] + gc[2]);
1910
1911    funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
1912    funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
1913    funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
1914    funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
1915    funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
1916    funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
1917    funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
1918    funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
1919    funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
1920    funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
1921    funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
1922    funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
1923    SHAPE_FUN_MACRO_END;
1924 }
1925
1926 /*!
1927  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
1928  * Case B.
1929  */
1930 void GaussInfo::hexa20bInit()
1931 {
1932   LOCAL_COORD_MACRO_BEGIN;
1933  case  0:
1934    coords[0] = -1.0;
1935    coords[1] = -1.0;
1936    coords[2] = -1.0;
1937    break;
1938  case  3:
1939    coords[0] =  1.0;
1940    coords[1] = -1.0;
1941    coords[2] = -1.0;
1942    break;
1943  case  2:
1944    coords[0] =  1.0;
1945    coords[1] =  1.0;
1946    coords[2] = -1.0;
1947    break;
1948  case  1:
1949    coords[0] = -1.0;
1950    coords[1] =  1.0;
1951    coords[2] = -1.0;
1952    break;
1953  case  4:
1954    coords[0] = -1.0;
1955    coords[1] = -1.0;
1956    coords[2] =  1.0;
1957    break;
1958  case  7:
1959    coords[0] =  1.0;
1960    coords[1] = -1.0;
1961    coords[2] =  1.0;
1962    break;
1963  case  6:
1964    coords[0] =  1.0;
1965    coords[1] =  1.0;
1966    coords[2] =  1.0;
1967    break;
1968  case  5:
1969    coords[0] = -1.0;
1970    coords[1] =  1.0;
1971    coords[2] =  1.0;
1972    break;
1973
1974  case 11:
1975    coords[0] =  0.0;
1976    coords[1] = -1.0;
1977    coords[2] = -1.0;
1978    break;
1979  case 10:
1980    coords[0] =  1.0;
1981    coords[1] =  0.0;
1982    coords[2] = -1.0;
1983    break;
1984  case  9:
1985    coords[0] =  0.0;
1986    coords[1] =  1.0;
1987    coords[2] = -1.0;
1988    break;
1989  case  8:
1990    coords[0] = -1.0;
1991    coords[1] =  0.0;
1992    coords[2] = -1.0;
1993    break;
1994  case 16:
1995    coords[0] = -1.0;
1996    coords[1] = -1.0;
1997    coords[2] =  0.0;
1998    break;
1999  case 19:
2000    coords[0] =  1.0;
2001    coords[1] = -1.0;
2002    coords[2] =  0.0;
2003    break;
2004  case 18:
2005    coords[0] =  1.0;
2006    coords[1] =  1.0;
2007    coords[2] =  0.0;
2008    break;
2009  case 17:
2010    coords[0] = -1.0;
2011    coords[1] =  1.0;
2012    coords[2] =  0.0;
2013    break;
2014  case 15:
2015    coords[0] =  0.0;
2016    coords[1] = -1.0;
2017    coords[2] =  1.0;
2018    break;
2019  case 14:
2020    coords[0] =  1.0;
2021    coords[1] =  0.0;
2022    coords[2] =  1.0;
2023    break;
2024  case 13:
2025    coords[0] =  0.0;
2026    coords[1] =  1.0;
2027    coords[2] =  1.0;
2028    break;
2029  case 12:
2030    coords[0] = -1.0;
2031    coords[1] =  0.0;
2032    coords[2] =  1.0;
2033    break;
2034    LOCAL_COORD_MACRO_END;
2035
2036    SHAPE_FUN_MACRO_BEGIN;
2037
2038    funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2039      (-2.0 - gc[0] - gc[1] - gc[2]);
2040    funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
2041      (-2.0 + gc[0] - gc[1] - gc[2]);
2042    funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2043      (-2.0 + gc[0] + gc[1] - gc[2]);
2044    funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
2045      (-2.0 - gc[0] + gc[1] - gc[2]);
2046    funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2047      (-2.0 - gc[0] - gc[1] + gc[2]);
2048    funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
2049      (-2.0 + gc[0] - gc[1] + gc[2]);
2050    funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2051      (-2.0 + gc[0] + gc[1] + gc[2]);
2052    funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
2053      (-2.0 - gc[0] + gc[1] + gc[2]);
2054
2055    funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
2056    funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
2057    funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
2058    funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
2059    funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
2060    funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
2061    funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
2062    funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
2063    funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
2064    funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
2065    funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
2066    funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
2067    SHAPE_FUN_MACRO_END;
2068 }
2069
2070
2071
2072 ////////////////////////////////////////////////////////////////////////////////////////////////
2073 //                                GAUSS COORD CLASS                                           //
2074 ////////////////////////////////////////////////////////////////////////////////////////////////
2075 /*!
2076  * Constructor
2077  */
2078 GaussCoords::GaussCoords()
2079 {
2080 }
2081
2082 /*!
2083  * Destructor
2084  */
2085 GaussCoords::~GaussCoords()
2086 {
2087   GaussInfoVector::iterator it = _my_gauss_info.begin();
2088   for( ; it != _my_gauss_info.end(); it++ )
2089     {
2090       if((*it) != NULL)
2091         delete (*it);
2092     }
2093 }
2094
2095 /*!
2096  * Add Gauss localization info
2097  */
2098 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
2099                                 int coordDim,
2100                                 const double* theGaussCoord,
2101                                 int theNbGauss,
2102                                 const double* theReferenceCoord,
2103                                 int theNbRef)
2104 {
2105   GaussInfoVector::iterator it = _my_gauss_info.begin();
2106   for( ; it != _my_gauss_info.end(); it++ )
2107     {
2108       if( (*it)->getCellType() == theGeometry )
2109         {
2110           break;
2111         }
2112     }
2113
2114   DataVector aGaussCoord;
2115   for(int i = 0 ; i < theNbGauss*coordDim; i++ )
2116     aGaussCoord.push_back(theGaussCoord[i]);
2117
2118   DataVector aReferenceCoord;
2119   for(int i = 0 ; i < theNbRef*coordDim; i++ )
2120     aReferenceCoord.push_back(theReferenceCoord[i]);
2121
2122
2123   GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
2124   info->initLocalInfo();
2125
2126   //If info with cell type doesn't exist add it
2127   if( it == _my_gauss_info.end() )
2128     {
2129       _my_gauss_info.push_back(info);
2130
2131       // If information exists, update it
2132     }
2133   else
2134     {
2135       int index = std::distance(_my_gauss_info.begin(),it);
2136       delete (*it);
2137       _my_gauss_info[index] = info;
2138     }
2139 }
2140
2141
2142 /*!
2143  * Calculate gauss points coordinates
2144  */
2145 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry,
2146                                       const double *theNodeCoords,
2147                                       const int theSpaceDim,
2148                                       const int *theIndex)
2149 {
2150   const GaussInfo *info = getInfoGivenCellType(theGeometry);
2151   int nbCoords = theSpaceDim * info->getNbGauss();
2152   double *aCoords = new double[nbCoords];
2153   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
2154   return aCoords;
2155 }
2156
2157
2158 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2159 {
2160   const GaussInfo *info = getInfoGivenCellType(theGeometry);
2161   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
2162 }
2163
2164 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const int *theIndex, double *result)
2165 {
2166   int aConn = info->getNbRef();
2167
2168   int nbCoords = theSpaceDim * info->getNbGauss();
2169   std::fill(result,result+nbCoords,0.);
2170
2171   for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ )
2172     {
2173       double *coord=result+gaussId*theSpaceDim;
2174       const double *function=info->getFunctionValues(gaussId);
2175       for ( int connId = 0; connId < aConn ; connId++ )
2176         {
2177           const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
2178           for( int dimId = 0; dimId < theSpaceDim; dimId++ )
2179             coord[dimId] += nodeCoord[dimId]*function[connId];
2180         }
2181     }
2182 }
2183
2184 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
2185 {
2186   GaussInfoVector::const_iterator it = _my_gauss_info.begin();
2187   //Try to find gauss localization info
2188   for( ; it != _my_gauss_info.end() ; it++ )
2189     if( (*it)->getCellType()==cellType)
2190       return (*it);
2191   throw INTERP_KERNEL::Exception("Can't find gauss localization information !");
2192 }