src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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