From: ageay <ageay>
Date: Wed, 23 Feb 2011 09:59:01 +0000 (+0000)
Subject: Thank you to Edward for its wonderful cooperation into code standardization.
X-Git-Tag: V6_main_FINAL~1054
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=2848eccd2c3064dd7b61091ccb3e11ec06370be8;p=tools%2Fmedcoupling.git

Thank you to Edward for its wonderful cooperation into code standardization.
---

diff --git a/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx
deleted file mode 100644
index 51ec9c344..000000000
--- a/src/INTERP_KERNEL/GaussPoints/GaussCoords.cxx
+++ /dev/null
@@ -1,2162 +0,0 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D
-//
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
-//
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
-//
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
-//
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-
-//Local includes
-#include "GaussCoords.hxx"
-#include "CellModel.hxx"
-using namespace INTERP_KERNEL;
-
-//STL includes
-#include <math.h>
-#include <algorithm>
-#include <sstream>
-
-//#define MYDEBUG
-
-#ifdef MYDEBUG
-#include <iostream>
-#endif
-
-//Define common part of the code in the MACRO
-//---------------------------------------------------------------
-#define LOCAL_COORD_MACRO_BEGIN                                         \
-  myLocalReferenceCoord.resize( myLocalRefDim*myLocalNbRef );           \
-  for( int refId = 0; refId < myLocalNbRef; refId++ )                   \
-    {                                                                   \
-      double* coords = &myLocalReferenceCoord[ refId*myLocalRefDim ];   \
-      switch(refId)                                                     \
-        {
-
-//---------------------------------------------------------------
-#define LOCAL_COORD_MACRO_END                   \
-  }                                             \
-}
-
-//---------------------------------------------------------------
-#define SHAPE_FUN_MACRO_BEGIN                                           \
-  for( int gaussId     = 0 ; gaussId < myNbGauss ; gaussId++ )          \
-    {                                                                   \
-      double* funValue =  &myFunctionValue[ gaussId * myNbRef ];        \
-      const double* gc = &myGaussCoord[ gaussId * GetGaussCoordDim() ];
-
-//---------------------------------------------------------------
-#define SHAPE_FUN_MACRO_END                     \
-  }
-
-#define CHECK_MACRO                                                        \
-  if( ! aSatify )                                                          \
-    {                                                                      \
-      std::ostringstream stream;                                           \
-      stream << "Error in the gauss localization for the cell with type "; \
-      stream << cellModel.getRepr();                                       \
-      stream << " !!!";                                                    \
-      throw INTERP_KERNEL::Exception(stream.str().c_str());                \
-    }
-
-
-//---------------------------------------------------------------
-static bool IsEqual(double theLeft, double theRight) 
-{
-  static double EPS = 1.0E-3;
-  if(fabs(theLeft) + fabs(theRight) > EPS)
-    return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
-  return true;
-}
-
-
-////////////////////////////////////////////////////////////////////////////////////////////////
-//                                GAUSS INFO CLASS                                            //
-////////////////////////////////////////////////////////////////////////////////////////////////
-
-/*!
- * Constructor of the GaussInfo
- */
-GaussInfo::GaussInfo( NormalizedCellType theGeometry,
-                      const DataVector& theGaussCoord,
-                      int theNbGauss,
-                      const DataVector& theReferenceCoord,
-                      int theNbRef ) :
-  myGeometry(theGeometry),
-  myNbGauss(theNbGauss),
-  myGaussCoord(theGaussCoord),
-  myNbRef(theNbRef),
-  myReferenceCoord(theReferenceCoord)
-{
-
-  //Allocate shape function values
-  myFunctionValue.resize( myNbGauss * myNbRef );
-}
-
-/*!
- * Destructor
- */
-GaussInfo::~GaussInfo()
-{
-}
-
-/*!
- * Return dimension of the gauss coordinates
- */
-int GaussInfo::GetGaussCoordDim() const 
-{
-  if( myNbGauss ) 
-    {
-      return myGaussCoord.size()/myNbGauss;
-    }
-  else 
-    {
-      return 0;
-    }
-}
-
-/*!
- * Return dimension of the reference coordinates
- */
-int GaussInfo::GetReferenceCoordDim() const 
-{
-  if( myNbRef ) 
-    {
-      return myReferenceCoord.size()/myNbRef;
-    }
-  else 
-    {
-      return 0;
-    }
-}
-
-/*!
- * Return type of the cell.
- */
-NormalizedCellType GaussInfo::GetCellType() const 
-{
-  return myGeometry;
-}
-
-/*!
- * Return Nb of the gauss points.
- */
-int GaussInfo::GetNbGauss() const 
-{
-  return myNbGauss;
-}
-
-/*!
- * Return Nb of the reference coordinates.
- */
-int GaussInfo::GetNbRef() const 
-{
-  return myNbRef;
-}
-
-/*!
- * Check coordinates
- */
-bool GaussInfo::isSatisfy() 
-{
-
-  bool anIsSatisfy = ((myLocalNbRef == myNbRef) && (myLocalRefDim == GetReferenceCoordDim()));
-  //Check coordinates
-  if(anIsSatisfy)
-    {
-      for( int refId = 0; refId < myLocalNbRef; refId++ ) 
-        {
-          double* refCoord = &myReferenceCoord[ refId*myLocalRefDim ];
-          double* localRefCoord = &myLocalReferenceCoord[ refId*myLocalRefDim ];
-          bool anIsEqual = false;
-          for( int dimId = 0; dimId < myLocalRefDim; dimId++ ) 
-            {
-              anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
-              if(!anIsEqual ) 
-                {
-                  return false;
-                }
-            }
-        }
-    }
-  return anIsSatisfy;
-}
-
-/*!
- * Initialize the internal vectors
- */
-void GaussInfo::InitLocalInfo() throw (INTERP_KERNEL::Exception) 
-{
-  bool aSatify = false;
-  const CellModel& cellModel=CellModel::getCellModel(myGeometry);
-  switch( myGeometry ) 
-    {
-    case NORM_SEG2:
-      myLocalRefDim = 1;
-      myLocalNbRef  = 2;
-      Seg2Init();
-      aSatify = isSatisfy();
-      CHECK_MACRO;
-      break;
-
-    case NORM_SEG3:
-      myLocalRefDim = 1;
-      myLocalNbRef  = 3;
-      Seg3Init();
-      aSatify = isSatisfy();
-      CHECK_MACRO;
-      break;
-
-    case NORM_TRI3:
-      myLocalRefDim = 2;
-      myLocalNbRef  = 3;
-      Tria3aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Tria3bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_TRI6:
-      myLocalRefDim = 2;
-      myLocalNbRef  = 6;
-      Tria6aInit();
-      aSatify = isSatisfy();
-      if(!aSatify)
-        {
-          Tria6bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_QUAD4:
-      myLocalRefDim = 2;
-      myLocalNbRef  = 4;
-      Quad4aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Quad4bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_QUAD8:
-      myLocalRefDim = 2;
-      myLocalNbRef  = 8;
-      Quad8aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Quad8bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_TETRA4:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 4;
-      Tetra4aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Tetra4bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_TETRA10:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 10;
-      Tetra10aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Tetra10bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_PYRA5:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 5;
-      Pyra5aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Pyra5bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_PYRA13:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 13;
-      Pyra13aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Pyra13bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_PENTA6:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 6;
-      Penta6aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Penta6bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_PENTA15:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 15;
-      Penta15aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Penta15bInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_HEXA8:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 8;
-      Hexa8aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Hexa8aInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    case NORM_HEXA20:
-      myLocalRefDim = 3;
-      myLocalNbRef  = 20;
-      Hexa20aInit();
-      aSatify = isSatisfy();
-
-      if(!aSatify)
-        {
-          Hexa20aInit();
-          aSatify = isSatisfy();
-          CHECK_MACRO;
-        }
-      break;
-
-    default:
-      throw INTERP_KERNEL::Exception("Bad Cell type !!!");
-      break;
-    }
-}
-
-/**
- * Return shape function value by node id
- */
-const double* GaussInfo::GetFunctionValues( const int theGaussId ) const 
-{
-  return &myFunctionValue[ myNbRef*theGaussId ];
-}
-
-/*!
- * Init Segment 2 Reference coordinates ans Shape function.
- */
-void GaussInfo::Seg2Init() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
-  case  0:
-    coords[0] = -1.0;
-    break;
-  case  1:
-    coords[0] =  1.0;
-    break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*(1.0 - gc[0]);
-   funValue[1] = 0.5*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Segment 3 Reference coordinates ans Shape function.
- */
-void GaussInfo::Seg3Init() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
-  case  0:
-    coords[0] = -1.0;
-    break;
-  case  1:
-    coords[0] =  1.0;
-    break;
-  case  2:
-    coords[0] =  0.0;
-    break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*(1.0 - gc[0])*gc[0];
-   funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
-   funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Triangle Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Tria3aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*(1.0 + gc[1]);
-   funValue[1] = -0.5*(gc[0] + gc[1]);
-   funValue[2] = 0.5*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Triangle Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Tria3bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   break;
- case  1:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   break;
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 1.0 - gc[0] - gc[1];
-   funValue[1] = gc[0];
-   funValue[2] = gc[1];
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Triangle Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Tria6aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
- case  3:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   break;
- case  5:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
-   funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
-   funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
-   funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
-   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
-   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Triangle Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Tria6bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   break;
-
- case  1:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   break;
-
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   break;
-
- case  3:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   break;
-
- case  4:
-   coords[0] =  0.5;
-   coords[1] =  0.5;
-   break;
-
- case  5:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   break;
-
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
-   funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
-   funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
-   funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
-   funValue[4] = 4.0*gc[0]*gc[1];
-   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadrangle Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Quad4aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   break;
-
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
-   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
-   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
-   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadrangle Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Quad4bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  1:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   break;
- case  3:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
-
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
-   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
-   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
-   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
-   SHAPE_FUN_MACRO_END;
-}
-
-
-/*!
- * Init Quadratic Quadrangle Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Quad8aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   break;
- case  4:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   break;
- case  5:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   break;
- case  6:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   break;
- case  7:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
-   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
-   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
-   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
-   funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
-   funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
-   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
-   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Quadrangle Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Quad8bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   break;
- case  1:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   break;
- case  3:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   break;
- case  6:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   break;
- case  7:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
-   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
-   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
-   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
-   funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
-   funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
-   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
-   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Tetrahedron Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Tetra4aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = gc[1];
-   funValue[1] = gc[2];
-   funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
-   funValue[3] = gc[0];
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Tetrahedron Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Tetra4bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = gc[1];
-   funValue[2] = gc[2];
-   funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
-   funValue[3] = gc[0];
-   SHAPE_FUN_MACRO_END;
-
-}
-
-/*!
- * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Tetra10aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case  5:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  6:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  7:
-   coords[0] =  0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  8:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  9:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
-   funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
-   funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
-   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
-   funValue[4] = 4.0*gc[1]*gc[2];
-   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
-   funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
-   funValue[7] = 4.0*gc[0]*gc[1];
-   funValue[8] = 4.0*gc[0]*gc[2];
-   funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Tetra10bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  6:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case  5:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  7:
-   coords[0] =  0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  9:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  8:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
-   funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
-   funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
-   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
-   funValue[6] = 4.0*gc[1]*gc[2];
-   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
-   funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
-   funValue[7] = 4.0*gc[0]*gc[1];
-   funValue[9] = 4.0*gc[0]*gc[2];
-   funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Pyramid Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Pyra5aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[4] = gc[2];
-   SHAPE_FUN_MACRO_END;
-}
-/*!
- * Init Pyramid Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Pyra5bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
-   funValue[4] = gc[2];
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Pyramid Reference coordinates ans Shape function.
- * Case A.
- */
-void GaussInfo::Pyra13aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
-
- case  5:
-   coords[0] =  0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  6:
-   coords[0] = -0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  7:
-   coords[0] = -0.5;
-   coords[1] = -0.5;
-   coords[2] =  0.0;
-   break;
- case  8:
-   coords[0] =  0.5;
-   coords[1] = -0.5;
-   coords[2] =  0.0;
-   break;
- case  9:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case 10:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case 11:
-   coords[0] = -0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case 12:
-   coords[0] =  0.0;
-   coords[1] = -0.5;
-   coords[2] =  0.5;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   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]);
-   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]);
-   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]);
-   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]);
-
-   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
-
-   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]);
-   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]);
-   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]);
-   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]);
-
-   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Quadratic Pyramid Reference coordinates ans Shape function.
- * Case B.
- */
-void GaussInfo::Pyra13bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  8:
-   coords[0] =  0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  7:
-   coords[0] = -0.5;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  6:
-   coords[0] = -0.5;
-   coords[1] = -0.5;
-   coords[2] =  0.0;
-   break;
- case  5:
-   coords[0] =  0.5;
-   coords[1] = -0.5;
-   coords[2] =  0.0;
-   break;
- case  9:
-   coords[0] =  0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case 12:
-   coords[0] =  0.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case 11:
-   coords[0] = -0.5;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case 10:
-   coords[0] =  0.0;
-   coords[1] = -0.5;
-   coords[2] =  0.5;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   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]);
-   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]);
-   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]);
-   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]);
-
-   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
-
-   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]);
-   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]);
-   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]);
-   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]);
-
-   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
-     (1.0 - gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-
-/*!
- * Init Pentahedron Reference coordinates and Shape function.
- * Case A.
- */
-void GaussInfo::Penta6aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -0.0;
-   coords[2] =  1.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
-   funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
-   funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-
-   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
-   funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
-   funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Pentahedron Reference coordinates and Shape function.
- * Case B.
- */
-void GaussInfo::Penta6bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] = -0.0;
-   coords[2] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  4:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
-   funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
-   funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
-   funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
-   funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-/*!
- * Init Pentahedron Reference coordinates and Shape function.
- * Case A.
- */
-void GaussInfo::Penta15aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] = -0.0;
-   coords[2] =  1.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  4:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-
- case  6:
-   coords[0] = -1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case  7:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  8:
-   coords[0] = -1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case  9:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 10:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case 11:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case 12:
-   coords[0] =  1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case 13:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case 14:
-   coords[0] =  1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
-   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
-   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
-
-   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
-   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
-   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
-
-   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
-   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-
-   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
-   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
-   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
-
-   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
-   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Qaudratic Pentahedron Reference coordinates and Shape function.
- * Case B.
- */
-void GaussInfo::Penta15bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  2:
-   coords[0] = -1.0;
-   coords[1] = -0.0;
-   coords[2] =  1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case  4:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
-
- case  8:
-   coords[0] = -1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case  7:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  6:
-   coords[0] = -1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
- case 12:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 14:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case 13:
-   coords[0] =  0.0;
-   coords[1] =  0.0;
-   coords[2] =  0.0;
-   break;
- case 11:
-   coords[0] =  1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.5;
-   break;
- case 10:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  0.5;
-   break;
- case  9:
-   coords[0] =  1.0;
-   coords[1] =  0.5;
-   coords[2] =  0.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
-   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
-   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
-
-   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
-   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
-   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
-
-   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
-   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
-
-   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
-   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
-   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
-
-   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
-   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Hehahedron Reference coordinates and Shape function.
- * Case A.
- */
-void GaussInfo::Hexa8aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  1:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  3:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  4:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  6:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case  7:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-
-   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Hehahedron Reference coordinates and Shape function.
- * Case B.
- */
-void GaussInfo::Hexa8bInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  4:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  7:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  6:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-
-   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Qaudratic Hehahedron Reference coordinates and Shape function.
- * Case A.
- */
-void GaussInfo::Hexa20aInit() 
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  1:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  3:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  4:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  6:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case  7:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
-
- case  8:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  9:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] = -1.0;
-   break;
- case 10:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case 11:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] = -1.0;
-   break;
- case 12:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case 13:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case 14:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 15:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 16:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case 17:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case 18:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case 19:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
-     (-2.0 - gc[0] - gc[1] - gc[2]);
-   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
-     (-2.0 + gc[0] - gc[1] - gc[2]);
-   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
-     (-2.0 + gc[0] + gc[1] - gc[2]);
-   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
-     (-2.0 - gc[0] + gc[1] - gc[2]);
-   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
-     (-2.0 - gc[0] - gc[1] + gc[2]);
-   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
-     (-2.0 + gc[0] - gc[1] + gc[2]);
-   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
-     (-2.0 + gc[0] + gc[1] + gc[2]);
-   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
-     (-2.0 - gc[0] + gc[1] + gc[2]);
-
-   funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
-   funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-   funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
-   funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
-   funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
-   funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
-   funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
-   funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
-   funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-/*!
- * Init Qaudratic Hehahedron Reference coordinates and Shape function.
- * Case B.
- */
-void GaussInfo::Hexa20bInit()
-{
-  LOCAL_COORD_MACRO_BEGIN;
- case  0:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  3:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case  2:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  1:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  4:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  7:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case  6:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case  5:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
-
- case 11:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] = -1.0;
-   break;
- case 10:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] = -1.0;
-   break;
- case  9:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] = -1.0;
-   break;
- case  8:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] = -1.0;
-   break;
- case 16:
-   coords[0] = -1.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case 19:
-   coords[0] =  1.0;
-   coords[1] = -1.0;
-   coords[2] =  0.0;
-   break;
- case 18:
-   coords[0] =  1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 17:
-   coords[0] = -1.0;
-   coords[1] =  1.0;
-   coords[2] =  0.0;
-   break;
- case 15:
-   coords[0] =  0.0;
-   coords[1] = -1.0;
-   coords[2] =  1.0;
-   break;
- case 14:
-   coords[0] =  1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
- case 13:
-   coords[0] =  0.0;
-   coords[1] =  1.0;
-   coords[2] =  1.0;
-   break;
- case 12:
-   coords[0] = -1.0;
-   coords[1] =  0.0;
-   coords[2] =  1.0;
-   break;
-   LOCAL_COORD_MACRO_END;
-
-   SHAPE_FUN_MACRO_BEGIN;
-
-   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
-     (-2.0 - gc[0] - gc[1] - gc[2]);
-   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
-     (-2.0 + gc[0] - gc[1] - gc[2]);
-   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
-     (-2.0 + gc[0] + gc[1] - gc[2]);
-   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
-     (-2.0 - gc[0] + gc[1] - gc[2]);
-   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
-     (-2.0 - gc[0] - gc[1] + gc[2]);
-   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
-     (-2.0 + gc[0] - gc[1] + gc[2]);
-   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
-     (-2.0 + gc[0] + gc[1] + gc[2]);
-   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
-     (-2.0 - gc[0] + gc[1] + gc[2]);
-
-   funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
-   funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
-   funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
-   funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
-   funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
-   funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
-   funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
-   funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
-   funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
-   funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
-   funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
-   funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
-   SHAPE_FUN_MACRO_END;
-}
-
-
-
-////////////////////////////////////////////////////////////////////////////////////////////////
-//                                GAUSS COORD CLASS                                           //
-////////////////////////////////////////////////////////////////////////////////////////////////
-/*!
- * Constructor
- */
-GaussCoords::GaussCoords()
-{
-}
-
-/*!
- * Destructor
- */
-GaussCoords::~GaussCoords()
-{
-  GaussInfoVector::iterator it = myGaussInfo.begin();
-  for( ; it != myGaussInfo.end(); it++ ) 
-    {
-      if((*it) != NULL)
-        delete (*it);
-    }
-}
-
-/*!
- * Add Gauss localization info 
- */
-void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
-                                int coordDim,
-                                const double* theGaussCoord,
-                                int theNbGauss,
-                                const double* theReferenceCoord,
-                                int theNbRef) throw (INTERP_KERNEL::Exception) 
-{
-  GaussInfoVector::iterator it = myGaussInfo.begin();
-  for( ; it != myGaussInfo.end(); it++ ) 
-    {
-      if( (*it)->GetCellType() == theGeometry ) 
-        {
-          break;
-        }
-    }
-
-  DataVector aGaussCoord;
-  for(int i = 0 ; i < theNbGauss*coordDim; i++ )
-    aGaussCoord.push_back(theGaussCoord[i]);
-
-  DataVector aReferenceCoord;
-  for(int i = 0 ; i < theNbRef*coordDim; i++ )
-    aReferenceCoord.push_back(theReferenceCoord[i]);
-
-
-  GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
-  info->InitLocalInfo();
-
-  //If info with cell type doesn't exist add it
-  if( it == myGaussInfo.end() ) 
-    {
-      myGaussInfo.push_back(info);
-
-      // If information exists, update it
-    } else 
-    {
-      int index = std::distance(myGaussInfo.begin(),it);
-      delete (*it);
-      myGaussInfo[index] = info;
-    }
-}
-
-
-/*!
- * Calculate gauss points coordinates
- */
-double* GaussCoords::CalculateCoords( NormalizedCellType theGeometry, 
-                                      const double* theNodeCoords, 
-                                      const int theSpaceDim,
-                                      const int* theIndex) 
-  throw (INTERP_KERNEL::Exception) 
-{
-  GaussInfoVector::const_iterator it = myGaussInfo.begin();
-  GaussInfoVector::const_iterator it_end = myGaussInfo.end();
-
-  //Try to find gauss localization info
-  for( ; it != it_end ; it++ ) 
-    {
-      if( (*it)->GetCellType() == theGeometry ) 
-        {
-          break;
-        }
-    }
-  if(it == it_end) 
-    {
-      throw INTERP_KERNEL::Exception("Can't find gauss localization information !!!");
-    }
-  const GaussInfo* info = (*it);
-  int aConn = info->GetNbRef();
-
-  int nbCoords = theSpaceDim * info->GetNbGauss();
-  double* aCoords = new double [nbCoords];
-  for( int i = 0; i < nbCoords; i++ )
-    aCoords[i] = 0.0;
-
-// #ifdef MYDEBUG
-//   std::cout<<"#######################################################"<<std::endl;
-//   std::cout<<"Cell type : "<<info->GetCellType()<<std::endl;
-// #endif
-
-  for( int gaussId = 0; gaussId < info->GetNbGauss(); gaussId++ ) 
-    {
-// #ifdef MYDEBUG
-//       std::cout<<"Gauss ID = "<<gaussId<<std::endl;
-// #endif
-
-      double* coord = &aCoords[ gaussId*theSpaceDim ];
-      const double* function = info->GetFunctionValues(gaussId);
-      for ( int connId = 0; connId < aConn ; connId++ ) 
-        {
-          const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
-
-// #ifdef MYDEBUG
-//           std::cout<<"Function Value["<<connId<<"] = "<<function[connId]<<std::endl;
-
-//           std::cout<<"Node #"<<connId <<" ";
-//           for( int dimId = 0; dimId < theSpaceDim; dimId++ ) 
-//             {
-//               switch(dimId)
-//                 {
-//                 case 0:
-//                   std::cout<<"( "<<nodeCoord[dimId];
-//                   break;
-//                 case 1:
-//                   std::cout<<", "<<nodeCoord[dimId];
-//                   break;
-//                 case 2:
-//                   std::cout<<", "<<nodeCoord[dimId]<<" )";
-//                   break;
-//                 }
-//             }
-//           std::cout<<std::endl;
-// #endif
-
-          for( int dimId = 0; dimId < theSpaceDim; dimId++ ) 
-            {
-              coord[dimId] += nodeCoord[dimId]*function[connId];
-            }
-        }
-    }
-  return aCoords;
-}
diff --git a/src/INTERP_KERNEL/GaussPoints/GaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/GaussCoords.hxx
deleted file mode 100644
index 7bca30eab..000000000
--- a/src/INTERP_KERNEL/GaussPoints/GaussCoords.hxx
+++ /dev/null
@@ -1,149 +0,0 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D
-//
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
-//
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
-//
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
-//
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-
-#ifndef __INTERPKERNELGAUSS_HXX__
-#define __INTERPKERNELGAUSS_HXX__
-
-#include "NormalizedUnstructuredMesh.hxx"
-#include "InterpKernelException.hxx"
-#include <vector>
-
-namespace INTERP_KERNEL 
-{
-  typedef std::vector<double> DataVector;
-  typedef std::vector<int>    IndexVector;
-
-  //Class to store Gauss Points information
-  class GaussInfo 
-  {
-  public:
-    GaussInfo( NormalizedCellType theGeometry,
-               const DataVector& theGaussCoord,
-               int theNbGauss,
-               const DataVector& theReferenceCoord,
-               int theNbRef
-               );
-    ~GaussInfo();
-
-    NormalizedCellType GetCellType() const;    
-
-    int GetGaussCoordDim() const;
-    int GetReferenceCoordDim() const;
-
-    int GetNbGauss() const;
-    int GetNbRef() const;
-
-    const double* GetFunctionValues( const int theGaussId ) const;
-
-    void InitLocalInfo() throw (INTERP_KERNEL::Exception);
-
-  protected:
-
-    bool isSatisfy();
-
-    //1D
-    void Seg2Init();
-    void Seg3Init();
-
-    //2D
-    void Tria3aInit();
-    void Tria3bInit();
-    void Tria6aInit();
-    void Tria6bInit();
-
-    void Quad4aInit();
-    void Quad4bInit();
-    void Quad8aInit();
-    void Quad8bInit();
-
-    //3D
-    void Tetra4aInit();
-    void Tetra4bInit();
-    void Tetra10aInit();
-    void Tetra10bInit();
-
-    void Pyra5aInit();
-    void Pyra5bInit();
-    void Pyra13aInit();
-    void Pyra13bInit();
-
-    void Penta6aInit();
-    void Penta6bInit();
-    void Penta15aInit();
-    void Penta15bInit();
-
-    void Hexa8aInit();
-    void Hexa8bInit();
-    void Hexa20aInit();
-    void Hexa20bInit();
-
-
-  private:
-    //INFORMATION from MEDMEM
-    NormalizedCellType myGeometry;               //Cell type
-
-    int                myNbGauss;                //Nb of the gauss points for element
-    DataVector         myGaussCoord;             //Gauss coordinates
-
-    int                myNbRef;                  //Nb of the nodes for element:
-                                                 //NORM_SEG2 - 2
-                                                 //NORM_SEG3 - 3
-                                                 //NORM_TRI3 - 3
-                                                 //.............
-
-    DataVector         myReferenceCoord;         //Reference coordinates
-
-    //LOCAL INFORMATION
-    DataVector         myLocalReferenceCoord;    //Vector to store reference coordinates
-    int                myLocalRefDim;            //Dimension of the local reference coordinates:
-                                                 // (x)       - 1D case
-                                                 // (x, y)    - 2D case
-                                                 // (x, y, z) - 3D case
-    int                myLocalNbRef;             //Nb of the local reference coordinates
-
-    DataVector         myFunctionValue;          //Shape Function values
-  };
-
-
-  //Class for calculation of the coordinates of the gauss points 
-  class GaussCoords 
-  {
-  public:
-
-    GaussCoords();
-    ~GaussCoords();
-
-    void addGaussInfo( NormalizedCellType theGeometry,
-                       int coordDim,
-                       const double* theGaussCoord,
-                       int theNbGauss,
-                       const double* theReferenceCoord,
-                       int theNbRef) throw (INTERP_KERNEL::Exception);
-
-    double* CalculateCoords( NormalizedCellType theGeometry, 
-                             const double* theNodeCoords, 
-                             const int theDimSpace,
-                             const int* theIndex
-                             ) throw(INTERP_KERNEL::Exception);
-  private:
-    typedef std::vector<GaussInfo*> GaussInfoVector;
-    GaussInfoVector myGaussInfo;
-  };
-}
-#endif //INTERPKERNELGAUSS
diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
new file mode 100644
index 000000000..57d274371
--- /dev/null
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -0,0 +1,2126 @@
+//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D
+//
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License.
+//
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+//Local includes
+#include "InterpKernelGaussCoords.hxx"
+#include "CellModel.hxx"
+
+//STL includes
+#include <math.h>
+#include <algorithm>
+#include <sstream>
+
+using namespace INTERP_KERNEL;
+
+//Define common part of the code in the MACRO
+//---------------------------------------------------------------
+#define LOCAL_COORD_MACRO_BEGIN                                         \
+  _my_local_reference_coord.resize( _my_local_ref_dim*_my_local_nb_ref );           \
+  for( int refId = 0; refId < _my_local_nb_ref; refId++ )                   \
+    {                                                                   \
+      double* coords = &_my_local_reference_coord[ refId*_my_local_ref_dim ];   \
+      switch(refId)                                                     \
+        {
+
+//---------------------------------------------------------------
+#define LOCAL_COORD_MACRO_END                   \
+  }                                             \
+}
+
+//---------------------------------------------------------------
+#define SHAPE_FUN_MACRO_BEGIN                                           \
+  for( int gaussId     = 0 ; gaussId < _my_nb_gauss ; gaussId++ )          \
+    {                                                                   \
+      double* funValue =  &_my_function_value[ gaussId * _my_nb_ref ];        \
+      const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ];
+
+//---------------------------------------------------------------
+#define SHAPE_FUN_MACRO_END                     \
+  }
+
+#define CHECK_MACRO                                                        \
+  if( ! aSatify )                                                          \
+    {                                                                      \
+      std::ostringstream stream;                                           \
+      stream << "Error in the gauss localization for the cell with type "; \
+      stream << cellModel.getRepr();                                       \
+      stream << " !!!";                                                    \
+      throw INTERP_KERNEL::Exception(stream.str().c_str());                \
+    }
+
+
+//---------------------------------------------------------------
+static bool IsEqual(double theLeft, double theRight) 
+{
+  static double EPS = 1.0E-3;
+  if(fabs(theLeft) + fabs(theRight) > EPS)
+    return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
+  return true;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+//                                GAUSS INFO CLASS                                            //
+////////////////////////////////////////////////////////////////////////////////////////////////
+
+/*!
+ * Constructor of the GaussInfo
+ */
+GaussInfo::GaussInfo( NormalizedCellType theGeometry,
+                      const DataVector& theGaussCoord,
+                      int theNbGauss,
+                      const DataVector& theReferenceCoord,
+                      int theNbRef ) :
+  _my_geometry(theGeometry),
+  _my_nb_gauss(theNbGauss),
+  _my_gauss_coord(theGaussCoord),
+  _my_nb_ref(theNbRef),
+  _my_reference_coord(theReferenceCoord)
+{
+
+  //Allocate shape function values
+  _my_function_value.resize( _my_nb_gauss * _my_nb_ref );
+}
+
+/*!
+ * Destructor
+ */
+GaussInfo::~GaussInfo()
+{
+}
+
+/*!
+ * Return dimension of the gauss coordinates
+ */
+int GaussInfo::getGaussCoordDim() const 
+{
+  if( _my_nb_gauss ) 
+    {
+      return _my_gauss_coord.size()/_my_nb_gauss;
+    }
+  else 
+    {
+      return 0;
+    }
+}
+
+/*!
+ * Return dimension of the reference coordinates
+ */
+int GaussInfo::getReferenceCoordDim() const 
+{
+  if( _my_nb_ref ) 
+    {
+      return _my_reference_coord.size()/_my_nb_ref;
+    }
+  else 
+    {
+      return 0;
+    }
+}
+
+/*!
+ * Return type of the cell.
+ */
+NormalizedCellType GaussInfo::getCellType() const 
+{
+  return _my_geometry;
+}
+
+/*!
+ * Return Nb of the gauss points.
+ */
+int GaussInfo::getNbGauss() const 
+{
+  return _my_nb_gauss;
+}
+
+/*!
+ * Return Nb of the reference coordinates.
+ */
+int GaussInfo::getNbRef() const 
+{
+  return _my_nb_ref;
+}
+
+/*!
+ * Check coordinates
+ */
+bool GaussInfo::isSatisfy() 
+{
+
+  bool anIsSatisfy = ((_my_local_nb_ref == _my_nb_ref) && (_my_local_ref_dim == getReferenceCoordDim()));
+  //Check coordinates
+  if(anIsSatisfy)
+    {
+      for( int refId = 0; refId < _my_local_nb_ref; refId++ ) 
+        {
+          double* refCoord = &_my_reference_coord[ refId*_my_local_ref_dim ];
+          double* localRefCoord = &_my_local_reference_coord[ refId*_my_local_ref_dim ];
+          bool anIsEqual = false;
+          for( int dimId = 0; dimId < _my_local_ref_dim; dimId++ ) 
+            {
+              anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
+              if(!anIsEqual ) 
+                {
+                  return false;
+                }
+            }
+        }
+    }
+  return anIsSatisfy;
+}
+
+/*!
+ * Initialize the internal vectors
+ */
+void GaussInfo::initLocalInfo() throw (INTERP_KERNEL::Exception) 
+{
+  bool aSatify = false;
+  const CellModel& cellModel=CellModel::getCellModel(_my_geometry);
+  switch( _my_geometry ) 
+    {
+    case NORM_SEG2:
+      _my_local_ref_dim = 1;
+      _my_local_nb_ref  = 2;
+      seg2Init();
+      aSatify = isSatisfy();
+      CHECK_MACRO;
+      break;
+
+    case NORM_SEG3:
+      _my_local_ref_dim = 1;
+      _my_local_nb_ref  = 3;
+      seg3Init();
+      aSatify = isSatisfy();
+      CHECK_MACRO;
+      break;
+
+    case NORM_TRI3:
+      _my_local_ref_dim = 2;
+      _my_local_nb_ref  = 3;
+      tria3aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          tria3bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_TRI6:
+      _my_local_ref_dim = 2;
+      _my_local_nb_ref  = 6;
+      tria6aInit();
+      aSatify = isSatisfy();
+      if(!aSatify)
+        {
+          tria6bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_QUAD4:
+      _my_local_ref_dim = 2;
+      _my_local_nb_ref  = 4;
+      quad4aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          quad4bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_QUAD8:
+      _my_local_ref_dim = 2;
+      _my_local_nb_ref  = 8;
+      quad8aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          quad8bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_TETRA4:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 4;
+      tetra4aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          tetra4bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_TETRA10:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 10;
+      tetra10aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          tetra10bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_PYRA5:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 5;
+      pyra5aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          pyra5bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_PYRA13:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 13;
+      pyra13aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          pyra13bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_PENTA6:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 6;
+      penta6aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          penta6bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_PENTA15:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 15;
+      penta15aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          penta15bInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_HEXA8:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 8;
+      hexa8aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          hexa8aInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    case NORM_HEXA20:
+      _my_local_ref_dim = 3;
+      _my_local_nb_ref  = 20;
+      hexa20aInit();
+      aSatify = isSatisfy();
+
+      if(!aSatify)
+        {
+          hexa20aInit();
+          aSatify = isSatisfy();
+          CHECK_MACRO;
+        }
+      break;
+
+    default:
+      throw INTERP_KERNEL::Exception("Not manged cell type !");
+      break;
+    }
+}
+
+/**
+ * Return shape function value by node id
+ */
+const double* GaussInfo::getFunctionValues( const int theGaussId ) const 
+{
+  return &_my_function_value[ _my_nb_ref*theGaussId ];
+}
+
+/*!
+ * Init Segment 2 Reference coordinates ans Shape function.
+ */
+void GaussInfo::seg2Init() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+  case  0:
+    coords[0] = -1.0;
+    break;
+  case  1:
+    coords[0] =  1.0;
+    break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 - gc[0]);
+   funValue[1] = 0.5*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Segment 3 Reference coordinates ans Shape function.
+ */
+void GaussInfo::seg3Init() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+  case  0:
+    coords[0] = -1.0;
+    break;
+  case  1:
+    coords[0] =  1.0;
+    break;
+  case  2:
+    coords[0] =  0.0;
+    break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 - gc[0])*gc[0];
+   funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
+   funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Triangle Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::tria3aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 + gc[1]);
+   funValue[1] = -0.5*(gc[0] + gc[1]);
+   funValue[2] = 0.5*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Triangle Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::tria3bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 1.0 - gc[0] - gc[1];
+   funValue[1] = gc[0];
+   funValue[2] = gc[1];
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Triangle Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::tria6aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
+   funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
+   funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
+   funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
+   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
+   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Triangle Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::tria6bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   break;
+
+ case  1:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   break;
+
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   break;
+
+ case  3:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   break;
+
+ case  4:
+   coords[0] =  0.5;
+   coords[1] =  0.5;
+   break;
+
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   break;
+
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
+   funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
+   funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
+   funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
+   funValue[4] = 4.0*gc[0]*gc[1];
+   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadrangle Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::quad4aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   break;
+
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
+   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
+   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
+   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadrangle Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::quad4bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
+   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
+   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
+   SHAPE_FUN_MACRO_END;
+}
+
+
+/*!
+ * Init Quadratic Quadrangle Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::quad8aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   break;
+ case  4:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   break;
+ case  6:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
+   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
+   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
+   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
+   funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
+   funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
+   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
+   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Quadrangle Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::quad8bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   break;
+ case  7:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
+   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
+   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
+   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
+   funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
+   funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
+   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
+   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Tetrahedron Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::tetra4aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = gc[1];
+   funValue[1] = gc[2];
+   funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
+   funValue[3] = gc[0];
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Tetrahedron Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::tetra4bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = gc[1];
+   funValue[2] = gc[2];
+   funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
+   funValue[3] = gc[0];
+   SHAPE_FUN_MACRO_END;
+
+}
+
+/*!
+ * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::tetra10aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  8:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  9:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
+   funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
+   funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
+   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
+   funValue[4] = 4.0*gc[1]*gc[2];
+   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
+   funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
+   funValue[7] = 4.0*gc[0]*gc[1];
+   funValue[8] = 4.0*gc[0]*gc[2];
+   funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::tetra10bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case  5:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] =  0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  9:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  8:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
+   funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
+   funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
+   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
+   funValue[6] = 4.0*gc[1]*gc[2];
+   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
+   funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
+   funValue[7] = 4.0*gc[0]*gc[1];
+   funValue[9] = 4.0*gc[0]*gc[2];
+   funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Pyramid Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::pyra5aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[1] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[3] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[4] = gc[2];
+   SHAPE_FUN_MACRO_END;
+}
+/*!
+ * Init Pyramid Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::pyra5bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.25*(-gc[0] + gc[1] - 1.0)*(-gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[3] = 0.25*(-gc[0] - gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[2] = 0.25*(+gc[0] + gc[1] - 1.0)*(+gc[0] - gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[1] = 0.25*(+gc[0] + gc[1] - 1.0)*(-gc[0] + gc[1] - 1.0)*(1.0 - gc[2]);
+   funValue[4] = gc[2];
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Pyramid Reference coordinates ans Shape function.
+ * Case A.
+ */
+void GaussInfo::pyra13aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+
+ case  5:
+   coords[0] =  0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] = -0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] = -0.5;
+   coords[1] = -0.5;
+   coords[2] =  0.0;
+   break;
+ case  8:
+   coords[0] =  0.5;
+   coords[1] = -0.5;
+   coords[2] =  0.0;
+   break;
+ case  9:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case 10:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case 11:
+   coords[0] = -0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case 12:
+   coords[0] =  0.0;
+   coords[1] = -0.5;
+   coords[2] =  0.5;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   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]);
+   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]);
+   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]);
+   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]);
+
+   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
+
+   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]);
+   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]);
+   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]);
+   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]);
+
+   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[10] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[12] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Quadratic Pyramid Reference coordinates ans Shape function.
+ * Case B.
+ */
+void GaussInfo::pyra13bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  8:
+   coords[0] =  0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  7:
+   coords[0] = -0.5;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  6:
+   coords[0] = -0.5;
+   coords[1] = -0.5;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  0.5;
+   coords[1] = -0.5;
+   coords[2] =  0.0;
+   break;
+ case  9:
+   coords[0] =  0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case 12:
+   coords[0] =  0.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case 11:
+   coords[0] = -0.5;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case 10:
+   coords[0] =  0.0;
+   coords[1] = -0.5;
+   coords[2] =  0.5;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   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]);
+   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]);
+   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]);
+   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]);
+
+   funValue[4] = 2.0*gc[2]*(gc[2] - 0.5);
+
+   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]);
+   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]);
+   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]);
+   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]);
+
+   funValue[9] = 0.5*gc[2]*(-gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] - gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[12] = 0.5*gc[2]*(-gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] - gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[11] = 0.5*gc[2]*(gc[0] - gc[1] + gc[2] - 1.0)*(gc[0] + gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   funValue[10] = 0.5*gc[2]*(gc[0] + gc[1] + gc[2] - 1.0)*(-gc[0] + gc[1] + gc[2] - 1.0)/
+     (1.0 - gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+
+/*!
+ * Init Pentahedron Reference coordinates and Shape function.
+ * Case A.
+ */
+void GaussInfo::penta6aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -0.0;
+   coords[2] =  1.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
+   funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
+   funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+
+   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
+   funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
+   funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Pentahedron Reference coordinates and Shape function.
+ * Case B.
+ */
+void GaussInfo::penta6bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] = -0.0;
+   coords[2] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  4:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
+   funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
+   funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
+   funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
+   funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+/*!
+ * Init Pentahedron Reference coordinates and Shape function.
+ * Case A.
+ */
+void GaussInfo::penta15aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] = -0.0;
+   coords[2] =  1.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  4:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+
+ case  6:
+   coords[0] = -1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case  7:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  8:
+   coords[0] = -1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case  9:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 10:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case 11:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case 12:
+   coords[0] =  1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case 13:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case 14:
+   coords[0] =  1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
+   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
+   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
+
+   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
+   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
+   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
+
+   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
+   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+
+   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
+   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
+   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
+
+   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
+   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Qaudratic Pentahedron Reference coordinates and Shape function.
+ * Case B.
+ */
+void GaussInfo::penta15bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  2:
+   coords[0] = -1.0;
+   coords[1] = -0.0;
+   coords[2] =  1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case  4:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+
+ case  8:
+   coords[0] = -1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case  7:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  6:
+   coords[0] = -1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+ case 12:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 14:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case 13:
+   coords[0] =  0.0;
+   coords[1] =  0.0;
+   coords[2] =  0.0;
+   break;
+ case 11:
+   coords[0] =  1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.5;
+   break;
+ case 10:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  0.5;
+   break;
+ case  9:
+   coords[0] =  1.0;
+   coords[1] =  0.5;
+   coords[2] =  0.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
+   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
+   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
+
+   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
+   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
+   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
+
+   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
+   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
+
+   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
+   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
+   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
+
+   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
+   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Hehahedron Reference coordinates and Shape function.
+ * Case A.
+ */
+void GaussInfo::hexa8aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  4:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  6:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case  7:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+
+   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Hehahedron Reference coordinates and Shape function.
+ * Case B.
+ */
+void GaussInfo::hexa8bInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  4:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  7:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  6:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+
+   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Qaudratic Hehahedron Reference coordinates and Shape function.
+ * Case A.
+ */
+void GaussInfo::hexa20aInit() 
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  1:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  3:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  4:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  6:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case  7:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+
+ case  8:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  9:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] = -1.0;
+   break;
+ case 10:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case 11:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] = -1.0;
+   break;
+ case 12:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case 13:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case 14:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 15:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 16:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case 17:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case 18:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case 19:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
+     (-2.0 - gc[0] - gc[1] - gc[2]);
+   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
+     (-2.0 + gc[0] - gc[1] - gc[2]);
+   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
+     (-2.0 + gc[0] + gc[1] - gc[2]);
+   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
+     (-2.0 - gc[0] + gc[1] - gc[2]);
+   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
+     (-2.0 - gc[0] - gc[1] + gc[2]);
+   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
+     (-2.0 + gc[0] - gc[1] + gc[2]);
+   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
+     (-2.0 + gc[0] + gc[1] + gc[2]);
+   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
+     (-2.0 - gc[0] + gc[1] + gc[2]);
+
+   funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
+   funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+   funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
+   funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
+   funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
+   funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
+   funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
+   funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+/*!
+ * Init Qaudratic Hehahedron Reference coordinates and Shape function.
+ * Case B.
+ */
+void GaussInfo::hexa20bInit()
+{
+  LOCAL_COORD_MACRO_BEGIN;
+ case  0:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  3:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case  2:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  1:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  4:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  7:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case  6:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case  5:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+
+ case 11:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] = -1.0;
+   break;
+ case 10:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] = -1.0;
+   break;
+ case  9:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] = -1.0;
+   break;
+ case  8:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] = -1.0;
+   break;
+ case 16:
+   coords[0] = -1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case 19:
+   coords[0] =  1.0;
+   coords[1] = -1.0;
+   coords[2] =  0.0;
+   break;
+ case 18:
+   coords[0] =  1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 17:
+   coords[0] = -1.0;
+   coords[1] =  1.0;
+   coords[2] =  0.0;
+   break;
+ case 15:
+   coords[0] =  0.0;
+   coords[1] = -1.0;
+   coords[2] =  1.0;
+   break;
+ case 14:
+   coords[0] =  1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+ case 13:
+   coords[0] =  0.0;
+   coords[1] =  1.0;
+   coords[2] =  1.0;
+   break;
+ case 12:
+   coords[0] = -1.0;
+   coords[1] =  0.0;
+   coords[2] =  1.0;
+   break;
+   LOCAL_COORD_MACRO_END;
+
+   SHAPE_FUN_MACRO_BEGIN;
+
+   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
+     (-2.0 - gc[0] - gc[1] - gc[2]);
+   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
+     (-2.0 + gc[0] - gc[1] - gc[2]);
+   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
+     (-2.0 + gc[0] + gc[1] - gc[2]);
+   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
+     (-2.0 - gc[0] + gc[1] - gc[2]);
+   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
+     (-2.0 - gc[0] - gc[1] + gc[2]);
+   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
+     (-2.0 + gc[0] - gc[1] + gc[2]);
+   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
+     (-2.0 + gc[0] + gc[1] + gc[2]);
+   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
+     (-2.0 - gc[0] + gc[1] + gc[2]);
+
+   funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
+   funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
+   funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
+   funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
+   funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
+   funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
+   funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
+   funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
+   funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
+   funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
+   funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
+   funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
+   SHAPE_FUN_MACRO_END;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+//                                GAUSS COORD CLASS                                           //
+////////////////////////////////////////////////////////////////////////////////////////////////
+/*!
+ * Constructor
+ */
+GaussCoords::GaussCoords()
+{
+}
+
+/*!
+ * Destructor
+ */
+GaussCoords::~GaussCoords()
+{
+  GaussInfoVector::iterator it = _my_gauss_info.begin();
+  for( ; it != _my_gauss_info.end(); it++ ) 
+    {
+      if((*it) != NULL)
+        delete (*it);
+    }
+}
+
+/*!
+ * Add Gauss localization info 
+ */
+void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
+                                int coordDim,
+                                const double* theGaussCoord,
+                                int theNbGauss,
+                                const double* theReferenceCoord,
+                                int theNbRef) throw (INTERP_KERNEL::Exception) 
+{
+  GaussInfoVector::iterator it = _my_gauss_info.begin();
+  for( ; it != _my_gauss_info.end(); it++ ) 
+    {
+      if( (*it)->getCellType() == theGeometry ) 
+        {
+          break;
+        }
+    }
+
+  DataVector aGaussCoord;
+  for(int i = 0 ; i < theNbGauss*coordDim; i++ )
+    aGaussCoord.push_back(theGaussCoord[i]);
+
+  DataVector aReferenceCoord;
+  for(int i = 0 ; i < theNbRef*coordDim; i++ )
+    aReferenceCoord.push_back(theReferenceCoord[i]);
+
+
+  GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, theNbGauss, aReferenceCoord, theNbRef);
+  info->initLocalInfo();
+
+  //If info with cell type doesn't exist add it
+  if( it == _my_gauss_info.end() ) 
+    {
+      _my_gauss_info.push_back(info);
+
+      // If information exists, update it
+    }
+  else 
+    {
+      int index = std::distance(_my_gauss_info.begin(),it);
+      delete (*it);
+      _my_gauss_info[index] = info;
+    }
+}
+
+
+/*!
+ * Calculate gauss points coordinates
+ */
+double* GaussCoords::calculateCoords( NormalizedCellType theGeometry, 
+                                      const double* theNodeCoords, 
+                                      const int theSpaceDim,
+                                      const int* theIndex) 
+  throw (INTERP_KERNEL::Exception) 
+{
+  GaussInfoVector::const_iterator it = _my_gauss_info.begin();
+  GaussInfoVector::const_iterator it_end = _my_gauss_info.end();
+
+  //Try to find gauss localization info
+  for( ; it != it_end ; it++ ) 
+    {
+      if( (*it)->getCellType() == theGeometry ) 
+        {
+          break;
+        }
+    }
+  if(it == it_end) 
+    {
+      throw INTERP_KERNEL::Exception("Can't find gauss localization information !!!");
+    }
+  const GaussInfo* info = (*it);
+  int aConn = info->getNbRef();
+
+  int nbCoords = theSpaceDim * info->getNbGauss();
+  double* aCoords = new double [nbCoords];
+  for( int i = 0; i < nbCoords; i++ )
+    aCoords[i] = 0.0;
+
+  for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ ) 
+    {
+      double* coord = &aCoords[ gaussId*theSpaceDim ];
+      const double* function = info->getFunctionValues(gaussId);
+      for ( int connId = 0; connId < aConn ; connId++ ) 
+        {
+          const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
+          for( int dimId = 0; dimId < theSpaceDim; dimId++ ) 
+            {
+              coord[dimId] += nodeCoord[dimId]*function[connId];
+            }
+        }
+    }
+  return aCoords;
+}
diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
new file mode 100644
index 000000000..7405cedc6
--- /dev/null
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
@@ -0,0 +1,150 @@
+//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D
+//
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License.
+//
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#ifndef __INTERPKERNELGAUSS_HXX__
+#define __INTERPKERNELGAUSS_HXX__
+
+#include "NormalizedUnstructuredMesh.hxx"
+#include "InterpKernelException.hxx"
+
+#include <vector>
+
+namespace INTERP_KERNEL 
+{
+  typedef std::vector<double> DataVector;
+  typedef std::vector<int>    IndexVector;
+
+  //Class to store Gauss Points information
+  class GaussInfo 
+  {
+  public:
+    GaussInfo( NormalizedCellType theGeometry,
+               const DataVector& theGaussCoord,
+               int theNbGauss,
+               const DataVector& theReferenceCoord,
+               int theNbRef
+               );
+    ~GaussInfo();
+
+    NormalizedCellType getCellType() const;    
+
+    int getGaussCoordDim() const;
+    int getReferenceCoordDim() const;
+
+    int getNbGauss() const;
+    int getNbRef() const;
+
+    const double* getFunctionValues( const int theGaussId ) const;
+
+    void initLocalInfo() throw (INTERP_KERNEL::Exception);
+
+  protected:
+
+    bool isSatisfy();
+
+    //1D
+    void seg2Init();
+    void seg3Init();
+
+    //2D
+    void tria3aInit();
+    void tria3bInit();
+    void tria6aInit();
+    void tria6bInit();
+
+    void quad4aInit();
+    void quad4bInit();
+    void quad8aInit();
+    void quad8bInit();
+
+    //3D
+    void tetra4aInit();
+    void tetra4bInit();
+    void tetra10aInit();
+    void tetra10bInit();
+
+    void pyra5aInit();
+    void pyra5bInit();
+    void pyra13aInit();
+    void pyra13bInit();
+
+    void penta6aInit();
+    void penta6bInit();
+    void penta15aInit();
+    void penta15bInit();
+
+    void hexa8aInit();
+    void hexa8bInit();
+    void hexa20aInit();
+    void hexa20bInit();
+
+
+  private:
+    //INFORMATION from MEDMEM
+    NormalizedCellType _my_geometry;               //Cell type
+
+    int                _my_nb_gauss;                //Nb of the gauss points for element
+    DataVector         _my_gauss_coord;             //Gauss coordinates
+
+    int                _my_nb_ref;                  //Nb of the nodes for element:
+                                                 //NORM_SEG2 - 2
+                                                 //NORM_SEG3 - 3
+                                                 //NORM_TRI3 - 3
+                                                 //.............
+
+    DataVector         _my_reference_coord;         //Reference coordinates
+
+    //LOCAL INFORMATION
+    DataVector         _my_local_reference_coord;    //Vector to store reference coordinates
+    int                _my_local_ref_dim;            //Dimension of the local reference coordinates:
+                                                 // (x)       - 1D case
+                                                 // (x, y)    - 2D case
+                                                 // (x, y, z) - 3D case
+    int                _my_local_nb_ref;             //Nb of the local reference coordinates
+
+    DataVector         _my_function_value;          //Shape Function values
+  };
+
+
+  //Class for calculation of the coordinates of the gauss points 
+  class GaussCoords 
+  {
+  public:
+
+    GaussCoords();
+    ~GaussCoords();
+
+    void addGaussInfo( NormalizedCellType theGeometry,
+                       int coordDim,
+                       const double* theGaussCoord,
+                       int theNbGauss,
+                       const double* theReferenceCoord,
+                       int theNbRef) throw (INTERP_KERNEL::Exception);
+
+    double* calculateCoords( NormalizedCellType theGeometry, 
+                             const double* theNodeCoords, 
+                             const int theDimSpace,
+                             const int* theIndex
+                             ) throw(INTERP_KERNEL::Exception);
+  private:
+    typedef std::vector<GaussInfo*> GaussInfoVector;
+    GaussInfoVector _my_gauss_info;
+  };
+}
+#endif //INTERPKERNELGAUSS
diff --git a/src/INTERP_KERNEL/GaussPoints/Makefile.am b/src/INTERP_KERNEL/GaussPoints/Makefile.am
index 0a95f62c3..c73afb735 100644
--- a/src/INTERP_KERNEL/GaussPoints/Makefile.am
+++ b/src/INTERP_KERNEL/GaussPoints/Makefile.am
@@ -26,11 +26,11 @@ include $(top_srcdir)/adm_local/unix/make_common_starter.am
 noinst_LTLIBRARIES = libinterpkernelgauss.la
 
 salomeinclude_HEADERS = \
-GaussCoords.hxx
+InterpKernelGaussCoords.hxx
 
 # Libraries targets
 dist_libinterpkernelgauss_la_SOURCES = \
-GaussCoords.cxx
+InterpKernelGaussCoords.cxx
 
 libinterpkernelgauss_la_CPPFLAGS=-I$(srcdir)/../Bases -I$(srcdir)/..