src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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