src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx

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