Salome HOME
fix: replace unordered_set/map with set/map
[tools/medcoupling.git] / src / INTERP_KERNEL / GaussPoints / InterpKernelGaussCoords.cxx
1 // Copyright (C) 2007-2024  CEA, EDF
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 //Local includes
21 #include "InterpKernelGaussCoords.hxx"
22 #include "CellModel.hxx"
23
24 //STL includes
25 #include <math.h>
26 #include <algorithm>
27 #include <sstream>
28 #include <cmath>
29
30 using namespace INTERP_KERNEL;
31
32 const double GaussInfo::SEG2A_REF[2]={-1.0, 1.0};
33
34 const double GaussInfo::SEG2B_REF[2]={0., 1.0};
35
36 const double GaussInfo::SEG3_REF[3]={-1.0, 1.0, 0.0};
37
38 const double GaussInfo::SEG4_REF[4]={-1.0, 1.0, -0.3333333333333333, 0.3333333333333333};
39
40 const double GaussInfo::TRIA3A_REF[6]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0};
41
42 const double GaussInfo::TRIA3B_REF[6]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
43
44 const double GaussInfo::TRIA6A_REF[12]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 0.0};
45
46 const double GaussInfo::TRIA6B_REF[12]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5};
47
48 const double GaussInfo::TRIA7A_REF[14]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.3333333333333333, 0.3333333333333333};
49
50 const double GaussInfo::QUAD4A_REF[8]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
51
52 const double GaussInfo::QUAD4B_REF[8]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
53
54 const double GaussInfo::QUAD8A_REF[16]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0};
55
56 const double GaussInfo::QUAD8B_REF[16]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0};
57
58 const double GaussInfo::QUAD9A_REF[18]={-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0};
59
60 const double GaussInfo::TETRA4A_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0};
61
62 const double GaussInfo::TETRA4B_REF[12]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0};
63
64 const double GaussInfo::TETRA10A_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.0};
65
66 const double GaussInfo::TETRA10B_REF[30]={0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5};
67
68 const double GaussInfo::PYRA5A_REF[15]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0};
69
70 const double GaussInfo::PYRA5B_REF[15]={1.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
71
72 const double GaussInfo::PYRA13A_REF[39]={1.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, -0.5, 0.5, 0.0, -0.5, -0.5, 0.0, 0.5, -0.5, 0.0, 0.5, 0.0, 0.5, 0.0, 0.5, 0.5, -0.5, 0.0, 0.5, 0.0, -0.5, 0.5};
73
74 const double GaussInfo::PYRA13B_REF[39]={1.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, -0.5, 0.0, -0.5, -0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5, 0.0, 0.5, 0.0, -0.5, 0.5, -0.5, 0.0, 0.5, 0.0, 0.5, 0.5};
75
76 const double GaussInfo::PENTA6A_REF[18]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0};
77
78 const double GaussInfo::PENTA6B_REF[18]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0};
79
80 const double GaussInfo::PENTA15A_REF[45]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, 0.5, 0.5, -1.0, 0.0, 0.5, -1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0};
81
82 const double GaussInfo::PENTA15B_REF[45]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -1.0, 0.5, 0.0, -1.0, 0.0, 0.5, -1.0, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0};
83
84 const double GaussInfo::PENTA18A_REF[54]={-1.0, 1.0, 0.0, -1.0, -0.0, 1.0, -1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, 0.5, 0.5, -1.0, 0.0, 0.5, -1.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0};
85
86 const double GaussInfo::PENTA18B_REF[54]={-1.0, 1.0, 0.0, -1.0, 0.0, 0.0, -1.0, -0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, -1.0, 0.5, 0.0, -1.0, 0.0, 0.5, -1.0, 0.5, 0.5, 1.0, 0.5, 0.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.5, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5};
87
88 const double GaussInfo::HEXA8A_REF[24]={-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
89
90 const double GaussInfo::HEXA8B_REF[24]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0};
91
92 const double GaussInfo::HEXA20A_REF[60]={-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 0.0, -1.0, -1.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, -1.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, -1.0, 0.0, 1.0};
93
94 const double GaussInfo::HEXA20B_REF[60]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 0.0, -1.0, 0.0, 1.0, -1.0, 1.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, -1.0, 1.0, -1.0, -1.0, 0.0, -1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, -1.0, 0.0};
95
96 const double GaussInfo::HEXA27A_REF[81]={-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 0.0, -1.0, 0.0, 1.0, -1.0, 1.0, 0.0, -1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 0.0, -1.0, 1.0, -1.0, -1.0, 0.0, -1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0};
97
98 //Define common part of the code in the MACRO
99 //---------------------------------------------------------------
100 #define LOCAL_COORD_MACRO_BEGIN                                         \
101   _my_local_reference_coord.resize( _my_local_ref_dim*_my_local_nb_ref );           \
102   for( int refId = 0; refId < _my_local_nb_ref; refId++ )                   \
103     {                                                                   \
104       double* coords = &_my_local_reference_coord[ refId*_my_local_ref_dim ];   \
105       switch(refId)                                                     \
106         {
107
108 //---------------------------------------------------------------
109 #define LOCAL_COORD_MACRO_END                   \
110   }                                             \
111 }
112
113 //---------------------------------------------------------------
114 #define SHAPE_FUN_MACRO_BEGIN                                              \
115   for( int gaussId     = 0 ; gaussId < _my_nb_gauss ; gaussId++ )          \
116     {                                                                      \
117       double* funValue =  &_my_function_value[ gaussId * _my_nb_ref ];     \
118       const double* gc = &_my_gauss_coord[ gaussId * getGaussCoordDim() ];
119
120 //---------------------------------------------------------------
121 #define SHAPE_FUN_MACRO_END                                                \
122   }
123
124 #define DEV_SHAPE_FUN_MACRO_BEGIN                                                                               \
125   for( int gaussId = 0 ; gaussId < _my_nb_gauss ; gaussId++ )                                                   \
126     {                                                                                                           \
127       double *devFunValue = _my_derivative_func_value.data() + gaussId * getReferenceCoordDim() * _my_nb_ref;   \
128       const double *gc = _my_gauss_coord.data() + gaussId * getGaussCoordDim();
129
130 #define DEV_SHAPE_FUN_MACRO_END                                                                            \
131     }
132
133 #define CHECK_MACRO                                                        \
134   if( ! aSatify )                                                          \
135     {                                                                      \
136       std::ostringstream stream;                                           \
137       stream << "Error in the gauss localization for the cell with type "; \
138       stream << cellModel.getRepr();                                       \
139       stream << " !!!";                                                    \
140       throw INTERP_KERNEL::Exception(stream.str().c_str());                \
141     }
142
143
144 //---------------------------------------------------------------
145 static bool IsEqual(double theLeft, double theRight) 
146 {
147   static double EPS = 1.0E-3;
148   if(fabs(theLeft) + fabs(theRight) > EPS)
149     return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
150   return true;
151 }
152
153
154 ////////////////////////////////////////////////////////////////////////////////////////////////
155 //                                GAUSS INFO CLASS                                            //
156 ////////////////////////////////////////////////////////////////////////////////////////////////
157
158 /*!
159  * Constructor of the GaussInfo
160  */
161 GaussInfo::GaussInfo( NormalizedCellType theGeometry,
162                       const DataVector& theGaussCoord,
163                       int theNbGauss,
164                       const DataVector& theReferenceCoord,
165                       int theNbRef ) :
166   _my_geometry(theGeometry),
167   _my_nb_gauss(theNbGauss),
168   _my_gauss_coord(theGaussCoord),
169   _my_nb_ref(theNbRef),
170   _my_reference_coord(theReferenceCoord)
171 {
172   //Allocate shape function values
173   _my_function_value.resize( _my_nb_gauss * _my_nb_ref );
174   _my_derivative_func_value.resize( _my_nb_gauss * _my_nb_ref * getReferenceCoordDim() );
175 }
176
177 /*!
178  * Destructor
179  */
180 GaussInfo::~GaussInfo()
181 {
182 }
183
184 /*!
185  * Return dimension of the gauss coordinates
186  */
187 int GaussInfo::getGaussCoordDim() const 
188 {
189   if( _my_nb_gauss ) 
190     {
191       return (int)_my_gauss_coord.size()/_my_nb_gauss;
192     }
193   else 
194     {
195       return 0;
196     }
197 }
198
199 /*!
200  * Return dimension of the reference coordinates
201  */
202 int GaussInfo::getReferenceCoordDim() const 
203 {
204   if( _my_nb_ref ) 
205     {
206       return (int)(_my_reference_coord.size()/_my_nb_ref);
207     }
208   else 
209     {
210       return 0;
211     }
212 }
213
214 /*!
215  * Return type of the cell.
216  */
217 NormalizedCellType GaussInfo::getCellType() const 
218 {
219   return _my_geometry;
220 }
221
222 /*!
223  * Return Nb of the gauss points.
224  */
225 int GaussInfo::getNbGauss() const 
226 {
227   return _my_nb_gauss;
228 }
229
230 /*!
231  * Return Nb of the reference coordinates.
232  */
233 int GaussInfo::getNbRef() const 
234 {
235   return _my_nb_ref;
236 }
237
238 GaussInfo GaussInfo::convertToLinear() const
239 {
240   switch(_my_geometry)
241     {
242     case NORM_SEG3:
243       {
244         std::vector<double> a(SEG3_REF,SEG3_REF+3);
245         if(IsSatisfy(a,_my_reference_coord))
246           {
247             std::vector<double> c(SEG2A_REF,SEG2A_REF+2);
248             return GaussInfo(NORM_SEG2,_my_gauss_coord,getNbGauss(),c,2);
249           }
250         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for SEG3 !");
251       }
252     case NORM_TRI6:
253       {
254         std::vector<double> a(TRIA6A_REF,TRIA6A_REF+12),b(TRIA6B_REF,TRIA6B_REF+12);
255         if(IsSatisfy(a,_my_reference_coord))
256           {
257             std::vector<double> c(TRIA3A_REF,TRIA3A_REF+6);
258             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
259           }
260         if(IsSatisfy(b,_my_reference_coord))
261           {
262             std::vector<double> c(TRIA3B_REF,TRIA3B_REF+6);
263             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
264           }
265         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TRI6 !");
266       }
267     case NORM_TRI7:
268       {
269         std::vector<double> a(TRIA7A_REF,TRIA7A_REF+14);
270         if(IsSatisfy(a,_my_reference_coord))
271           {
272             std::vector<double> c(TRIA3B_REF,TRIA3B_REF+6);
273             return GaussInfo(NORM_TRI3,_my_gauss_coord,getNbGauss(),c,3);
274           }
275         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TRI7 !");
276       }
277     case NORM_QUAD8:
278       {
279         std::vector<double> a(QUAD8A_REF,QUAD8A_REF+16),b(QUAD8B_REF,QUAD8B_REF+16);
280         if(IsSatisfy(a,_my_reference_coord))
281           {
282             std::vector<double> c(QUAD4A_REF,QUAD4A_REF+8);
283             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
284           }
285         if(IsSatisfy(b,_my_reference_coord))
286           {
287             std::vector<double> c(QUAD4B_REF,QUAD4B_REF+8);
288             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
289           }
290         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for QUAD8 !");
291       }
292     case NORM_QUAD9:
293       {
294         std::vector<double> a(QUAD9A_REF,QUAD9A_REF+18);
295         if(IsSatisfy(a,_my_reference_coord))
296           {
297             std::vector<double> c(QUAD4B_REF,QUAD4B_REF+8);
298             return GaussInfo(NORM_QUAD4,_my_gauss_coord,getNbGauss(),c,4);
299           }
300         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for QUAD9 !");
301       }
302     case NORM_TETRA10:
303       {
304         std::vector<double> a(TETRA10A_REF,TETRA10A_REF+30),b(TETRA10B_REF,TETRA10B_REF+30);
305         if(IsSatisfy(a,_my_reference_coord))
306           {
307             std::vector<double> c(TETRA4A_REF,TETRA4A_REF+12);
308             return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
309           }
310         if(IsSatisfy(b,_my_reference_coord))
311           {
312             std::vector<double> c(TETRA4B_REF,TETRA4B_REF+12);
313             return GaussInfo(NORM_TETRA4,_my_gauss_coord,getNbGauss(),c,4);
314           }
315         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for TETRA10 !");
316       }
317     case NORM_PYRA13:
318       {
319         std::vector<double> a(PYRA13A_REF,PYRA13A_REF+39),b(PYRA13B_REF,PYRA13B_REF+39);
320         if(IsSatisfy(a,_my_reference_coord))
321           {
322             std::vector<double> c(PYRA5A_REF,PYRA5A_REF+15);
323             return GaussInfo(NORM_PYRA5,_my_gauss_coord,getNbGauss(),c,5);
324           }
325         if(IsSatisfy(b,_my_reference_coord))
326           {
327             std::vector<double> c(PYRA5B_REF,PYRA5B_REF+15);
328             return GaussInfo(NORM_PYRA5,_my_gauss_coord,getNbGauss(),c,5);
329           }
330         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PYRA13 !");
331       }
332     case NORM_PENTA15:
333       {
334         std::vector<double> a(PENTA15A_REF,PENTA15A_REF+45),b(PENTA15B_REF,PENTA15B_REF+45);
335         if(IsSatisfy(a,_my_reference_coord))
336           {
337             std::vector<double> c(PENTA6A_REF,PENTA6A_REF+18);
338             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
339           }
340         if(IsSatisfy(b,_my_reference_coord))
341           {
342             std::vector<double> c(PENTA6B_REF,PENTA6B_REF+18);
343             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
344           }
345         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PENTA15 !");
346       }
347     case NORM_PENTA18:
348       {
349         std::vector<double> a(PENTA18A_REF,PENTA18A_REF+54),b(PENTA18B_REF,PENTA18B_REF+54);
350         if(IsSatisfy(a,_my_reference_coord))
351           {
352             std::vector<double> c(PENTA6A_REF,PENTA6A_REF+18);
353             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
354           }
355         if(IsSatisfy(b,_my_reference_coord))
356           {
357             std::vector<double> c(PENTA6B_REF,PENTA6B_REF+18);
358             return GaussInfo(NORM_PENTA6,_my_gauss_coord,getNbGauss(),c,6);
359           }
360         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for PENTA18 !");
361       }
362     case NORM_HEXA20:
363       {
364         std::vector<double> a(HEXA20A_REF,HEXA20A_REF+60),b(HEXA20B_REF,HEXA20B_REF+60);
365         if(IsSatisfy(a,_my_reference_coord))
366           {
367             std::vector<double> c(HEXA8A_REF,HEXA8A_REF+24);
368             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
369           }
370         if(IsSatisfy(b,_my_reference_coord))
371           {
372             std::vector<double> c(HEXA8B_REF,HEXA8B_REF+24);
373             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
374           }
375         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for HEXA20 !");
376       }
377     case NORM_HEXA27:
378       {
379         std::vector<double> a(HEXA27A_REF,HEXA27A_REF+81);
380         if(IsSatisfy(a,_my_reference_coord))
381           {
382             std::vector<double> c(HEXA8B_REF,HEXA8B_REF+24);
383             return GaussInfo(NORM_HEXA8,_my_gauss_coord,getNbGauss(),c,8);
384           }
385         throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not recognized pattern for HEXA27 !");
386       }
387     default:
388       throw INTERP_KERNEL::Exception("GaussInfo::convertToLinear : not implemented yet for other types than TETRA10, HEXA20, HEXA27, TRI3, QUAD8, QUAD8, PYRA13, PENTA15 !");
389     }
390 }
391
392
393 bool GaussInfo::IsSatisfy(const std::vector<double>& ref1, const std::vector<double>& ref2)
394 {
395   std::size_t sz(ref1.size());
396   if(sz!=ref2.size())
397     return false;
398   for(std::size_t i=0;i<sz;i++)
399     if(!IsEqual(ref1[i],ref2[i]))
400       return false;
401   return true;
402 }
403
404 /*!
405  * Check coordinates
406  */
407 bool GaussInfo::isSatisfy() 
408 {
409
410   bool anIsSatisfy = ((_my_local_nb_ref == _my_nb_ref) && (_my_local_ref_dim == getReferenceCoordDim()));
411   //Check coordinates
412   if(anIsSatisfy)
413     {
414       for( int refId = 0; refId < _my_local_nb_ref; refId++ ) 
415         {
416           double* refCoord = &_my_reference_coord[ refId*_my_local_ref_dim ];
417           double* localRefCoord = &_my_local_reference_coord[ refId*_my_local_ref_dim ];
418           bool anIsEqual = false;
419           for( int dimId = 0; dimId < _my_local_ref_dim; dimId++ ) 
420             {
421               anIsEqual = IsEqual( localRefCoord[dimId], refCoord[dimId]);
422               if(!anIsEqual ) 
423                 {
424                   return false;
425                 }
426             }
427         }
428     }
429   return anIsSatisfy;
430 }
431
432 std::vector<double> GaussInfo::NormalizeCoordinatesIfNecessary(NormalizedCellType ct, int inputDim, const std::vector<double>& inputArray)
433 {
434   std::size_t sz(inputArray.size()),dim((std::size_t)inputDim);
435   if(dim==0)
436     throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid dimension ! Must be !=0 !");
437   if(sz%dim!=0)
438     throw INTERP_KERNEL::Exception("GaussInfo::NormalizeCoordinatesIfNecessary : invalid input array ! Inconsistent with the given dimension !");
439   const CellModel& cm(CellModel::GetCellModel(ct));
440   std::size_t baseDim((std::size_t)cm.getDimension());
441   if(baseDim==dim)
442     return inputArray;
443   std::size_t nbOfItems(sz/dim);
444   std::vector<double> ret(nbOfItems*baseDim);
445   if(baseDim>dim)
446     {
447       for(std::size_t i=0;i<nbOfItems;i++)
448         {
449           std::size_t j=0;
450           for(;j<dim;j++)
451             ret[i*baseDim+j]=inputArray[i*dim+j];
452           for(;j<baseDim;j++)
453             ret[i*baseDim+j]=0.;
454         }
455     }
456   else
457     {
458       for(std::size_t i=0;i<nbOfItems;i++)
459         {
460           std::size_t j=0;
461           for(;j<baseDim;j++)
462             ret[i*baseDim+j]=inputArray[i*dim+j];
463         }
464     }
465   return ret;
466 }
467
468 std::vector<double> GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellType ct)
469 {
470   switch(ct)
471   {
472     case INTERP_KERNEL::NORM_SEG2:
473       return std::vector<double>(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double));
474     case INTERP_KERNEL::NORM_SEG3:
475       return std::vector<double>(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double));
476     case INTERP_KERNEL::NORM_SEG4:
477       return std::vector<double>(SEG4_REF,SEG4_REF+sizeof(SEG4_REF)/sizeof(double));
478     case INTERP_KERNEL::NORM_TRI3:
479       return std::vector<double>(TRIA3A_REF,TRIA3A_REF+sizeof(TRIA3A_REF)/sizeof(double));
480     case INTERP_KERNEL::NORM_TRI6:
481       return std::vector<double>(TRIA6A_REF,TRIA6A_REF+sizeof(TRIA6A_REF)/sizeof(double));
482     case INTERP_KERNEL::NORM_TRI7:
483       return std::vector<double>(TRIA7A_REF,TRIA7A_REF+sizeof(TRIA7A_REF)/sizeof(double));
484     case INTERP_KERNEL::NORM_QUAD4:
485       return std::vector<double>(QUAD4A_REF,QUAD4A_REF+sizeof(QUAD4A_REF)/sizeof(double));
486     case INTERP_KERNEL::NORM_QUAD8:
487       return std::vector<double>(QUAD8A_REF,QUAD8A_REF+sizeof(QUAD8A_REF)/sizeof(double));
488     case INTERP_KERNEL::NORM_QUAD9:
489       return std::vector<double>(QUAD9A_REF,QUAD9A_REF+sizeof(QUAD9A_REF)/sizeof(double));
490     case INTERP_KERNEL::NORM_TETRA4:
491       return std::vector<double>(TETRA4A_REF,TETRA4A_REF+sizeof(TETRA4A_REF)/sizeof(double));
492     case INTERP_KERNEL::NORM_TETRA10:
493       return std::vector<double>(TETRA10A_REF,TETRA10A_REF+sizeof(TETRA10A_REF)/sizeof(double));
494     case INTERP_KERNEL::NORM_PYRA5:
495       return std::vector<double>(PYRA5A_REF,PYRA5A_REF+sizeof(PYRA5A_REF)/sizeof(double));
496     case INTERP_KERNEL::NORM_PYRA13:
497       return std::vector<double>(PYRA13A_REF,PYRA13A_REF+sizeof(PYRA13A_REF)/sizeof(double));
498     case INTERP_KERNEL::NORM_PENTA6:
499       return std::vector<double>(PENTA6A_REF,PENTA6A_REF+sizeof(PENTA6A_REF)/sizeof(double));
500     case INTERP_KERNEL::NORM_PENTA15:
501       return std::vector<double>(PENTA15A_REF,PENTA15A_REF+sizeof(PENTA15A_REF)/sizeof(double));
502     case INTERP_KERNEL::NORM_PENTA18:
503       return std::vector<double>(PENTA18A_REF,PENTA18A_REF+sizeof(PENTA18A_REF)/sizeof(double));
504     case INTERP_KERNEL::NORM_HEXA8:
505       return std::vector<double>(HEXA8A_REF,HEXA8A_REF+sizeof(HEXA8A_REF)/sizeof(double));
506     case INTERP_KERNEL::NORM_HEXA20:
507       return std::vector<double>(HEXA20A_REF,HEXA20A_REF+sizeof(HEXA20A_REF)/sizeof(double));
508     case INTERP_KERNEL::NORM_HEXA27:
509       return std::vector<double>(HEXA27A_REF,HEXA27A_REF+sizeof(HEXA27A_REF)/sizeof(double));
510     default:
511       THROW_IK_EXCEPTION("Input type " << ct << "is not managed by GetDefaultReferenceCoordinatesOf")
512   }
513 }
514
515 /*!
516  * Returns true if \a ptInRefCoo is in reference cell of type \a ct (regarding GetDefaultReferenceCoordinatesOf)
517  * \sa GetDefaultReferenceCoordinatesOf
518  */
519 bool GaussInfo::IsInOrOutForReference(NormalizedCellType ct, const double *ptInRefCoo, double eps)
520 {
521   switch(ct)
522   {
523     case INTERP_KERNEL::NORM_SEG2:
524     case INTERP_KERNEL::NORM_SEG3:
525     {
526       return std::fabs(ptInRefCoo[0]) < 1.0+eps;
527     }
528     case INTERP_KERNEL::NORM_QUAD4:
529     case INTERP_KERNEL::NORM_QUAD8:
530     case INTERP_KERNEL::NORM_QUAD9:
531     {
532       return std::find_if(ptInRefCoo,ptInRefCoo+2,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+2;
533     }
534     case INTERP_KERNEL::NORM_HEXA8:
535     case INTERP_KERNEL::NORM_HEXA20:
536     case INTERP_KERNEL::NORM_HEXA27:
537     {
538       return std::find_if(ptInRefCoo,ptInRefCoo+3,[eps](double v){ return std::fabs(v) > 1.0+eps; }) == ptInRefCoo+3;
539     }
540     default:
541       THROW_IK_EXCEPTION("IsInOrOutForReference : not implemented for this geo type !");
542   }
543 }
544
545 typedef void (*MapToShapeFunction)(GaussInfo& obj);
546
547 /*!
548  * Initialize the internal vectors
549  */
550 void GaussInfo::initLocalInfo()
551 {
552   bool aSatify = false;
553   const CellModel& cellModel(CellModel::GetCellModel(_my_geometry));
554   switch( _my_geometry ) 
555     {
556     case NORM_POINT1:
557       _my_local_ref_dim = 0;
558       _my_local_nb_ref  = 1;
559       point1Init();
560       aSatify = isSatisfy();
561       CHECK_MACRO;
562       break;
563
564     case NORM_SEG2:
565       _my_local_ref_dim = 1;
566       _my_local_nb_ref  = 2;
567       seg2aInit();
568       aSatify = isSatisfy();
569       if(!aSatify)
570         {
571           seg2bInit();
572           aSatify = isSatisfy();
573           CHECK_MACRO;
574         }
575       break;
576
577     case NORM_SEG3:
578       _my_local_ref_dim = 1;
579       _my_local_nb_ref  = 3;
580       seg3Init();
581       aSatify = isSatisfy();
582       CHECK_MACRO;
583       break;
584
585     case NORM_SEG4:
586       _my_local_ref_dim = 1;
587       _my_local_nb_ref  = 4;
588       seg4Init();
589       aSatify = isSatisfy();
590       CHECK_MACRO;
591       break;
592
593     case NORM_TRI3:
594       _my_local_ref_dim = 2;
595       _my_local_nb_ref  = 3;
596       tria3aInit();
597       aSatify = isSatisfy();
598
599       if(!aSatify)
600         {
601           tria3bInit();
602           aSatify = isSatisfy();
603           CHECK_MACRO;
604         }
605       break;
606
607     case NORM_TRI6:
608       _my_local_ref_dim = 2;
609       _my_local_nb_ref  = 6;
610       tria6aInit();
611       aSatify = isSatisfy();
612       if(!aSatify)
613         {
614           tria6bInit();
615           aSatify = isSatisfy();
616           CHECK_MACRO;
617         }
618       break;
619
620     case NORM_TRI7:
621       _my_local_ref_dim = 2;
622       _my_local_nb_ref  = 7;
623       tria7aInit();
624       aSatify = isSatisfy();
625       CHECK_MACRO;
626       break;
627
628     case NORM_QUAD4:
629       {
630         _my_local_ref_dim = 2;
631         _my_local_nb_ref  = 4;
632         MapToShapeFunction QUAD4PTR[]={Quad4aInit,Quad4bInit,Quad4cInit,Quad4DegSeg2Init};
633         std::size_t NB_OF_QUAD4PTR(sizeof(QUAD4PTR)/sizeof(MapToShapeFunction));
634         for(std::size_t i=0;i<NB_OF_QUAD4PTR && !aSatify;i++)
635           {
636             (QUAD4PTR[i])(*this);
637             aSatify = isSatisfy();
638           }
639         CHECK_MACRO;
640         break;
641       }
642       break;
643
644     case NORM_QUAD8:
645       _my_local_ref_dim = 2;
646       _my_local_nb_ref  = 8;
647       quad8aInit();
648       aSatify = isSatisfy();
649
650       if(!aSatify)
651         {
652           quad8bInit();
653           aSatify = isSatisfy();
654           CHECK_MACRO;
655         }
656       break;
657
658     case NORM_QUAD9:
659       _my_local_ref_dim = 2;
660       _my_local_nb_ref  = 9;
661       quad9aInit();
662       aSatify = isSatisfy();
663       CHECK_MACRO;
664       break;
665
666     case NORM_TETRA4:
667       _my_local_ref_dim = 3;
668       _my_local_nb_ref  = 4;
669       tetra4aInit();
670       aSatify = isSatisfy();
671
672       if(!aSatify)
673         {
674           tetra4bInit();
675           aSatify = isSatisfy();
676           CHECK_MACRO;
677         }
678       break;
679
680     case NORM_TETRA10:
681       _my_local_ref_dim = 3;
682       _my_local_nb_ref  = 10;
683       tetra10aInit();
684       aSatify = isSatisfy();
685
686       if(!aSatify)
687         {
688           tetra10bInit();
689           aSatify = isSatisfy();
690           CHECK_MACRO;
691         }
692       break;
693
694     case NORM_PYRA5:
695       _my_local_ref_dim = 3;
696       _my_local_nb_ref  = 5;
697       pyra5aInit();
698       aSatify = isSatisfy();
699
700       if(!aSatify)
701         {
702           pyra5bInit();
703           aSatify = isSatisfy();
704           CHECK_MACRO;
705         }
706       break;
707
708     case NORM_PYRA13:
709       _my_local_ref_dim = 3;
710       _my_local_nb_ref  = 13;
711       pyra13aInit();
712       aSatify = isSatisfy();
713
714       if(!aSatify)
715         {
716           pyra13bInit();
717           aSatify = isSatisfy();
718           CHECK_MACRO;
719         }
720       break;
721
722     case NORM_PENTA6:
723       {
724         _my_local_ref_dim = 3;
725         _my_local_nb_ref  = 6;
726         MapToShapeFunction PENTA6PTR[]={Penta6aInit,Penta6bInit,Penta6DegTria3aInit,Penta6DegTria3bInit};
727         std::size_t NB_OF_PENTA6PTR(sizeof(PENTA6PTR)/sizeof(MapToShapeFunction));
728         for(std::size_t i=0;i<NB_OF_PENTA6PTR && !aSatify;i++)
729           {
730             (PENTA6PTR[i])(*this);
731             aSatify = isSatisfy();
732           }
733         CHECK_MACRO;
734         break;
735       }
736
737 /*      _my_local_ref_dim = 3;
738       _my_local_nb_ref  = 6;
739       penta6aInit();
740       aSatify = isSatisfy();
741
742       if(!aSatify)
743         {
744           penta6bInit();
745           aSatify = isSatisfy();
746           CHECK_MACRO;
747         }
748       break;
749 */
750     case NORM_PENTA15:
751       {
752         _my_local_ref_dim = 3;
753         _my_local_nb_ref  = 15;
754         MapToShapeFunction PENTA15PTR[]={Penta15aInit,Penta15bInit};
755         std::size_t NB_OF_PENTA15PTR(sizeof(PENTA15PTR)/sizeof(MapToShapeFunction));
756         for(std::size_t i=0;i<NB_OF_PENTA15PTR && !aSatify;i++)
757           {
758             (PENTA15PTR[i])(*this);
759             aSatify = isSatisfy();
760           }
761         CHECK_MACRO;
762         break;
763       }
764
765     case NORM_PENTA18:
766       {
767         _my_local_ref_dim = 3;
768         _my_local_nb_ref  = 18;
769         MapToShapeFunction PENTA18PTR[]={Penta18aInit,Penta18bInit};
770         std::size_t NB_OF_PENTA18PTR(sizeof(PENTA18PTR)/sizeof(MapToShapeFunction));
771         for(std::size_t i=0;i<NB_OF_PENTA18PTR && !aSatify;i++)
772           {
773             (PENTA18PTR[i])(*this);
774             aSatify = isSatisfy();
775           }
776         CHECK_MACRO;
777         break;
778       }
779
780     case NORM_HEXA8:
781       {
782         _my_local_ref_dim = 3;
783         _my_local_nb_ref  = 8;
784         MapToShapeFunction HEXA8PTR[]={Hexa8aInit,Hexa8bInit,Hexa8DegQuad4aInit,Hexa8DegQuad4bInit,Hexa8DegQuad4cInit};
785         std::size_t NB_OF_HEXA8PTR(sizeof(HEXA8PTR)/sizeof(MapToShapeFunction));
786         for(std::size_t i=0;i<NB_OF_HEXA8PTR && !aSatify;i++)
787           {
788             (HEXA8PTR[i])(*this);
789             aSatify = isSatisfy();
790           }
791         CHECK_MACRO;
792         break;
793       }
794
795     case NORM_HEXA20:
796       _my_local_ref_dim = 3;
797       _my_local_nb_ref  = 20;
798       hexa20aInit();
799       aSatify = isSatisfy();
800
801       if(!aSatify)
802         {
803           hexa20bInit();
804           aSatify = isSatisfy();
805           CHECK_MACRO;
806         }
807       break;
808
809     case NORM_HEXA27:
810       _my_local_ref_dim = 3;
811       _my_local_nb_ref  = 27;
812       hexa27aInit();
813       aSatify = isSatisfy();
814       CHECK_MACRO
815       break;
816
817     default:
818       throw INTERP_KERNEL::Exception("Not managed cell type !");
819       break;
820     }
821 }
822
823 /**
824  * Return shape function value by node id
825  */
826 const double *GaussInfo::getFunctionValues( const int theGaussId ) const 
827 {
828   return _my_function_value.data() + _my_nb_ref*theGaussId ;
829 }
830
831 /**
832  * Return the derivative of shape function value by node id
833  */
834 const double *GaussInfo::getDerivativeOfShapeFunctionAt( const int theGaussId ) const 
835 {
836   return _my_derivative_func_value.data() + _my_nb_ref*getReferenceCoordDim()*theGaussId;
837 }
838
839 void GaussInfo::point1Init()
840 {
841   double *funValue(&_my_function_value[0]);
842   funValue[0] = 1. ;
843 }
844
845 /*!
846  * Init Segment 2 Reference coordinates ans Shape function.
847  */
848 void GaussInfo::seg2aInit() 
849 {
850   LOCAL_COORD_MACRO_BEGIN;
851   case 0:
852     coords[0] = SEG2A_REF[0];
853     break;
854   case 1:
855     coords[0] = SEG2A_REF[1];
856     break;
857   LOCAL_COORD_MACRO_END;
858   
859   SHAPE_FUN_MACRO_BEGIN;
860   funValue[0] = 0.5*(1.0 - gc[0]);
861   funValue[1] = 0.5*(1.0 + gc[0]);
862   SHAPE_FUN_MACRO_END;
863   
864   DEV_SHAPE_FUN_MACRO_BEGIN;
865
866   devFunValue[0] = -0.5;
867
868   devFunValue[1] = 0.5;
869
870   DEV_SHAPE_FUN_MACRO_END;
871 }
872
873 void GaussInfo::seg2bInit() 
874 {
875   LOCAL_COORD_MACRO_BEGIN;
876   case 0:
877     coords[0] = SEG2B_REF[0];
878     break;
879   case 1:
880     coords[0] = SEG2B_REF[1];
881     break;
882   LOCAL_COORD_MACRO_END;
883   
884   SHAPE_FUN_MACRO_BEGIN;
885   funValue[0] = 1.0 - gc[0];
886   funValue[1] = gc[0];
887   SHAPE_FUN_MACRO_END;
888   
889   DEV_SHAPE_FUN_MACRO_BEGIN;
890
891   devFunValue[0] = -1.0 ;
892
893   devFunValue[1] = 1.0 ;
894
895   DEV_SHAPE_FUN_MACRO_END;
896 }
897
898 /*!
899  * Init Segment 3 Reference coordinates ans Shape function.
900  */
901 void GaussInfo::seg3Init() 
902 {
903   LOCAL_COORD_MACRO_BEGIN;
904   case 0:
905     coords[0] = SEG3_REF[0];
906     break;
907   case 1:
908     coords[0] = SEG3_REF[1];
909     break;
910   case 2:
911     coords[0] = SEG3_REF[2];
912   LOCAL_COORD_MACRO_END;
913   
914   SHAPE_FUN_MACRO_BEGIN;
915   funValue[0] = -0.5*(1.0 - gc[0])*gc[0];
916   funValue[1] = 0.5*(1.0 + gc[0])*gc[0];
917   funValue[2] = (1.0 + gc[0])*(1.0 - gc[0]);
918   SHAPE_FUN_MACRO_END;
919
920   DEV_SHAPE_FUN_MACRO_BEGIN;
921   devFunValue[0] = -0.5*(1-2*gc[0]);
922   devFunValue[1] = 0.5*(2*gc[0]+1);
923   devFunValue[2] = -2*gc[0];
924   DEV_SHAPE_FUN_MACRO_END;
925 }
926
927 /*!
928  * Init Segment 4 Reference coordinates ans Shape function.
929  */
930 void GaussInfo::seg4Init() 
931 {
932   LOCAL_COORD_MACRO_BEGIN;
933   case 0:
934     coords[0] = SEG4_REF[0];
935     break;
936   case 1:
937     coords[0] = SEG4_REF[1];
938     break;
939   case 2:
940     coords[0] = SEG4_REF[2];
941     break;
942   case 3:
943     coords[0] = SEG4_REF[3];
944   LOCAL_COORD_MACRO_END;
945   
946   SHAPE_FUN_MACRO_BEGIN;
947   funValue[0] = 9.0/16.0 * (1-gc[0])*(gc[0]+1.0/3.0)*(gc[0]-1.0/3.0);
948   funValue[1] = -9.0/16.0 * (1+gc[0])*(1.0/3.0-gc[0])*(gc[0]+1.0/3.0);
949   funValue[2] = 27.0/16.0 * (gc[0]-1)*(gc[0]+1)*(gc[0]-1.0/3.0);
950   funValue[3] = -27.0/16.0 * (gc[0]-1)*(gc[0]+1)*(gc[0]+1.0/3.0);
951   SHAPE_FUN_MACRO_END;
952
953   DEV_SHAPE_FUN_MACRO_BEGIN;
954   devFunValue[0] = 9.0/16.0 * (-3.0*gc[0]*gc[0]+2.0*gc[0]+1.0/9.0);
955   devFunValue[1] = 9.0/16.0 * (3.0*gc[0]*gc[0]+2.0*gc[0]-1.0/9.0);
956   devFunValue[2] = 27.0/16.0 * (3.0*gc[0]*gc[0]-2.0/3.0*gc[0]-1.0);
957   devFunValue[3] = -27.0/16.0 * (3.0*gc[0]*gc[0]+2.0/3.0*gc[0]-1.0);
958   DEV_SHAPE_FUN_MACRO_END;
959 }
960
961 /*!
962  * Init Triangle Reference coordinates ans Shape function.
963  * Case A.
964  */
965 void GaussInfo::tria3aInit() 
966 {
967   LOCAL_COORD_MACRO_BEGIN;
968   case 0:
969     coords[0] = TRIA3A_REF[0];
970     coords[1] = TRIA3A_REF[1];
971     break;
972   case 1:
973     coords[0] = TRIA3A_REF[2];
974     coords[1] = TRIA3A_REF[3];
975     break;
976   case 2:
977     coords[0] = TRIA3A_REF[4];
978     coords[1] = TRIA3A_REF[5];
979     break;
980   LOCAL_COORD_MACRO_END;
981   
982   SHAPE_FUN_MACRO_BEGIN;
983   funValue[0] = 0.5*(1.0 + gc[1]);
984   funValue[1] = -0.5*(gc[0] + gc[1]);
985   funValue[2] = 0.5*(1.0 + gc[0]);
986   SHAPE_FUN_MACRO_END;
987   
988   DEV_SHAPE_FUN_MACRO_BEGIN;
989
990   devFunValue[0] = 0.0 ;
991   devFunValue[1] = 0.5 ;
992
993   devFunValue[2] = -0.5;
994   devFunValue[3] = -0.5;
995
996   devFunValue[4] = 0.5;
997   devFunValue[5] = 0.0;
998
999   DEV_SHAPE_FUN_MACRO_END;
1000 }
1001
1002 /*!
1003  * Init Triangle Reference coordinates ans Shape function.
1004  * Case B.
1005  */
1006 void GaussInfo::tria3bInit() 
1007 {
1008   LOCAL_COORD_MACRO_BEGIN;
1009   case 0:
1010     coords[0] = TRIA3B_REF[0];
1011     coords[1] = TRIA3B_REF[1];
1012     break;
1013   case 1:
1014     coords[0] = TRIA3B_REF[2];
1015     coords[1] = TRIA3B_REF[3];
1016     break;
1017   case 2:
1018     coords[0] = TRIA3B_REF[4];
1019     coords[1] = TRIA3B_REF[5];
1020     break;
1021   LOCAL_COORD_MACRO_END;
1022   
1023   SHAPE_FUN_MACRO_BEGIN;
1024   funValue[0] = 1.0 - gc[0] - gc[1];
1025   funValue[1] = gc[0];
1026   funValue[2] = gc[1];
1027   SHAPE_FUN_MACRO_END;
1028
1029   DEV_SHAPE_FUN_MACRO_BEGIN;
1030
1031   devFunValue[0] = -1.0 ;
1032   devFunValue[1] = -1.0 ;
1033
1034   devFunValue[2] = 1.0 ;
1035   devFunValue[3] = 0.0 ;
1036
1037   devFunValue[4] = 0.0 ;
1038   devFunValue[5] = 1.0 ;
1039
1040   DEV_SHAPE_FUN_MACRO_END;
1041 }
1042
1043 /*!
1044  * Init Quadratic Triangle Reference coordinates ans Shape function.
1045  * Case A.
1046  */
1047 void GaussInfo::tria6aInit() 
1048 {
1049   LOCAL_COORD_MACRO_BEGIN;
1050   case 0:
1051     coords[0] = TRIA6A_REF[0];
1052     coords[1] = TRIA6A_REF[1];
1053     break;
1054   case 1:
1055     coords[0] = TRIA6A_REF[2];
1056     coords[1] = TRIA6A_REF[3];
1057     break;
1058   case 2:
1059     coords[0] = TRIA6A_REF[4];
1060     coords[1] = TRIA6A_REF[5];
1061     break;
1062   case 3:
1063     coords[0] = TRIA6A_REF[6];
1064     coords[1] = TRIA6A_REF[7];
1065     break;
1066   case 4:
1067     coords[0] = TRIA6A_REF[8];
1068     coords[1] = TRIA6A_REF[9];
1069     break;
1070   case 5:
1071     coords[0] = TRIA6A_REF[10];
1072     coords[1] = TRIA6A_REF[11];
1073     break;
1074   LOCAL_COORD_MACRO_END;
1075   
1076   SHAPE_FUN_MACRO_BEGIN;
1077   funValue[0] = 0.5*(1.0 + gc[1])*gc[1];
1078   funValue[1] = 0.5*(gc[0] + gc[1])*(gc[0] + gc[1] + 1);
1079   funValue[2] = 0.5*(1.0 + gc[0])*gc[0];
1080   funValue[3] = -1.0*(1.0 + gc[1])*(gc[0] + gc[1]);
1081   funValue[4] = -1.0*(1.0 + gc[0])*(gc[0] + gc[1]);
1082   funValue[5] = (1.0 + gc[1])*(1.0 + gc[1]);
1083   SHAPE_FUN_MACRO_END;
1084   
1085   DEV_SHAPE_FUN_MACRO_BEGIN;
1086
1087   devFunValue[0] = 0.0;
1088   devFunValue[1] = 0.5*( 2*gc[1] + 1.0 );
1089
1090   devFunValue[2] = 0.5*( 2*gc[0] + 2.0*gc[1] + 1.0);
1091   devFunValue[3] = 0.5*( 2*gc[1] + 2.0*gc[0] + 1.0);
1092
1093   devFunValue[4] = gc[0] + 0.5;
1094   devFunValue[5] = 0.0;
1095
1096   devFunValue[6] = -1.0*(1.0 + gc[1]);
1097   devFunValue[7] = -1.0*(2*gc[1]+gc[0]+1.0);
1098
1099   devFunValue[8] = -1.0*(2*gc[0]+gc[1]+1.0);
1100   devFunValue[9] = -1.0*(1.0 + gc[0]);
1101
1102   devFunValue[10] = 0.0;
1103   devFunValue[11] = (2*gc[1]+2.0);
1104
1105   DEV_SHAPE_FUN_MACRO_END;
1106 }
1107
1108 /*!
1109  * Init Quadratic Triangle Reference coordinates ans Shape function.
1110  * Case B.
1111  */
1112 void GaussInfo::tria6bInit() 
1113 {
1114   LOCAL_COORD_MACRO_BEGIN;
1115   case 0:
1116     coords[0] = TRIA6B_REF[0];
1117     coords[1] = TRIA6B_REF[1];
1118     break;
1119   case 1:
1120     coords[0] = TRIA6B_REF[2];
1121     coords[1] = TRIA6B_REF[3];
1122     break;
1123   case 2:
1124     coords[0] = TRIA6B_REF[4];
1125     coords[1] = TRIA6B_REF[5];
1126     break;
1127   case 3:
1128     coords[0] = TRIA6B_REF[6];
1129     coords[1] = TRIA6B_REF[7];
1130     break;
1131   case 4:
1132     coords[0] = TRIA6B_REF[8];
1133     coords[1] = TRIA6B_REF[9];
1134     break;
1135   case 5:
1136     coords[0] = TRIA6B_REF[10];
1137     coords[1] = TRIA6B_REF[11];
1138     break;
1139   LOCAL_COORD_MACRO_END;
1140   
1141   SHAPE_FUN_MACRO_BEGIN;
1142   funValue[0] = (1.0 - gc[0] - gc[1])*(1.0 - 2.0*gc[0] - 2.0*gc[1]);
1143   funValue[1] = gc[0]*(2.0*gc[0] - 1.0);
1144   funValue[2] = gc[1]*(2.0*gc[1] - 1.0);
1145   funValue[3] = 4.0*gc[0]*(1.0 - gc[0] - gc[1]);
1146   funValue[4] = 4.0*gc[0]*gc[1];
1147   funValue[5] = 4.0*gc[1]*(1.0 - gc[0] - gc[1]);
1148   SHAPE_FUN_MACRO_END;
1149
1150   DEV_SHAPE_FUN_MACRO_BEGIN;
1151
1152   devFunValue[0] = 4*gc[0] + 4*gc[1] - 3.0 ;
1153   devFunValue[1] = 4*gc[1] + 4*gc[0] - 3.0 ;
1154
1155   devFunValue[2] = 4*gc[0] - 1.0 ;
1156   devFunValue[3] = 0.0 ;
1157
1158   devFunValue[4] = 0.0 ;
1159   devFunValue[5] = 4*gc[1] - 1.0 ;
1160
1161   devFunValue[6] = -8.0*gc[0] - 4.0 * gc[1] + 4.0 ;
1162   devFunValue[7] = -4.0*gc[0] ; 
1163
1164   devFunValue[8] = 4.0*gc[1];
1165   devFunValue[9] = 4.0*gc[0];
1166
1167   devFunValue[10] = -4.0*gc[1] ;
1168   devFunValue[11] = -8.0*gc[1] - 4.0*gc[0] + 4.0 ; 
1169
1170   DEV_SHAPE_FUN_MACRO_END;
1171 }
1172
1173 void GaussInfo::tria7aInit()
1174 {
1175   LOCAL_COORD_MACRO_BEGIN;
1176   case 0:
1177     coords[0] = TRIA7A_REF[0];
1178     coords[1] = TRIA7A_REF[1];
1179     break;
1180   case 1:
1181     coords[0] = TRIA7A_REF[2];
1182     coords[1] = TRIA7A_REF[3];
1183     break;
1184   case 2:
1185     coords[0] = TRIA7A_REF[4];
1186     coords[1] = TRIA7A_REF[5];
1187     break;
1188   case 3:
1189     coords[0] = TRIA7A_REF[6];
1190     coords[1] = TRIA7A_REF[7];
1191     break;
1192   case 4:
1193     coords[0] = TRIA7A_REF[8];
1194     coords[1] = TRIA7A_REF[9];
1195     break;
1196   case 5:
1197     coords[0] = TRIA7A_REF[10];
1198     coords[1] = TRIA7A_REF[11];
1199     break;
1200   case 6:
1201     coords[0] = TRIA7A_REF[12];
1202     coords[1] = TRIA7A_REF[13];
1203     break;
1204   LOCAL_COORD_MACRO_END;
1205
1206   SHAPE_FUN_MACRO_BEGIN;
1207   funValue[0]=1-3*(gc[0]+gc[1])+2*(gc[0]*gc[0]+gc[1]*gc[1])+7*gc[0]*gc[1]-3*gc[0]*gc[1]*(gc[0]+gc[1]);
1208   funValue[1]=gc[0]*(-1+2*gc[0]+3*gc[1]-3*gc[1]*(gc[0]+gc[1]));
1209   funValue[2]=gc[1]*(-1.+3.*gc[0]+2.*gc[1]-3.*gc[0]*(gc[0]+gc[1]));
1210   funValue[3]=4*gc[0]*(1-gc[0]-4*gc[1]+3*gc[1]*(gc[0]+gc[1]));
1211   funValue[4]=4*gc[0]*gc[1]*(-2+3*(gc[0]+gc[1]));
1212   funValue[5]=4*gc[1]*(1-4*gc[0]-gc[1]+3*gc[0]*(gc[0]+gc[1]));
1213   funValue[6]=27*gc[0]*gc[1]*(1-gc[0]-gc[1]);
1214   SHAPE_FUN_MACRO_END;
1215   
1216   DEV_SHAPE_FUN_MACRO_BEGIN;
1217   devFunValue[0] = 7*gc[1]-3*gc[1]*(2*gc[0]+gc[1])+4*gc[0]-3;
1218   devFunValue[1] = -3+4.0*gc[1]+7.0*gc[0]-6.0*gc[0]*gc[1]-3.0*gc[0]*gc[0];
1219   
1220   devFunValue[2] = -6*gc[0]*gc[1]+4*gc[0]-3*gc[1]*gc[1]+3*gc[1]-1;
1221   devFunValue[3] = 3.0*gc[0]*(1-2.0*gc[1]-gc[0]);
1222   
1223   devFunValue[4] = 3.0*gc[1]*(1-2.0*gc[0]-gc[1]);
1224   devFunValue[5] = -1+4.0*gc[1]+3.0*gc[0]-6.0*gc[0]*gc[1]-3.0*gc[0]*gc[0];
1225  
1226   devFunValue[6] = 4*(6*gc[0]*gc[1]-2*gc[0]+3*gc[1]*gc[1]-4*gc[1]+1);
1227   devFunValue[7] = 4.0*gc[0]*(-4.0+6.0*gc[1]+3.0*gc[0]);
1228   
1229   devFunValue[8] = 4*gc[1]*(6*gc[0]+3*gc[1]-2);
1230   devFunValue[9] = 4.0*gc[0]*(-2.0+6.0*gc[1]+3.0*gc[0]);
1231   
1232   devFunValue[10] = 4*gc[1]*(3*(gc[1]+2*gc[0])-4);
1233   devFunValue[11] = 4.0*(1-2.0*gc[1]-4.0*gc[0]+6.0*gc[0]*gc[1]+3.0*gc[0]*gc[0]);
1234   
1235   devFunValue[12] = 27*gc[1]*(-2*gc[0]-gc[1]+1);
1236   devFunValue[13] = 27.*gc[0]*(1-2.0*gc[1]-gc[0]);
1237   DEV_SHAPE_FUN_MACRO_END;
1238 }
1239
1240 /*!
1241  * Init Quadrangle Reference coordinates ans Shape function.
1242  * Case A.
1243  */
1244 void GaussInfo::quad4aInit() 
1245 {
1246   LOCAL_COORD_MACRO_BEGIN;
1247   case 0:
1248     coords[0] = QUAD4A_REF[0];
1249     coords[1] = QUAD4A_REF[1];
1250     break;
1251   case 1:
1252     coords[0] = QUAD4A_REF[2];
1253     coords[1] = QUAD4A_REF[3];
1254     break;
1255   case 2:
1256     coords[0] = QUAD4A_REF[4];
1257     coords[1] = QUAD4A_REF[5];
1258     break;
1259   case 3:
1260     coords[0] = QUAD4A_REF[6];
1261     coords[1] = QUAD4A_REF[7];
1262     break;
1263   LOCAL_COORD_MACRO_END;
1264   
1265   SHAPE_FUN_MACRO_BEGIN;
1266   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
1267   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
1268   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
1269   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0]);
1270   SHAPE_FUN_MACRO_END;
1271   
1272   DEV_SHAPE_FUN_MACRO_BEGIN;
1273
1274   devFunValue[0] = -0.25*(1.0 + gc[1]);
1275   devFunValue[1] = 0.25*(1.0 - gc[0]);
1276
1277   devFunValue[2] = -0.25*(1.0 - gc[1]);
1278   devFunValue[3] = -0.25*(1.0 - gc[0]);
1279
1280   devFunValue[4] = 0.25*(1.0 - gc[1]);
1281   devFunValue[5] = -0.25*(1.0 + gc[0]);
1282
1283   devFunValue[6] = 0.25*(1.0 + gc[1]);
1284   devFunValue[7] = 0.25*(1.0 + gc[0]);
1285   
1286   DEV_SHAPE_FUN_MACRO_END;
1287 }
1288
1289 /*!
1290  * Init Quadrangle Reference coordinates ans Shape function.
1291  * Case B.
1292  */
1293 void GaussInfo::quad4bInit() 
1294 {
1295   LOCAL_COORD_MACRO_BEGIN;
1296   case 0:
1297     coords[0] = QUAD4B_REF[0];
1298     coords[1] = QUAD4B_REF[1];
1299     break;
1300   case 1:
1301     coords[0] = QUAD4B_REF[2];
1302     coords[1] = QUAD4B_REF[3];
1303     break;
1304   case 2:
1305     coords[0] = QUAD4B_REF[4];
1306     coords[1] = QUAD4B_REF[5];
1307     break;
1308   case 3:
1309     coords[0] = QUAD4B_REF[6];
1310     coords[1] = QUAD4B_REF[7];
1311     break;
1312   LOCAL_COORD_MACRO_END;
1313   
1314   SHAPE_FUN_MACRO_BEGIN;
1315   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1316   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1317   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1318   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1319   SHAPE_FUN_MACRO_END;
1320   
1321   DEV_SHAPE_FUN_MACRO_BEGIN;
1322
1323   devFunValue[0] = -0.25*(1.0 - gc[1]);
1324   devFunValue[1] = -0.25*(1.0 - gc[0]);
1325
1326   devFunValue[2] = 0.25*(1.0 - gc[1]);
1327   devFunValue[3] = -0.25*(1.0 + gc[0]);
1328
1329   devFunValue[4] = 0.25*(1.0 + gc[1]);
1330   devFunValue[5] = 0.25*(1.0 + gc[0]);
1331
1332   devFunValue[6] = -0.25*(1.0 + gc[1]);
1333   devFunValue[7] = 0.25*(1.0 - gc[0]);
1334
1335   DEV_SHAPE_FUN_MACRO_END;
1336 }
1337
1338 void GaussInfo::quad4cInit() 
1339 {
1340   LOCAL_COORD_MACRO_BEGIN;
1341  case  0:
1342    coords[0] = -1.0;
1343    coords[1] = -1.0;
1344    break;
1345  case  1:
1346    coords[0] = -1.0;
1347    coords[1] =  1.0;
1348    break;
1349  case  2:
1350    coords[0] =  1.0;
1351    coords[1] =  1.0;
1352    break;
1353  case  3:
1354    coords[0] =  1.0;
1355    coords[1] = -1.0;
1356    break;
1357
1358    LOCAL_COORD_MACRO_END;
1359
1360    SHAPE_FUN_MACRO_BEGIN;
1361    funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
1362    funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
1363    funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
1364    funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
1365    SHAPE_FUN_MACRO_END;
1366    
1367   DEV_SHAPE_FUN_MACRO_BEGIN;
1368
1369   devFunValue[0] = -0.25*(1.0 - gc[1]);
1370   devFunValue[1] = -0.25*(1.0 - gc[0]);
1371
1372   devFunValue[2] = -0.25*(1.0 + gc[1]);
1373   devFunValue[3] = 0.25*(1.0 - gc[0]);
1374
1375   devFunValue[4] = 0.25*(1.0 + gc[1]);
1376   devFunValue[5] = 0.25*(1.0 + gc[0]);
1377
1378   devFunValue[6] = 0.25*(1.0 - gc[1]);
1379   devFunValue[7] = -0.25*(1.0 + gc[0]);
1380
1381   DEV_SHAPE_FUN_MACRO_END;
1382 }
1383
1384 /*!
1385  * This shapefunc map is same as degenerated seg2aInit
1386  */
1387 void GaussInfo::quad4DegSeg2Init()
1388 {
1389   LOCAL_COORD_MACRO_BEGIN;
1390  case  0:
1391    coords[0] = -1.0;
1392    coords[1] =  0.0;
1393    break;
1394  case  1:
1395    coords[0] =  1.0;
1396    coords[1] =  0.0;
1397    break;
1398  case  2:
1399    coords[0] =  0.0;
1400    coords[1] =  0.0;
1401    break;
1402  case  3:
1403    coords[0] =  0.0;
1404    coords[1] =  0.0;
1405    break;
1406    LOCAL_COORD_MACRO_END;
1407
1408    SHAPE_FUN_MACRO_BEGIN;
1409    funValue[0] = 0.5*(1.0 - gc[0]);
1410    funValue[1] = 0.5*(1.0 + gc[0]);
1411    funValue[2] = 0.;
1412    funValue[3] = 0.;
1413    SHAPE_FUN_MACRO_END;
1414    
1415   DEV_SHAPE_FUN_MACRO_BEGIN;
1416
1417   devFunValue[0] = -0.5;
1418   devFunValue[1] = 0.0;
1419
1420   devFunValue[2] = 0.5;
1421   devFunValue[3] = 0.0;
1422
1423   devFunValue[4] = 0.0;
1424   devFunValue[5] = 0.0;
1425
1426   devFunValue[6] = 0.0;
1427   devFunValue[7] = 0.0;
1428
1429   DEV_SHAPE_FUN_MACRO_END;
1430 }
1431
1432 /*!
1433  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1434  * Case A.
1435  */
1436 void GaussInfo::quad8aInit() 
1437 {
1438   LOCAL_COORD_MACRO_BEGIN;
1439   case 0:
1440     coords[0] = QUAD8A_REF[0];
1441     coords[1] = QUAD8A_REF[1];
1442     break;
1443   case 1:
1444     coords[0] = QUAD8A_REF[2];
1445     coords[1] = QUAD8A_REF[3];
1446     break;
1447   case 2:
1448     coords[0] = QUAD8A_REF[4];
1449     coords[1] = QUAD8A_REF[5];
1450     break;
1451   case 3:
1452     coords[0] = QUAD8A_REF[6];
1453     coords[1] = QUAD8A_REF[7];
1454     break;
1455   case 4:
1456     coords[0] = QUAD8A_REF[8];
1457     coords[1] = QUAD8A_REF[9];
1458     break;
1459   case 5:
1460     coords[0] = QUAD8A_REF[10];
1461     coords[1] = QUAD8A_REF[11];
1462     break;
1463   case 6:
1464     coords[0] = QUAD8A_REF[12];
1465     coords[1] = QUAD8A_REF[13];
1466     break;
1467   case 7:
1468     coords[0] = QUAD8A_REF[14];
1469     coords[1] = QUAD8A_REF[15];
1470     break;
1471   LOCAL_COORD_MACRO_END;
1472   
1473   SHAPE_FUN_MACRO_BEGIN;
1474   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0])*(gc[1] - gc[0] - 1.0);
1475   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0])*(-gc[1] - gc[0] - 1.0);
1476   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0])*(-gc[1] + gc[0] - 1.0);
1477   funValue[3] = 0.25*(1.0 + gc[1])*(1.0 + gc[0])*(gc[1] + gc[0] - 1.0);
1478   funValue[4] = 0.5*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1479   funValue[5] = 0.5*(1.0 - gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1480   funValue[6] = 0.5*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[1]);
1481   funValue[7] = 0.5*(1.0 + gc[1])*(1.0 - gc[0])*(1.0 + gc[0]);
1482   SHAPE_FUN_MACRO_END;
1483
1484   DEV_SHAPE_FUN_MACRO_BEGIN;
1485
1486   devFunValue[0] = 0.25*(1.0 + gc[1])*(2*gc[0]-gc[1]);
1487   devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1]-gc[0]);
1488
1489   devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0]+gc[1]);
1490   devFunValue[3] = 0.25*(1.0 - gc[0])*(2*gc[1]+gc[0]);
1491
1492   devFunValue[4] = 0.25*(1.0 - gc[1])*(2*gc[0]-gc[1]);
1493   devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1]-gc[0]);
1494
1495   devFunValue[6] = 0.25*(1.0 + gc[1])*(2*gc[0]+gc[1]);
1496   devFunValue[7] = 0.25*(1.0 + gc[0])*(2*gc[1]+gc[0]);
1497
1498   devFunValue[8] = -0.5*(1.0 - gc[1])*(1.0 + gc[1]);
1499   devFunValue[9] = 0.5*(1.0 - gc[0])*(-2*gc[1]);
1500
1501   devFunValue[10] = 0.5*(1.0 - gc[1])*(-2*gc[0]);
1502   devFunValue[11] = -0.5*(1.0 - gc[0])*(1.0 + gc[0]);
1503
1504   devFunValue[12] = 0.5*(1.0 - gc[1])*(1.0 + gc[1]);
1505   devFunValue[13] = 0.5*(1.0 + gc[0])*(-2*gc[1]);
1506
1507   devFunValue[14] = 0.5*(1.0 + gc[1])*(-2*gc[0]);
1508   devFunValue[15] = 0.5*(1.0 - gc[0])*(1.0 + gc[0]);
1509
1510   DEV_SHAPE_FUN_MACRO_END;
1511 }
1512
1513 /*!
1514  * Init Quadratic Quadrangle Reference coordinates ans Shape function.
1515  * Case B.
1516  */
1517 void GaussInfo::quad8bInit() 
1518 {
1519   LOCAL_COORD_MACRO_BEGIN;
1520   case 0:
1521     coords[0] = QUAD8B_REF[0];
1522     coords[1] = QUAD8B_REF[1];
1523     break;
1524   case 1:
1525     coords[0] = QUAD8B_REF[2];
1526     coords[1] = QUAD8B_REF[3];
1527     break;
1528   case 2:
1529     coords[0] = QUAD8B_REF[4];
1530     coords[1] = QUAD8B_REF[5];
1531     break;
1532   case 3:
1533     coords[0] = QUAD8B_REF[6];
1534     coords[1] = QUAD8B_REF[7];
1535     break;
1536   case 4:
1537     coords[0] = QUAD8B_REF[8];
1538     coords[1] = QUAD8B_REF[9];
1539     break;
1540   case 5:
1541     coords[0] = QUAD8B_REF[10];
1542     coords[1] = QUAD8B_REF[11];
1543     break;
1544   case 6:
1545     coords[0] = QUAD8B_REF[12];
1546     coords[1] = QUAD8B_REF[13];
1547     break;
1548   case 7:
1549     coords[0] = QUAD8B_REF[14];
1550     coords[1] = QUAD8B_REF[15];
1551     break;
1552   LOCAL_COORD_MACRO_END;
1553   
1554   SHAPE_FUN_MACRO_BEGIN;
1555   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1]);
1556   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1]);
1557   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1]);
1558   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1]);
1559   funValue[4] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
1560   funValue[5] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
1561   funValue[6] = 0.5*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
1562   funValue[7] = 0.5*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
1563   SHAPE_FUN_MACRO_END;
1564   
1565   DEV_SHAPE_FUN_MACRO_BEGIN;
1566
1567   devFunValue[0] = 0.25*(1.0 - gc[1])*(2*gc[0] + gc[1] );
1568   devFunValue[1] = 0.25*(1.0 - gc[0])*(2*gc[1] + gc[0] );
1569
1570   devFunValue[2] = 0.25*(1.0 - gc[1])*(2*gc[0] - gc[1] );
1571   devFunValue[3] = 0.25*(1.0 + gc[0])*(2*gc[1] - gc[0] );
1572
1573   devFunValue[4] = 0.25*(1.0 + gc[1])*(2*gc[0] + gc[1]);
1574   devFunValue[5] = 0.25*(1.0 + gc[0])*(2*gc[1] + gc[0]);
1575
1576   devFunValue[6] = 0.25*( 1.0 + gc[1])*(2*gc[0] - gc[1] );
1577   devFunValue[7] = 0.25*(1.0 - gc[0])*(2*gc[1] - gc[0]);
1578
1579   devFunValue[8] = -(1.0 - gc[1])*gc[0];
1580   devFunValue[9] = -0.5*(1.0 - gc[0]*gc[0]);
1581
1582   devFunValue[10] = 0.5*(1.0 - gc[1]*gc[1]);
1583   devFunValue[11] = -(1.0 + gc[0])*gc[1];
1584
1585   devFunValue[12] = -(1.0 + gc[1])*gc[0];
1586   devFunValue[13] = 0.5*(1.0 - gc[0]*gc[0]);
1587
1588   devFunValue[14] = -0.5*(1.0 - gc[1]*gc[1]);
1589   devFunValue[15] = -(1.0 - gc[0])*gc[1];
1590
1591   DEV_SHAPE_FUN_MACRO_END;
1592 }
1593
1594 void GaussInfo::quad9aInit()
1595 {
1596   LOCAL_COORD_MACRO_BEGIN;
1597   case 0:
1598     coords[0] = QUAD9A_REF[0];
1599     coords[1] = QUAD9A_REF[1];
1600     break;
1601   case 1:
1602     coords[0] = QUAD9A_REF[2];
1603     coords[1] = QUAD9A_REF[3];
1604     break;
1605   case 2:
1606     coords[0] = QUAD9A_REF[4];
1607     coords[1] = QUAD9A_REF[5];
1608     break;
1609   case 3:
1610     coords[0] = QUAD9A_REF[6];
1611     coords[1] = QUAD9A_REF[7];
1612     break;
1613   case 4:
1614     coords[0] = QUAD9A_REF[8];
1615     coords[1] = QUAD9A_REF[9];
1616     break;
1617   case 5:
1618     coords[0] = QUAD9A_REF[10];
1619     coords[1] = QUAD9A_REF[11];
1620     break;
1621   case 6:
1622     coords[0] = QUAD9A_REF[12];
1623     coords[1] = QUAD9A_REF[13];
1624     break;
1625   case 7:
1626     coords[0] = QUAD9A_REF[14];
1627     coords[1] = QUAD9A_REF[15];
1628     break;
1629   case 8:
1630     coords[0] = QUAD9A_REF[16];
1631     coords[1] = QUAD9A_REF[17];
1632     break;
1633   LOCAL_COORD_MACRO_END;
1634   
1635   SHAPE_FUN_MACRO_BEGIN;
1636   funValue[0] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]-1.);
1637   funValue[1] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]-1.);
1638   funValue[2] = 0.25*gc[0]*gc[1]*(gc[0]+1.)*(gc[1]+1.);
1639   funValue[3] = 0.25*gc[0]*gc[1]*(gc[0]-1.)*(gc[1]+1.);
1640   funValue[4] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.);
1641   funValue[5] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1]);
1642   funValue[6] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.);
1643   funValue[7] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1]);
1644   funValue[8] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1]);
1645   SHAPE_FUN_MACRO_END;
1646   
1647   DEV_SHAPE_FUN_MACRO_BEGIN;
1648   devFunValue[0] = 0.5*(2*gc[0]-1.0)*0.5*gc[1]*(gc[1]-1.0);
1649   devFunValue[1] = 0.5*gc[0]*(gc[0]-1.0)*0.5*(2*gc[1]-1.0);;
1650   
1651   devFunValue[2] = 0.5*(2.0*gc[0]+1.0)*0.5*gc[1]*(gc[1]-1.0);
1652   devFunValue[3] = 0.5*gc[0]*(1.0+gc[0])*0.5*(2.0*gc[1]-1.0);
1653   
1654   devFunValue[4] = 0.5*(2.0*gc[0]+1.0)*0.5*gc[1]*(1.0+gc[1]);
1655   devFunValue[5] = 0.5*gc[0]*(1.0+gc[0])*0.5*(2.0*gc[1]+1.0);
1656   
1657   devFunValue[6] = 0.5*(2.0*gc[0]-1.0)*0.5*gc[1]*(1.0+gc[1]);
1658   devFunValue[7] = 0.5*gc[0]*(gc[0]-1.0)*0.5*(2.0*gc[1]+1.0);
1659   
1660   devFunValue[8] = -2.0*gc[0]*0.5*gc[1]*(gc[1]-1.0);
1661   devFunValue[9] = (1.0+gc[0])*(1.0-gc[0])*0.5*(2.0*gc[1]-1.0);
1662   
1663   devFunValue[10] = 0.5*(2.0*gc[0]+1.0)*(1.0+gc[1])*(1.0-gc[1]);
1664   devFunValue[11] = 0.5*gc[0]*(1+gc[0])*-2.0*gc[1];
1665   
1666   devFunValue[12] = -2.0*gc[0]*0.5*gc[1]*(1.0+gc[1]);
1667   devFunValue[13] = (1.0+gc[0])*(1.0-gc[0])*0.5*(2.0*gc[1]+1.0);
1668   
1669   devFunValue[14] = 0.5*(2.0*gc[0]-1.0)*(1.0+gc[1])*(1.0-gc[1]);
1670   devFunValue[15] = 0.5*gc[0]*(gc[0]-1.0)*-2.0*gc[1];
1671   
1672   devFunValue[16] = -2.0*gc[0]*(1.0+gc[1])*(1.0-gc[1]);
1673   devFunValue[17] = (1.0+gc[0])*(1.0-gc[0])*-2.0*gc[1];
1674   DEV_SHAPE_FUN_MACRO_END;
1675 }
1676
1677 /*!
1678  * Init Tetrahedron Reference coordinates ans Shape function.
1679  * Case A.
1680  */
1681 void GaussInfo::tetra4aInit() 
1682 {
1683   LOCAL_COORD_MACRO_BEGIN;
1684   case 0:
1685     coords[0] = TETRA4A_REF[0];
1686     coords[1] = TETRA4A_REF[1];
1687     coords[2] = TETRA4A_REF[2];
1688     break;
1689   case 1:
1690     coords[0] = TETRA4A_REF[3];
1691     coords[1] = TETRA4A_REF[4];
1692     coords[2] = TETRA4A_REF[5];
1693     break;
1694   case 2:
1695     coords[0] = TETRA4A_REF[6];
1696     coords[1] = TETRA4A_REF[7];
1697     coords[2] = TETRA4A_REF[8];
1698     break;
1699   case 3:
1700     coords[0] = TETRA4A_REF[9];
1701     coords[1] = TETRA4A_REF[10];
1702     coords[2] = TETRA4A_REF[11];
1703     break;
1704   LOCAL_COORD_MACRO_END;
1705   
1706   SHAPE_FUN_MACRO_BEGIN;
1707   funValue[0] = gc[1];
1708   funValue[1] = gc[2];
1709   funValue[2] = 1.0 - gc[0] - gc[1] - gc[2];
1710   funValue[3] = gc[0];
1711   SHAPE_FUN_MACRO_END;
1712   
1713   DEV_SHAPE_FUN_MACRO_BEGIN;
1714   devFunValue[0] = 0.0;
1715   devFunValue[1] = 1.0;
1716   devFunValue[2] = 0.0;
1717   
1718   devFunValue[3] = 0.0;
1719   devFunValue[4] = 0.0;
1720   devFunValue[5] = 1.0;
1721   
1722   devFunValue[6] = -1.0;
1723   devFunValue[7] = -1.0;
1724   devFunValue[8] = -1.0;
1725   
1726   devFunValue[9] = 1.0;
1727   devFunValue[10] = 0.0;
1728   devFunValue[11] = 0.0;
1729   DEV_SHAPE_FUN_MACRO_END;
1730 }
1731
1732 /*!
1733  * Init Tetrahedron Reference coordinates ans Shape function.
1734  * Case B.
1735  */
1736 void GaussInfo::tetra4bInit() 
1737 {
1738   LOCAL_COORD_MACRO_BEGIN;
1739   case 0:
1740     coords[0] = TETRA4B_REF[0];
1741     coords[1] = TETRA4B_REF[1];
1742     coords[2] = TETRA4B_REF[2];
1743     break;
1744   case 1:
1745     coords[0] = TETRA4B_REF[3];
1746     coords[1] = TETRA4B_REF[4];
1747     coords[2] = TETRA4B_REF[5];
1748     break;
1749   case 2:
1750     coords[0] = TETRA4B_REF[6];
1751     coords[1] = TETRA4B_REF[7];
1752     coords[2] = TETRA4B_REF[8];
1753     break;
1754   case 3:
1755     coords[0] = TETRA4B_REF[9];
1756     coords[1] = TETRA4B_REF[10];
1757     coords[2] = TETRA4B_REF[11];
1758     break;
1759   LOCAL_COORD_MACRO_END;
1760   
1761   SHAPE_FUN_MACRO_BEGIN;
1762   funValue[0] = gc[1];
1763   funValue[2] = gc[2];
1764   funValue[1] = 1.0 - gc[0] - gc[1] - gc[2];
1765   funValue[3] = gc[0];
1766   SHAPE_FUN_MACRO_END;
1767   
1768   DEV_SHAPE_FUN_MACRO_BEGIN;
1769   devFunValue[0] = 0.0;
1770   devFunValue[1] = 1.0;
1771   devFunValue[2] = 0.0;
1772   
1773   devFunValue[3] = -1.0;
1774   devFunValue[4] = -1.0;
1775   devFunValue[5] = -1.0;
1776   
1777   devFunValue[6] = 0.0;
1778   devFunValue[7] = 0.0;
1779   devFunValue[8] = 1.0;
1780   
1781   devFunValue[9] = 1.0;
1782   devFunValue[10] = 0.0;
1783   devFunValue[11] = 0.0;
1784   DEV_SHAPE_FUN_MACRO_END;
1785 }
1786
1787 /*!
1788  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1789  * Case A.
1790  */
1791 void GaussInfo::tetra10aInit() 
1792 {
1793   LOCAL_COORD_MACRO_BEGIN;
1794   case 0:
1795     coords[0] = TETRA10A_REF[0];
1796     coords[1] = TETRA10A_REF[1];
1797     coords[2] = TETRA10A_REF[2];
1798     break;
1799   case 1:
1800     coords[0] = TETRA10A_REF[3];
1801     coords[1] = TETRA10A_REF[4];
1802     coords[2] = TETRA10A_REF[5];
1803     break;
1804   case 2:
1805     coords[0] = TETRA10A_REF[6];
1806     coords[1] = TETRA10A_REF[7];
1807     coords[2] = TETRA10A_REF[8];
1808     break;
1809   case 3:
1810     coords[0] = TETRA10A_REF[9];
1811     coords[1] = TETRA10A_REF[10];
1812     coords[2] = TETRA10A_REF[11];
1813     break;
1814   case 4:
1815     coords[0] = TETRA10A_REF[12];
1816     coords[1] = TETRA10A_REF[13];
1817     coords[2] = TETRA10A_REF[14];
1818     break;
1819   case 5:
1820     coords[0] = TETRA10A_REF[15];
1821     coords[1] = TETRA10A_REF[16];
1822     coords[2] = TETRA10A_REF[17];
1823     break;
1824   case 6:
1825     coords[0] = TETRA10A_REF[18];
1826     coords[1] = TETRA10A_REF[19];
1827     coords[2] = TETRA10A_REF[20];
1828     break;
1829   case 7:
1830     coords[0] = TETRA10A_REF[21];
1831     coords[1] = TETRA10A_REF[22];
1832     coords[2] = TETRA10A_REF[23];
1833     break;
1834   case 8:
1835     coords[0] = TETRA10A_REF[24];
1836     coords[1] = TETRA10A_REF[25];
1837     coords[2] = TETRA10A_REF[26];
1838     break;
1839   case 9:
1840     coords[0] = TETRA10A_REF[27];
1841     coords[1] = TETRA10A_REF[28];
1842     coords[2] = TETRA10A_REF[29];
1843     break;
1844   LOCAL_COORD_MACRO_END;
1845
1846   SHAPE_FUN_MACRO_BEGIN;
1847   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1848   funValue[1] = gc[2]*(2.0*gc[2] - 1.0);
1849   funValue[2] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1850   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1851   funValue[4] = 4.0*gc[1]*gc[2];
1852   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1853   funValue[6] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1854   funValue[7] = 4.0*gc[0]*gc[1];
1855   funValue[8] = 4.0*gc[0]*gc[2];
1856   funValue[9] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1857   SHAPE_FUN_MACRO_END;
1858   
1859   DEV_SHAPE_FUN_MACRO_BEGIN;
1860   devFunValue[0] = 0.0;
1861   devFunValue[1] = 4.0*gc[1]-1.0;
1862   devFunValue[2] = 0.0;
1863   
1864   devFunValue[3] = 0.0;
1865   devFunValue[4] = 0.0;
1866   devFunValue[5] = 4.0*gc[2]-1.0;
1867   
1868   devFunValue[6] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1869   devFunValue[7] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1870   devFunValue[8] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1871   
1872   devFunValue[9] = 4.0*gc[0]-1.0;
1873   devFunValue[10] = 0.0;
1874   devFunValue[11] = 0.0;
1875   
1876   devFunValue[12] = 0.0;
1877   devFunValue[13] = 4.0*gc[2];
1878   devFunValue[14] = 4.0*gc[1];
1879   
1880   devFunValue[15] = -4.0*gc[2];
1881   devFunValue[16] = -4.0*gc[2];
1882   devFunValue[17] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[2]);
1883   
1884   devFunValue[18] = -4.0*gc[1];
1885   devFunValue[19] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[1]);;
1886   devFunValue[20] = -4.0*gc[1];
1887   
1888   devFunValue[21] = 4.0*gc[1];
1889   devFunValue[22] = 4.0*gc[0];
1890   devFunValue[23] = 0.0;
1891   
1892   devFunValue[24] = 4.0*gc[2];
1893   devFunValue[25] = 0.0;
1894   devFunValue[26] = 4.0*gc[0];
1895   
1896   devFunValue[27] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[0]);
1897   devFunValue[28] = -4.0*gc[0];
1898   devFunValue[29] = -4.0*gc[0];
1899   DEV_SHAPE_FUN_MACRO_END;
1900 }
1901
1902 /*!
1903  * Init Quadratic Tetrahedron Reference coordinates ans Shape function.
1904  * Case B.
1905  */
1906 void GaussInfo::tetra10bInit() 
1907 {
1908   LOCAL_COORD_MACRO_BEGIN;
1909   case 0:
1910     coords[0] = TETRA10B_REF[0];
1911     coords[1] = TETRA10B_REF[1];
1912     coords[2] = TETRA10B_REF[2];
1913     break;
1914   case 1:
1915     coords[0] = TETRA10B_REF[3];
1916     coords[1] = TETRA10B_REF[4];
1917     coords[2] = TETRA10B_REF[5];
1918     break;
1919   case 2:
1920     coords[0] = TETRA10B_REF[6];
1921     coords[1] = TETRA10B_REF[7];
1922     coords[2] = TETRA10B_REF[8];
1923     break;
1924   case 3:
1925     coords[0] = TETRA10B_REF[9];
1926     coords[1] = TETRA10B_REF[10];
1927     coords[2] = TETRA10B_REF[11];
1928     break;
1929   case 4:
1930     coords[0] = TETRA10B_REF[12];
1931     coords[1] = TETRA10B_REF[13];
1932     coords[2] = TETRA10B_REF[14];
1933     break;
1934   case 5:
1935     coords[0] = TETRA10B_REF[15];
1936     coords[1] = TETRA10B_REF[16];
1937     coords[2] = TETRA10B_REF[17];
1938     break;
1939   case 6:
1940     coords[0] = TETRA10B_REF[18];
1941     coords[1] = TETRA10B_REF[19];
1942     coords[2] = TETRA10B_REF[20];
1943     break;
1944   case 7:
1945     coords[0] = TETRA10B_REF[21];
1946     coords[1] = TETRA10B_REF[22];
1947     coords[2] = TETRA10B_REF[23];
1948     break;
1949   case 8:
1950     coords[0] = TETRA10B_REF[24];
1951     coords[1] = TETRA10B_REF[25];
1952     coords[2] = TETRA10B_REF[26];
1953     break;
1954   case 9:
1955     coords[0] = TETRA10B_REF[27];
1956     coords[1] = TETRA10B_REF[28];
1957     coords[2] = TETRA10B_REF[29];
1958     break;
1959   LOCAL_COORD_MACRO_END;
1960   SHAPE_FUN_MACRO_BEGIN;
1961   funValue[0] = gc[1]*(2.0*gc[1] - 1.0);
1962   funValue[2] = gc[2]*(2.0*gc[2] - 1.0);
1963   funValue[1] = (1.0 - gc[0] - gc[1] - gc[2])*(1.0 - 2.0*gc[0] - 2.0*gc[1] - 2.0*gc[2]);
1964   funValue[3] = gc[0]*(2.0*gc[0] - 1.0);
1965   funValue[6] = 4.0*gc[1]*gc[2];
1966   funValue[5] = 4.0*gc[2]*(1.0 - gc[0] - gc[1] - gc[2]);
1967   funValue[4] = 4.0*gc[1]*(1.0 - gc[0] - gc[1] - gc[2]);
1968   funValue[7] = 4.0*gc[0]*gc[1];
1969   funValue[9] = 4.0*gc[0]*gc[2];
1970   funValue[8] = 4.0*gc[0]*(1.0 - gc[0] - gc[1] - gc[2]);
1971   SHAPE_FUN_MACRO_END;
1972   
1973   DEV_SHAPE_FUN_MACRO_BEGIN;
1974   devFunValue[0] = 0.0;
1975   devFunValue[1] = 4.0*gc[1]-1.0;
1976   devFunValue[2] = 0.0;
1977   
1978   devFunValue[3] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1979   devFunValue[4] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1980   devFunValue[5] = 1.0-4.0*(1-gc[0]-gc[1]-gc[2]);
1981   
1982   devFunValue[6] = 0.0;
1983   devFunValue[7] = 0.0;
1984   devFunValue[8] = 4.0*gc[2]-1.0;
1985   
1986   devFunValue[9] = 4.0*gc[0]-1.0;
1987   devFunValue[10] = 0.0;
1988   devFunValue[11] = 0.0;
1989   
1990   devFunValue[12] = -4.0*gc[1];
1991   devFunValue[13] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[1]);;
1992   devFunValue[14] = -4.0*gc[1];
1993   
1994   devFunValue[15] = -4.0*gc[2];
1995   devFunValue[16] = -4.0*gc[2];
1996   devFunValue[17] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[2]);
1997   
1998   devFunValue[18] = 0.0;
1999   devFunValue[19] = 4.0*gc[2];
2000   devFunValue[20] = 4.0*gc[1];
2001   
2002   devFunValue[21] = 4.0*gc[1];
2003   devFunValue[22] = 4.0*gc[0];
2004   devFunValue[23] = 0.0;
2005   
2006   devFunValue[24] = 4.0*((1-gc[0]-gc[1]-gc[2])-gc[0]);
2007   devFunValue[25] = -4.0*gc[0];
2008   devFunValue[26] = -4.0*gc[0];
2009   
2010   devFunValue[27] = 4.0*gc[2];
2011   devFunValue[28] = 0.0;
2012   devFunValue[29] = 4.0*gc[0];
2013   DEV_SHAPE_FUN_MACRO_END;
2014 }
2015
2016 /*!
2017  * Init Pyramid Reference coordinates ans Shape function.
2018  * Case A.
2019  */
2020 void GaussInfo::pyra5aInit() 
2021 {
2022   LOCAL_COORD_MACRO_BEGIN;
2023   case 0:
2024     coords[0] = PYRA5A_REF[0];
2025     coords[1] = PYRA5A_REF[1];
2026     coords[2] = PYRA5A_REF[2];
2027     break;
2028   case 1:
2029     coords[0] = PYRA5A_REF[3];
2030     coords[1] = PYRA5A_REF[4];
2031     coords[2] = PYRA5A_REF[5];
2032     break;
2033   case 2:
2034     coords[0] = PYRA5A_REF[6];
2035     coords[1] = PYRA5A_REF[7];
2036     coords[2] = PYRA5A_REF[8];
2037     break;
2038   case 3:
2039     coords[0] = PYRA5A_REF[9];
2040     coords[1] = PYRA5A_REF[10];
2041     coords[2] = PYRA5A_REF[11];
2042     break;
2043   case 4:
2044     coords[0] = PYRA5A_REF[12];
2045     coords[1] = PYRA5A_REF[13];
2046     coords[2] = PYRA5A_REF[14];
2047     break;
2048   LOCAL_COORD_MACRO_END;
2049   
2050   SHAPE_FUN_MACRO_BEGIN;
2051   funValue[0] = 0.25*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2052   funValue[1] = 0.25*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2053   funValue[2] = 0.25*(gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2054   funValue[3] = 0.25*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2055   funValue[4] = gc[2];
2056   SHAPE_FUN_MACRO_END;
2057
2058   DEV_SHAPE_FUN_MACRO_BEGIN;
2059
2060   devFunValue[0] = (-(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2061   devFunValue[1] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2062   devFunValue[2] = ((-gc[0]+gc[1]+gc[2]-1.0)+(-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]))/((1.0-gc[2])*4.0);
2063   
2064   devFunValue[3] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2065   devFunValue[4] = (-(-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2066   devFunValue[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)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]))/((1.0-gc[2])*4.0);
2067
2068   devFunValue[6] = ((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2069   devFunValue[7] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2070   devFunValue[8] = ((gc[0]-gc[1]+gc[2]-1.0)+(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]))/((1.0-gc[2])*4.0);
2071   
2072   devFunValue[9] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2073   devFunValue[10] = ((gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2074   devFunValue[11] = ((gc[0]+gc[1]+gc[2]-1.0)+(-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]))/((1.0-gc[2])*4.0);
2075   
2076   devFunValue[12] = 0.0;
2077   devFunValue[13] = 0.0;
2078   devFunValue[14] = 1.0;
2079
2080   DEV_SHAPE_FUN_MACRO_END;
2081 }
2082
2083 /*!
2084  * Init Pyramid Reference coordinates ans Shape function.
2085  * Case B.
2086  */
2087 void GaussInfo::pyra5bInit() 
2088 {
2089   LOCAL_COORD_MACRO_BEGIN;
2090   case 0:
2091     coords[0] = PYRA5B_REF[0];
2092     coords[1] = PYRA5B_REF[1];
2093     coords[2] = PYRA5B_REF[2];
2094     break;
2095   case 1:
2096     coords[0] = PYRA5B_REF[3];
2097     coords[1] = PYRA5B_REF[4];
2098     coords[2] = PYRA5B_REF[5];
2099     break;
2100   case 2:
2101     coords[0] = PYRA5B_REF[6];
2102     coords[1] = PYRA5B_REF[7];
2103     coords[2] = PYRA5B_REF[8];
2104     break;
2105   case 3:
2106     coords[0] = PYRA5B_REF[9];
2107     coords[1] = PYRA5B_REF[10];
2108     coords[2] = PYRA5B_REF[11];
2109     break;
2110   case 4:
2111     coords[0] = PYRA5B_REF[12];
2112     coords[1] = PYRA5B_REF[13];
2113     coords[2] = PYRA5B_REF[14];
2114     break;
2115   LOCAL_COORD_MACRO_END;
2116   
2117   SHAPE_FUN_MACRO_BEGIN;
2118   funValue[0] = 0.25*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2119   funValue[3] = 0.25*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2120   funValue[2] = 0.25*(gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2121   funValue[1] = 0.25*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2122   funValue[4] = gc[2];
2123   SHAPE_FUN_MACRO_END;
2124   
2125   DEV_SHAPE_FUN_MACRO_BEGIN;
2126   devFunValue[0] = (-(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2127   devFunValue[1] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2128   devFunValue[2] = ((-gc[0]+gc[1]+gc[2]-1.0)+(-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]))/((1.0-gc[2])*4.0);
2129
2130   devFunValue[3] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2131   devFunValue[4] = ((gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2132   devFunValue[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)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]))/((1.0-gc[2])*4.0);
2133   
2134   devFunValue[6] = ((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2135   devFunValue[7] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2136   devFunValue[8] = ((gc[0]-gc[1]+gc[2]-1.0)+(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]))/((1.0-gc[2])*4.0);
2137   
2138   devFunValue[9] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2139   devFunValue[10] = (-(-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))/((1.0-gc[2])*4.0);
2140   devFunValue[11] = ((-gc[0]-gc[1]+gc[2]-1.0)+(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]))/((1.0-gc[2])*4.0);
2141   
2142   devFunValue[12] = 0.0;
2143   devFunValue[13] = 0.0;
2144   devFunValue[14] = 1.0;
2145   DEV_SHAPE_FUN_MACRO_END;
2146 }
2147
2148 /*!
2149  * Init Quadratic Pyramid Reference coordinates ans Shape function.
2150  * Case A.
2151  */
2152 void GaussInfo::pyra13aInit() 
2153 {
2154   LOCAL_COORD_MACRO_BEGIN;
2155   case 0:
2156     coords[0] = PYRA13A_REF[0];
2157     coords[1] = PYRA13A_REF[1];
2158     coords[2] = PYRA13A_REF[2];
2159     break;
2160   case 1:
2161     coords[0] = PYRA13A_REF[3];
2162     coords[1] = PYRA13A_REF[4];
2163     coords[2] = PYRA13A_REF[5];
2164     break;
2165   case 2:
2166     coords[0] = PYRA13A_REF[6];
2167     coords[1] = PYRA13A_REF[7];
2168     coords[2] = PYRA13A_REF[8];
2169     break;
2170   case 3:
2171     coords[0] = PYRA13A_REF[9];
2172     coords[1] = PYRA13A_REF[10];
2173     coords[2] = PYRA13A_REF[11];
2174     break;
2175   case 4:
2176     coords[0] = PYRA13A_REF[12];
2177     coords[1] = PYRA13A_REF[13];
2178     coords[2] = PYRA13A_REF[14];
2179     break;
2180   case 5:
2181     coords[0] = PYRA13A_REF[15];
2182     coords[1] = PYRA13A_REF[16];
2183     coords[2] = PYRA13A_REF[17];
2184     break;
2185   case 6:
2186     coords[0] = PYRA13A_REF[18];
2187     coords[1] = PYRA13A_REF[19];
2188     coords[2] = PYRA13A_REF[20];
2189     break;
2190   case 7:
2191     coords[0] = PYRA13A_REF[21];
2192     coords[1] = PYRA13A_REF[22];
2193     coords[2] = PYRA13A_REF[23];
2194     break;
2195   case 8:
2196     coords[0] = PYRA13A_REF[24];
2197     coords[1] = PYRA13A_REF[25];
2198     coords[2] = PYRA13A_REF[26];
2199     break;
2200   case 9:
2201     coords[0] = PYRA13A_REF[27];
2202     coords[1] = PYRA13A_REF[28];
2203     coords[2] = PYRA13A_REF[29];
2204     break;
2205   case 10:
2206     coords[0] = PYRA13A_REF[30];
2207     coords[1] = PYRA13A_REF[31];
2208     coords[2] = PYRA13A_REF[32];
2209     break;
2210   case 11:
2211     coords[0] = PYRA13A_REF[33];
2212     coords[1] = PYRA13A_REF[34];
2213     coords[2] = PYRA13A_REF[35];
2214     break;
2215   case 12:
2216     coords[0] = PYRA13A_REF[36];
2217     coords[1] = PYRA13A_REF[37];
2218     coords[2] = PYRA13A_REF[38];
2219     break;
2220   LOCAL_COORD_MACRO_END;
2221   
2222   SHAPE_FUN_MACRO_BEGIN;
2223   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]);
2224   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]);
2225   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]);
2226   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]);
2227   
2228   funValue[4]=2.0*gc[2]*(gc[2]-0.5);
2229   
2230   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]);
2231   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]);
2232   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]);
2233   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]);
2234   
2235   funValue[9]=gc[2]*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2236   funValue[10]=gc[2]*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2237   funValue[11]=gc[2]*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2238   funValue[12]=gc[2]*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2239   SHAPE_FUN_MACRO_END;
2240
2241   DEV_SHAPE_FUN_MACRO_BEGIN;
2242   devFunValue[0] = ((-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-((-gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0))*(gc[0]-0.5))/(2.0*(1.0-gc[2]));
2243   devFunValue[1] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*(gc[0]-0.5)/(2.0*(1.0-gc[2]));
2244   devFunValue[2] = ((-gc[0]+gc[1]+gc[2]-1.0)+(-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]))*(gc[0]-0.5)/(2.0*(1.0-gc[2]));
2245   
2246   devFunValue[3] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))*(gc[1]-0.5)/(2.0*(1.0-gc[2]));
2247   devFunValue[4] = ((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)-((-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*(gc[1]-0.5))/(2.0*(1.0-gc[2]));
2248   devFunValue[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)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]))*(gc[1]-0.5)/(2.0*(1.0-gc[2]));
2249   
2250   devFunValue[6] = (((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*(-gc[0]-0.5)-(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2251   devFunValue[7] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*(-gc[0]-0.5)/(2.0*(1.0-gc[2]));
2252   devFunValue[8] = ((gc[0]+gc[1]+gc[2]-1.0)+(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]))*(-gc[0]-0.5)/(2.0*(1.0-gc[2]));
2253   
2254   devFunValue[9] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*(-gc[1]-0.5)/(2.0*(1.0-gc[2]));
2255   devFunValue[10] = (((-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0))*(-gc[1]-0.5)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2256   devFunValue[11] = ((-gc[0]+gc[1]+gc[2]-1.0)+(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]))*(-gc[1]-0.5)/(2.0*(1.0-gc[2]));
2257   
2258   
2259   devFunValue[12] = 0.0;
2260   devFunValue[13] = 0.0;
2261   devFunValue[14] = 4.0*gc[2]-1.0;
2262   
2263   devFunValue[15] = -((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2264   devFunValue[16] = -((-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2265   devFunValue[17] = ((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(-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]))/(2.0*(1.0-gc[2]));
2266   
2267   devFunValue[18] = -((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2268   devFunValue[19] = -((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2269   devFunValue[20] = ((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-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]))/(2.0*(1.0-gc[2]));
2270   
2271   devFunValue[21] = -((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2272   devFunValue[22] = -((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2273   devFunValue[23] = ((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(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]))/(2.0*(1.0-gc[2]));
2274   
2275   devFunValue[24] = -((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2276   devFunValue[25] = -((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2277   devFunValue[26] = ((-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(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]))/(2.0*(1.0-gc[2]));
2278   
2279   
2280   devFunValue[27] =(-(-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2281   devFunValue[28] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2282   devFunValue[29] = (-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((-gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2283
2284   devFunValue[30] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2285   devFunValue[31] = (-(gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2286   devFunValue[32] = (-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2287   
2288   devFunValue[33] = ((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2289   devFunValue[34] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2290   devFunValue[35] = (gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2291   
2292   devFunValue[36] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2293   devFunValue[37] = ((-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2294   devFunValue[38] = (gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2295   DEV_SHAPE_FUN_MACRO_END;
2296 }
2297
2298 /*!
2299  * Init Quadratic Pyramid Reference coordinates ans Shape function.
2300  * Case B.
2301  */
2302 void GaussInfo::pyra13bInit() 
2303 {
2304   LOCAL_COORD_MACRO_BEGIN;
2305   case 0:
2306     coords[0] = PYRA13B_REF[0];
2307     coords[1] = PYRA13B_REF[1];
2308     coords[2] = PYRA13B_REF[2];
2309     break;
2310   case 1:
2311     coords[0] = PYRA13B_REF[3];
2312     coords[1] = PYRA13B_REF[4];
2313     coords[2] = PYRA13B_REF[5];
2314     break;
2315   case 2:
2316     coords[0] = PYRA13B_REF[6];
2317     coords[1] = PYRA13B_REF[7];
2318     coords[2] = PYRA13B_REF[8];
2319     break;
2320   case 3:
2321     coords[0] = PYRA13B_REF[9];
2322     coords[1] = PYRA13B_REF[10];
2323     coords[2] = PYRA13B_REF[11];
2324     break;
2325   case 4:
2326     coords[0] = PYRA13B_REF[12];
2327     coords[1] = PYRA13B_REF[13];
2328     coords[2] = PYRA13B_REF[14];
2329     break;
2330   case 5:
2331     coords[0] = PYRA13B_REF[15];
2332     coords[1] = PYRA13B_REF[16];
2333     coords[2] = PYRA13B_REF[17];
2334     break;
2335   case 6:
2336     coords[0] = PYRA13B_REF[18];
2337     coords[1] = PYRA13B_REF[19];
2338     coords[2] = PYRA13B_REF[20];
2339     break;
2340   case 7:
2341     coords[0] = PYRA13B_REF[21];
2342     coords[1] = PYRA13B_REF[22];
2343     coords[2] = PYRA13B_REF[23];
2344     break;
2345   case 8:
2346     coords[0] = PYRA13B_REF[24];
2347     coords[1] = PYRA13B_REF[25];
2348     coords[2] = PYRA13B_REF[26];
2349     break;
2350   case 9:
2351     coords[0] = PYRA13B_REF[27];
2352     coords[1] = PYRA13B_REF[28];
2353     coords[2] = PYRA13B_REF[29];
2354     break;
2355   case 10:
2356     coords[0] = PYRA13B_REF[30];
2357     coords[1] = PYRA13B_REF[31];
2358     coords[2] = PYRA13B_REF[32];
2359     break;
2360   case 11:
2361     coords[0] = PYRA13B_REF[33];
2362     coords[1] = PYRA13B_REF[34];
2363     coords[2] = PYRA13B_REF[35];
2364     break;
2365   case 12:
2366     coords[0] = PYRA13B_REF[36];
2367     coords[1] = PYRA13B_REF[37];
2368     coords[2] = PYRA13B_REF[38];
2369     break;
2370   LOCAL_COORD_MACRO_END;
2371   
2372   SHAPE_FUN_MACRO_BEGIN;
2373   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]);
2374   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]);
2375   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]);
2376   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]);
2377   
2378   funValue[4] =2.0*gc[2]*(gc[2]-0.5);
2379   
2380   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]);
2381   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]);
2382   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]);
2383   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]);
2384   
2385   funValue[9] =gc[2]*(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2386   funValue[10]=gc[2]*(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2387   funValue[11]=gc[2]*(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2388   funValue[12]=gc[2]*(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2]);
2389   
2390   SHAPE_FUN_MACRO_END;
2391   
2392   DEV_SHAPE_FUN_MACRO_BEGIN;
2393   devFunValue[0] = ((-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-((-gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0))*(gc[0]-0.5))/(2.0*(1.0-gc[2]));
2394   devFunValue[1] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*(gc[0]-0.5)/(2.0*(1.0-gc[2]));
2395   devFunValue[2] = ((-gc[0]+gc[1]+gc[2]-1.0)+(-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]))*(gc[0]-0.5)/(2.0*(1.0-gc[2]));
2396   
2397   devFunValue[9] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))*(gc[1]-0.5)/(2.0*(1.0-gc[2]));
2398   devFunValue[10] = ((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)-((-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*(gc[1]-0.5))/(2.0*(1.0-gc[2]));
2399   devFunValue[11] = ((-gc[0]-gc[1]+gc[2]-1.0)+(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]))*(gc[1]-0.5)/(2.0*(1.0-gc[2]));
2400   
2401   devFunValue[6] = (((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*(-gc[0]-0.5)-(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2402   devFunValue[7] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*(-gc[0]-0.5)/(2.0*(1.0-gc[2]));
2403   devFunValue[8] = ((gc[0]+gc[1]+gc[2]-1.0)+(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]))*(-gc[0]-0.5)/(2.0*(1.0-gc[2]));
2404   
2405   devFunValue[3] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*(-gc[1]-0.5)/(2.0*(1.0-gc[2]));
2406   devFunValue[4] = (((-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0))*(-gc[1]-0.5)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2407   devFunValue[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)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2]))*(-gc[1]-0.5)/(2.0*(1.0-gc[2]));
2408   
2409   
2410   devFunValue[12] = 0.0;
2411   devFunValue[13] = 0.0;
2412   devFunValue[14] = 4.0*gc[2]-1.0;
2413   
2414   devFunValue[24] = ((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2415   devFunValue[25] = ((-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2416   devFunValue[26] = -((-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(-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]))/(2.0*(1.0-gc[2]));
2417   
2418   devFunValue[21] = ((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2419   devFunValue[22] = ((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2420   devFunValue[23] = -((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)+(-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]))/(2.0*(1.0-gc[2]));
2421   
2422   devFunValue[18] = ((gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2423   devFunValue[19] = ((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2424   devFunValue[20] = -((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)+(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]))/(2.0*(1.0-gc[2]));
2425   
2426   devFunValue[15] = ((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2427   devFunValue[16] = ((gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0))/(2.0*(1.0-gc[2]));
2428   devFunValue[17] = -((-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)+(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]))/(2.0*(1.0-gc[2]));
2429   
2430   
2431   devFunValue[27] =(-(-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2432   devFunValue[28] = ((-gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2433   devFunValue[29] = (-gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((-gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2434
2435   devFunValue[36] = ((-gc[0]-gc[1]+gc[2]-1.0)-(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2436   devFunValue[37] = (-(gc[0]-gc[1]+gc[2]-1.0)-(-gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2437   devFunValue[38] = (-gc[0]-gc[1]+gc[2]-1.0)*(gc[0]-gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]-gc[1]+gc[2]-1.0)+(-gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2438   
2439   devFunValue[33] = ((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2440   devFunValue[34] = ((gc[0]-gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2441   devFunValue[35] = (gc[0]-gc[1]+gc[2]-1.0)*(gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]+gc[1]+gc[2]-1.0)+(gc[0]-gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2442   
2443   devFunValue[30] = ((-gc[0]+gc[1]+gc[2]-1.0)-(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2444   devFunValue[31] = ((-gc[0]+gc[1]+gc[2]-1.0)+(gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2445   devFunValue[32] = (gc[0]+gc[1]+gc[2]-1.0)*(-gc[0]+gc[1]+gc[2]-1.0)/(1.0-gc[2])/(1.0-gc[2])+((gc[0]+gc[1]+gc[2]-1.0)+(-gc[0]+gc[1]+gc[2]-1.0))*gc[2]/(1.0-gc[2]);
2446   DEV_SHAPE_FUN_MACRO_END;
2447 }
2448
2449
2450 /*!
2451  * Init Pentahedron Reference coordinates and Shape function.
2452  * Case A.
2453  */
2454 void GaussInfo::penta6aInit() 
2455 {
2456   LOCAL_COORD_MACRO_BEGIN;
2457   case 0:
2458     coords[0] = PENTA6A_REF[0];
2459     coords[1] = PENTA6A_REF[1];
2460     coords[2] = PENTA6A_REF[2];
2461     break;
2462   case 1:
2463     coords[0] = PENTA6A_REF[3];
2464     coords[1] = PENTA6A_REF[4];
2465     coords[2] = PENTA6A_REF[5];
2466     break;
2467   case 2:
2468     coords[0] = PENTA6A_REF[6];
2469     coords[1] = PENTA6A_REF[7];
2470     coords[2] = PENTA6A_REF[8];
2471     break;
2472   case 3:
2473     coords[0] = PENTA6A_REF[9];
2474     coords[1] = PENTA6A_REF[10];
2475     coords[2] = PENTA6A_REF[11];
2476     break;
2477   case 4:
2478     coords[0] = PENTA6A_REF[12];
2479     coords[1] = PENTA6A_REF[13];
2480     coords[2] = PENTA6A_REF[14];
2481     break;
2482   case 5:
2483     coords[0] = PENTA6A_REF[15];
2484     coords[1] = PENTA6A_REF[16];
2485     coords[2] = PENTA6A_REF[17];
2486     break;
2487   LOCAL_COORD_MACRO_END;
2488   
2489   SHAPE_FUN_MACRO_BEGIN;
2490   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
2491   funValue[1] = 0.5*gc[2]*(1.0 - gc[0]);
2492   funValue[2] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2493   
2494   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
2495   funValue[4] = 0.5*gc[2]*(gc[0] + 1.0);
2496   funValue[5] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2497   SHAPE_FUN_MACRO_END;
2498
2499   DEV_SHAPE_FUN_MACRO_BEGIN;
2500
2501   devFunValue[0] = 0.5*gc[1]*(-1.0);
2502   devFunValue[1] = 0.5*(1.0 - gc[0]);
2503   devFunValue[2] = 0.0;
2504
2505   devFunValue[3] = 0.5*gc[2]*(-1.0);
2506   devFunValue[4] = 0.0;
2507   devFunValue[5] = 0.5*(1.0 - gc[0]);
2508
2509   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(-1.0);
2510   devFunValue[7] = 0.5*(-1.0)*(1.0 - gc[0]);
2511   devFunValue[8] = 0.5*(-1.0)*(1.0 - gc[0]);
2512
2513   devFunValue[9] = 0.5*gc[1];
2514   devFunValue[10] = 0.5*(gc[0] + 1.0);
2515   devFunValue[11] = 0.0;
2516
2517   devFunValue[12] = 0.5*gc[2];
2518   devFunValue[13] = 0.0;
2519   devFunValue[14] = 0.5*(gc[0] + 1.0);
2520
2521   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2]);
2522   devFunValue[16] = 0.5*(-1.0)*(1.0 + gc[0]);
2523   devFunValue[17] = 0.5*(-1.0)*(1.0 + gc[0]);
2524
2525   DEV_SHAPE_FUN_MACRO_END;
2526 }
2527
2528 /*!
2529  * Init Pentahedron Reference coordinates and Shape function.
2530  * Case B.
2531  */
2532 void GaussInfo::penta6bInit() 
2533 {
2534   LOCAL_COORD_MACRO_BEGIN;
2535   case 0:
2536     coords[0] = PENTA6B_REF[0];
2537     coords[1] = PENTA6B_REF[1];
2538     coords[2] = PENTA6B_REF[2];
2539     break;
2540   case 1:
2541     coords[0] = PENTA6B_REF[3];
2542     coords[1] = PENTA6B_REF[4];
2543     coords[2] = PENTA6B_REF[5];
2544     break;
2545   case 2:
2546     coords[0] = PENTA6B_REF[6];
2547     coords[1] = PENTA6B_REF[7];
2548     coords[2] = PENTA6B_REF[8];
2549     break;
2550   case 3:
2551     coords[0] = PENTA6B_REF[9];
2552     coords[1] = PENTA6B_REF[10];
2553     coords[2] = PENTA6B_REF[11];
2554     break;
2555   case 4:
2556     coords[0] = PENTA6B_REF[12];
2557     coords[1] = PENTA6B_REF[13];
2558     coords[2] = PENTA6B_REF[14];
2559     break;
2560   case 5:
2561     coords[0] = PENTA6B_REF[15];
2562     coords[1] = PENTA6B_REF[16];
2563     coords[2] = PENTA6B_REF[17];
2564     break;
2565   LOCAL_COORD_MACRO_END;
2566   
2567   SHAPE_FUN_MACRO_BEGIN;
2568   funValue[0] = 0.5*gc[1]*(1.0 - gc[0]);
2569   funValue[2] = 0.5*gc[2]*(1.0 - gc[0]);
2570   funValue[1] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2571   funValue[3] = 0.5*gc[1]*(gc[0] + 1.0);
2572   funValue[5] = 0.5*gc[2]*(gc[0] + 1.0);
2573   funValue[4] = 0.5*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2574   SHAPE_FUN_MACRO_END;
2575   
2576   DEV_SHAPE_FUN_MACRO_BEGIN;
2577   devFunValue[0] = 0.5*gc[1]*(-1.0);
2578   devFunValue[1] = 0.5*(1.0 - gc[0]);
2579   devFunValue[2] = 0.0;
2580
2581   devFunValue[6] = 0.5*gc[2]*(-1.0);
2582   devFunValue[7] = 0.0;
2583   devFunValue[8] = 0.5*(1.0 - gc[0]);
2584
2585   devFunValue[3] = 0.5*(1.0 - gc[1] - gc[2])*(-1.0);
2586   devFunValue[4] = 0.5*(-1.0)*(1.0 - gc[0]);
2587   devFunValue[5] = 0.5*(-1.0)*(1.0 - gc[0]);
2588
2589   devFunValue[9] = 0.5*gc[1];
2590   devFunValue[10] = 0.5*(gc[0] + 1.0);
2591   devFunValue[11] = 0.0;
2592
2593   devFunValue[15] = 0.5*gc[2];
2594   devFunValue[16] = 0.0;
2595   devFunValue[17] = 0.5*(gc[0] + 1.0);
2596
2597   devFunValue[12] = 0.5*(1.0 - gc[1] - gc[2]);
2598   devFunValue[13] = 0.5*(-1.0)*(1.0 + gc[0]);
2599   devFunValue[14] = 0.5*(-1.0)*(1.0 + gc[0]);
2600   DEV_SHAPE_FUN_MACRO_END;
2601 }
2602
2603 /*!
2604  * This shapefunc map is same as degenerated tria3aInit
2605  */
2606 void GaussInfo::penta6DegTria3aInit() 
2607 {
2608   LOCAL_COORD_MACRO_BEGIN;
2609  case  0:
2610    coords[0] = -1.0;
2611    coords[1] =  1.0;
2612    coords[2] =  0.0;
2613    break;
2614  case  1:
2615    coords[0] = -1.0;
2616    coords[1] = -1.0;
2617    coords[2] =  0.0;
2618    break;
2619  case  2:
2620    coords[0] =  1.0;
2621    coords[1] = -1.0;
2622    coords[2] =  0.0;
2623    break;
2624  case  3:
2625    coords[0] =  0.0;
2626    coords[1] =  0.0;
2627    coords[2] =  0.0;
2628    break;
2629  case  4:
2630    coords[0] =  0.0;
2631    coords[1] =  0.0;
2632    coords[2] =  0.0;
2633    break;
2634  case  5:
2635    coords[0] =  0.0;
2636    coords[1] =  0.0;
2637    coords[2] =  0.0;
2638    break;
2639    LOCAL_COORD_MACRO_END;
2640
2641    SHAPE_FUN_MACRO_BEGIN;
2642    funValue[0] = 0.5*(1.0 + gc[1]);
2643    funValue[1] = -0.5*(gc[0] + gc[1]);
2644    funValue[2] = 0.5*(1.0 + gc[0]);
2645    funValue[3] = 0.;
2646    funValue[4] = 0.;
2647    funValue[5] = 0.;
2648    SHAPE_FUN_MACRO_END;
2649    
2650   DEV_SHAPE_FUN_MACRO_BEGIN;
2651   devFunValue[0] = 0.0;
2652   devFunValue[1] = 0.5;
2653   devFunValue[2] = 0.0;
2654   
2655   devFunValue[3] = -0.5;
2656   devFunValue[4] = -0.5;
2657   devFunValue[5] = 0.0;
2658   
2659   devFunValue[6] = 0.5;
2660   devFunValue[7] = 0.0;
2661   devFunValue[8] = 0.0;
2662   
2663   devFunValue[9] = 0.0;
2664   devFunValue[10] = 0.0;
2665   devFunValue[11] = 0.0;
2666   
2667   devFunValue[12] = 0.0;
2668   devFunValue[13] = 0.0;
2669   devFunValue[14] = 0.0;
2670   
2671   devFunValue[15] = 0.0;
2672   devFunValue[16] = 0.0;
2673   devFunValue[17] = 0.0;
2674   DEV_SHAPE_FUN_MACRO_END;
2675 }
2676
2677 /*!
2678  * This shapefunc map is same as degenerated tria3bInit
2679  */
2680 void GaussInfo::penta6DegTria3bInit() 
2681 {
2682   LOCAL_COORD_MACRO_BEGIN;
2683  case  0:
2684    coords[0] =  0.0;
2685    coords[1] =  0.0;
2686    coords[2] =  0.0;
2687    break;
2688  case  1:
2689    coords[0] =  1.0;
2690    coords[1] =  0.0;
2691    coords[2] =  0.0;
2692    break;
2693  case  2:
2694    coords[0] =  0.0;
2695    coords[1] =  1.0;
2696    coords[2] =  0.0;
2697    break;
2698  case  3:
2699    coords[0] =  0.0;
2700    coords[1] =  0.0;
2701    coords[2] =  0.0;
2702    break;
2703  case  4:
2704    coords[0] =  0.0;
2705    coords[1] =  0.0;
2706    coords[2] =  0.0;
2707    break;
2708  case  5:
2709    coords[0] =  0.0;
2710    coords[1] =  0.0;
2711    coords[2] =  0.0;
2712    break;
2713    LOCAL_COORD_MACRO_END;
2714
2715    SHAPE_FUN_MACRO_BEGIN;
2716    funValue[0] = 1.0 - gc[0] - gc[1];
2717    funValue[1] = gc[0];
2718    funValue[2] = gc[1];
2719    funValue[3] = 0.;
2720    funValue[4] = 0.;
2721    funValue[5] = 0.;
2722    SHAPE_FUN_MACRO_END;
2723    
2724   DEV_SHAPE_FUN_MACRO_BEGIN;
2725   devFunValue[0] = -1.0;
2726   devFunValue[1] = -1.0;
2727   devFunValue[2] = 0.0;
2728   
2729   devFunValue[3] = 1.0;
2730   devFunValue[4] = 0.0;
2731   devFunValue[5] = 0.0;
2732   
2733   devFunValue[6] = 0.0;
2734   devFunValue[7] = 1.0;
2735   devFunValue[8] = 0.0;
2736   
2737   devFunValue[9] = 0.0;
2738   devFunValue[10] = 0.0;
2739   devFunValue[11] = 0.0;
2740   
2741   devFunValue[12] = 0.0;
2742   devFunValue[13] = 0.0;
2743   devFunValue[14] = 0.0;
2744   
2745   devFunValue[15] = 0.0;
2746   devFunValue[16] = 0.0;
2747   devFunValue[17] = 0.0;
2748   DEV_SHAPE_FUN_MACRO_END;
2749 }
2750
2751 /*!
2752  * Init Pentahedron Reference coordinates and Shape function.
2753  * Case A.
2754  */
2755 void GaussInfo::penta15aInit() 
2756 {
2757   LOCAL_COORD_MACRO_BEGIN;
2758   case 0:
2759     coords[0] = PENTA15A_REF[0];
2760     coords[1] = PENTA15A_REF[1];
2761     coords[2] = PENTA15A_REF[2];
2762     break;
2763   case 1:
2764     coords[0] = PENTA15A_REF[3];
2765     coords[1] = PENTA15A_REF[4];
2766     coords[2] = PENTA15A_REF[5];
2767     break;
2768   case 2:
2769     coords[0] = PENTA15A_REF[6];
2770     coords[1] = PENTA15A_REF[7];
2771     coords[2] = PENTA15A_REF[8];
2772     break;
2773   case 3:
2774     coords[0] = PENTA15A_REF[9];
2775     coords[1] = PENTA15A_REF[10];
2776     coords[2] = PENTA15A_REF[11];
2777     break;
2778   case 4:
2779     coords[0] = PENTA15A_REF[12];
2780     coords[1] = PENTA15A_REF[13];
2781     coords[2] = PENTA15A_REF[14];
2782     break;
2783   case 5:
2784     coords[0] = PENTA15A_REF[15];
2785     coords[1] = PENTA15A_REF[16];
2786     coords[2] = PENTA15A_REF[17];
2787     break;
2788   case 6:
2789     coords[0] = PENTA15A_REF[18];
2790     coords[1] = PENTA15A_REF[19];
2791     coords[2] = PENTA15A_REF[20];
2792     break;
2793   case 7:
2794     coords[0] = PENTA15A_REF[21];
2795     coords[1] = PENTA15A_REF[22];
2796     coords[2] = PENTA15A_REF[23];
2797     break;
2798   case 8:
2799     coords[0] = PENTA15A_REF[24];
2800     coords[1] = PENTA15A_REF[25];
2801     coords[2] = PENTA15A_REF[26];
2802     break;
2803   case 9:
2804     coords[0] = PENTA15A_REF[27];
2805     coords[1] = PENTA15A_REF[28];
2806     coords[2] = PENTA15A_REF[29];
2807     break;
2808   case 10:
2809     coords[0] = PENTA15A_REF[30];
2810     coords[1] = PENTA15A_REF[31];
2811     coords[2] = PENTA15A_REF[32];
2812     break;
2813   case 11:
2814     coords[0] = PENTA15A_REF[33];
2815     coords[1] = PENTA15A_REF[34];
2816     coords[2] = PENTA15A_REF[35];
2817     break;
2818   case 12:
2819     coords[0] = PENTA15A_REF[36];
2820     coords[1] = PENTA15A_REF[37];
2821     coords[2] = PENTA15A_REF[38];
2822     break;
2823   case 13:
2824     coords[0] = PENTA15A_REF[39];
2825     coords[1] = PENTA15A_REF[40];
2826     coords[2] = PENTA15A_REF[41];
2827     break;
2828   case 14:
2829     coords[0] = PENTA15A_REF[42];
2830     coords[1] = PENTA15A_REF[43];
2831     coords[2] = PENTA15A_REF[44];
2832     break;
2833   LOCAL_COORD_MACRO_END;
2834   
2835   SHAPE_FUN_MACRO_BEGIN;
2836   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
2837   funValue[1] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
2838   funValue[2] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2839   
2840   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
2841   funValue[4] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
2842   funValue[5] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
2843   
2844   funValue[6] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
2845   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2846   funValue[8] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
2847   
2848   funValue[9] = gc[1]*(1.0 - gc[0]*gc[0]);
2849   funValue[10] = gc[2]*(1.0 - gc[0]*gc[0]);
2850   funValue[11] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
2851   
2852   funValue[12] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
2853   funValue[13] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2854   funValue[14] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
2855   SHAPE_FUN_MACRO_END;
2856   
2857   DEV_SHAPE_FUN_MACRO_BEGIN;
2858
2859   devFunValue[0] = 0.5*gc[1]*(2 * gc[0] - 2 * gc[1] + 1.0 );
2860   devFunValue[1] = (1.0 - gc[0])*(2.0*gc[1] -1 - 0.5*gc[0]);
2861   devFunValue[2] = 0.0;
2862
2863   devFunValue[3] = 0.5*gc[2]*(2 * gc[0] - 2 * gc[2] + 1.0 );
2864   devFunValue[4] = 0.0;
2865   devFunValue[5] = (1.0 - gc[0])*(2.0*gc[2] -1 - 0.5*gc[0]);
2866
2867   devFunValue[6] = 0.5*(1.0 - gc[1] - gc[2])*(2*gc[0] -1.0 + 2*gc[1] + 2*gc[2]);
2868   devFunValue[7] = 0.5*(gc[0] - 1.0)* ( -4*gc[1] -gc[0] -4*gc[2] + 2.0);
2869   devFunValue[8] = 0.5*(gc[0] - 1.0)* ( -4*gc[2] -4*gc[1] -gc[0] + 2.0);
2870
2871   devFunValue[9] = 0.5*gc[1]*( 2*gc[0] + 2*gc[1] -1.0 );
2872   devFunValue[10] = (1.0 + gc[0])*(2.0*gc[1] - 1.0 + 0.5*gc[0]);
2873   devFunValue[11] = 0.0;
2874
2875   devFunValue[12] = 0.5*gc[2]*( 2*gc[0] + 2*gc[2] -1.0 );
2876   devFunValue[13] = 0.0;
2877   devFunValue[14] = (1.0 + gc[0])*(2.0*gc[2] - 1.0 + 0.5*gc[0]);
2878
2879   devFunValue[15] = 0.5*(1.0 - gc[1] - gc[2])*(2 * gc[0] - 2.0 * gc[1] - 2.0 * gc[2] + 1.0 );
2880   devFunValue[16] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[1] + gc[0] - 4.0 * gc[2] + 2.0) ;
2881   devFunValue[17] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[2] - 4*gc[1] + gc[0] + 2.0 );
2882
2883   devFunValue[18] = - 2.0*gc[1]*gc[2];
2884   devFunValue[19] = 2.0*gc[2]*(1.0 - gc[0]);
2885   devFunValue[20] = 2.0*gc[1]*(1.0 - gc[0]);
2886
2887   devFunValue[21] = -2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2888   devFunValue[22] = -2.0*gc[2]*(1.0 - gc[0]);
2889   devFunValue[23] = (-4*gc[2]-2*gc[1]+2.0)*(1.0 - gc[0]);
2890
2891   devFunValue[24] = -2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2892   devFunValue[25] = (-4*gc[1]-2*gc[2]+2)*(1.0 - gc[0]);
2893   devFunValue[26] = -2.0*gc[1]*(1.0 - gc[0]);
2894
2895   devFunValue[27] = gc[1]*(- 2.0 * gc[0]);
2896   devFunValue[28] = (1.0 - gc[0]*gc[0]);
2897   devFunValue[29] = 0.0;
2898
2899   devFunValue[30] = -2.0 * gc[2] *gc[0];
2900   devFunValue[31] = 0.0;
2901   devFunValue[32] = (1.0 - gc[0]*gc[0]);
2902
2903   devFunValue[33] = -2.0 * (1.0 - gc[1] - gc[2]) * gc[0];
2904   devFunValue[34] = -1.0*(1.0 - gc[0]*gc[0]);
2905   devFunValue[35] = -1.0*(1.0 - gc[0]*gc[0]);
2906
2907   devFunValue[36] = 2.0*gc[1]*gc[2];
2908   devFunValue[37] = 2.0*gc[2]*(1.0 + gc[0]);
2909   devFunValue[38] = 2.0*gc[1]*(1.0 + gc[0]);
2910
2911   devFunValue[39] = 2.0*gc[2]*(1.0 - gc[1] - gc[2]);
2912   devFunValue[40] = -2.0*gc[2]*(1.0 + gc[0]);
2913   devFunValue[41] = (2.0 - 2.0 * gc[1] - 4.0 * gc[2])*(1.0 + gc[0]);
2914
2915   devFunValue[42] = 2.0*gc[1]*(1.0 - gc[1] - gc[2]);
2916   devFunValue[43] = (2.0 - 4*gc[1] - 2*gc[2])*(1.0 + gc[0]);
2917   devFunValue[44] = -2.0*gc[1]*(1.0 + gc[0]);
2918
2919   DEV_SHAPE_FUN_MACRO_END;
2920 }
2921
2922 /*!
2923  * Init Qaudratic Pentahedron Reference coordinates and Shape function.
2924  * Case B.
2925  */
2926 void GaussInfo::penta15bInit() 
2927 {
2928   LOCAL_COORD_MACRO_BEGIN;
2929   case 0:
2930     coords[0] = PENTA15B_REF[0];
2931     coords[1] = PENTA15B_REF[1];
2932     coords[2] = PENTA15B_REF[2];
2933     break;
2934   case 1:
2935     coords[0] = PENTA15B_REF[3];
2936     coords[1] = PENTA15B_REF[4];
2937     coords[2] = PENTA15B_REF[5];
2938     break;
2939   case 2:
2940     coords[0] = PENTA15B_REF[6];
2941     coords[1] = PENTA15B_REF[7];
2942     coords[2] = PENTA15B_REF[8];
2943     break;
2944   case 3:
2945     coords[0] = PENTA15B_REF[9];
2946     coords[1] = PENTA15B_REF[10];
2947     coords[2] = PENTA15B_REF[11];
2948     break;
2949   case 4:
2950     coords[0] = PENTA15B_REF[12];
2951     coords[1] = PENTA15B_REF[13];
2952     coords[2] = PENTA15B_REF[14];
2953     break;
2954   case 5:
2955     coords[0] = PENTA15B_REF[15];
2956     coords[1] = PENTA15B_REF[16];
2957     coords[2] = PENTA15B_REF[17];
2958     break;
2959   case 6:
2960     coords[0] = PENTA15B_REF[18];
2961     coords[1] = PENTA15B_REF[19];
2962     coords[2] = PENTA15B_REF[20];
2963     break;
2964   case 7:
2965     coords[0] = PENTA15B_REF[21];
2966     coords[1] = PENTA15B_REF[22];
2967     coords[2] = PENTA15B_REF[23];
2968     break;
2969   case 8:
2970     coords[0] = PENTA15B_REF[24];
2971     coords[1] = PENTA15B_REF[25];
2972     coords[2] = PENTA15B_REF[26];
2973     break;
2974   case 9:
2975     coords[0] = PENTA15B_REF[27];
2976     coords[1] = PENTA15B_REF[28];
2977     coords[2] = PENTA15B_REF[29];
2978     break;
2979   case 10:
2980     coords[0] = PENTA15B_REF[30];
2981     coords[1] = PENTA15B_REF[31];
2982     coords[2] = PENTA15B_REF[32];
2983     break;
2984   case 11:
2985     coords[0] = PENTA15B_REF[33];
2986     coords[1] = PENTA15B_REF[34];
2987     coords[2] = PENTA15B_REF[35];
2988     break;
2989   case 12:
2990     coords[0] = PENTA15B_REF[36];
2991     coords[1] = PENTA15B_REF[37];
2992     coords[2] = PENTA15B_REF[38];
2993     break;
2994   case 13:
2995     coords[0] = PENTA15B_REF[39];
2996     coords[1] = PENTA15B_REF[40];
2997     coords[2] = PENTA15B_REF[41];
2998     break;
2999   case 14:
3000     coords[0] = PENTA15B_REF[42];
3001     coords[1] = PENTA15B_REF[43];
3002     coords[2] = PENTA15B_REF[44];
3003     break;
3004   LOCAL_COORD_MACRO_END;
3005   
3006   SHAPE_FUN_MACRO_BEGIN;
3007   funValue[0] = 0.5*gc[1]*(1.0 - gc[0])*(2.0*gc[1] - 2.0 - gc[0]);
3008   funValue[2] = 0.5*gc[2]*(1.0 - gc[0])*(2.0*gc[2] - 2.0 - gc[0]);
3009   funValue[1] = 0.5*(gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(gc[0] + 2.0*gc[1] + 2.0*gc[2]);
3010   
3011   funValue[3] = 0.5*gc[1]*(1.0 + gc[0])*(2.0*gc[1] - 2.0 + gc[0]);
3012   funValue[5] = 0.5*gc[2]*(1.0 + gc[0])*(2.0*gc[2] - 2.0 + gc[0]);
3013   funValue[4] = 0.5*(-gc[0] - 1.0)*(1.0 - gc[1] - gc[2])*(-gc[0] + 2.0*gc[1] + 2.0*gc[2]);
3014   
3015   funValue[8] = 2.0*gc[1]*gc[2]*(1.0 - gc[0]);
3016   funValue[7] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
3017   funValue[6] = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 - gc[0]);
3018   
3019   funValue[12] = gc[1]*(1.0 - gc[0]*gc[0]);
3020   funValue[14] = gc[2]*(1.0 - gc[0]*gc[0]);
3021   funValue[13] = (1.0 - gc[1] - gc[2])*(1.0 - gc[0]*gc[0]);
3022   
3023   funValue[11] = 2.0*gc[1]*gc[2]*(1.0 + gc[0]);
3024   funValue[10] = 2.0*gc[2]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
3025   funValue[9]  = 2.0*gc[1]*(1.0 - gc[1] - gc[2])*(1.0 + gc[0]);
3026   SHAPE_FUN_MACRO_END;
3027   
3028   DEV_SHAPE_FUN_MACRO_BEGIN;
3029   devFunValue[0] = 0.5*gc[1]*(2 * gc[0] - 2 * gc[1] + 1.0 );
3030   devFunValue[1] = (1.0 - gc[0])*(2.0*gc[1] -1 - 0.5*gc[0]);
3031   devFunValue[2] = 0.0;
3032
3033   devFunValue[6] = 0.5*gc[2]*(2 * gc[0] - 2 * gc[2] + 1.0 );
3034   devFunValue[7] = 0.0;
3035   devFunValue[8] = (1.0 - gc[0])*(2.0*gc[2] -1 - 0.5*gc[0]);
3036
3037   devFunValue[3] = 0.5*(1.0 - gc[1] - gc[2])*(2*gc[0] -1.0 + 2*gc[1] + 2*gc[2]);
3038   devFunValue[4] = 0.5*(gc[0] - 1.0)* ( -4*gc[1] -gc[0] -4*gc[2] + 2.0);
3039   devFunValue[5] = 0.5*(gc[0] - 1.0)* ( -4*gc[2] -4*gc[1] -gc[0] + 2.0);
3040
3041   devFunValue[9] = 0.5*gc[1]*( 2*gc[0] + 2*gc[1] -1.0 );
3042   devFunValue[10] = (1.0 + gc[0])*(2.0*gc[1] - 1.0 + 0.5*gc[0]);
3043   devFunValue[11] = 0.0;
3044
3045   devFunValue[15] = 0.5*gc[2]*( 2*gc[0] + 2*gc[2] -1.0 );
3046   devFunValue[16] = 0.0;
3047   devFunValue[17] = (1.0 + gc[0])*(2.0*gc[2] - 1.0 + 0.5*gc[0]);
3048
3049   devFunValue[12] = 0.5*(1.0 - gc[1] - gc[2])*(2 * gc[0] - 2.0 * gc[1] - 2.0 * gc[2] + 1.0 );
3050   devFunValue[13] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[1] + gc[0] - 4.0 * gc[2] + 2.0) ;
3051   devFunValue[14] = 0.5*(-gc[0] - 1.0) * (-4.0*gc[2] - 4*gc[1] + gc[0] + 2.0 );
3052
3053   devFunValue[24] = - 2.0*gc[1]*gc[2];
3054   devFunValue[25] = 2.0*gc[2]*(1.0 - gc[0]);
3055   devFunValue[26] = 2.0*gc[1]*(1.0 - gc[0]);
3056
3057   devFunValue[21] = -2.0*gc[2]*(1.0 - gc[1] - gc[2]);
3058   devFunValue[22] = -2.0*gc[2]*(1.0 - gc[0]);
3059   devFunValue[23] = (-4*gc[2]-2*gc[1]+2.0)*(1.0 - gc[0]);
3060
3061   devFunValue[18] = -2.0*gc[1]*(1.0 - gc[1] - gc[2]);
3062   devFunValue[19] = (-4*gc[1]-2*gc[2]+2)*(1.0 - gc[0]);
3063   devFunValue[20] = -2.0*gc[1]*(1.0 - gc[0]);
3064
3065
3066   devFunValue[36] = gc[1]*(- 2.0 * gc[0]);
3067   devFunValue[37] = (1.0 - gc[0]*gc[0]);
3068   devFunValue[38] = 0.0;
3069
3070   devFunValue[42] = -2.0 * gc[2] *gc[0];
3071   devFunValue[43] = 0.0;
3072   devFunValue[44] = (1.0 - gc[0]*gc[0]);
3073
3074   devFunValue[39] = -2.0 * (1.0 - gc[1] - gc[2]) * gc[0];
3075   devFunValue[40] = -1.0*(1.0 - gc[0]*gc[0]);
3076   devFunValue[41] = -1.0*(1.0 - gc[0]*gc[0]);
3077
3078   devFunValue[33] = 2.0*gc[1]*gc[2];
3079   devFunValue[34] = 2.0*gc[2]*(1.0 + gc[0]);
3080   devFunValue[35] = 2.0*gc[1]*(1.0 + gc[0]);
3081
3082   devFunValue[30] = 2.0*gc[2]*(1.0 - gc[1] - gc[2]);
3083   devFunValue[31] = -2.0*gc[2]*(1.0 + gc[0]);
3084   devFunValue[32] = (2.0 - 2.0 * gc[1] - 4.0 * gc[2])*(1.0 + gc[0]);
3085
3086   devFunValue[27] = 2.0*gc[1]*(1.0 - gc[1] - gc[2]);
3087   devFunValue[28] = (2.0 - 4*gc[1] - 2*gc[2])*(1.0 + gc[0]);
3088   devFunValue[29] = -2.0*gc[1]*(1.0 + gc[0]);
3089   
3090   DEV_SHAPE_FUN_MACRO_END;
3091 }
3092
3093 void GaussInfo::penta18aInit() 
3094 {
3095   LOCAL_COORD_MACRO_BEGIN;
3096   case 0:
3097     coords[0] = PENTA18A_REF[0];
3098     coords[1] = PENTA18A_REF[1];
3099     coords[2] = PENTA18A_REF[2];
3100     break;
3101   case 1:
3102     coords[0] = PENTA18A_REF[3];
3103     coords[1] = PENTA18A_REF[4];
3104     coords[2] = PENTA18A_REF[5];
3105     break;
3106   case 2:
3107     coords[0] = PENTA18A_REF[6];
3108     coords[1] = PENTA18A_REF[7];
3109     coords[2] = PENTA18A_REF[8];
3110     break;
3111   case 3:
3112     coords[0] = PENTA18A_REF[9];
3113     coords[1] = PENTA18A_REF[10];
3114     coords[2] = PENTA18A_REF[11];
3115     break;
3116   case 4:
3117     coords[0] = PENTA18A_REF[12];
3118     coords[1] = PENTA18A_REF[13];
3119     coords[2] = PENTA18A_REF[14];
3120     break;
3121   case 5:
3122     coords[0] = PENTA18A_REF[15];
3123     coords[1] = PENTA18A_REF[16];
3124     coords[2] = PENTA18A_REF[17];
3125     break;
3126   case 6:
3127     coords[0] = PENTA18A_REF[18];
3128     coords[1] = PENTA18A_REF[19];
3129     coords[2] = PENTA18A_REF[20];
3130     break;
3131   case 7:
3132     coords[0] = PENTA18A_REF[21];
3133     coords[1] = PENTA18A_REF[22];
3134     coords[2] = PENTA18A_REF[23];
3135     break;
3136   case 8:
3137     coords[0] = PENTA18A_REF[24];
3138     coords[1] = PENTA18A_REF[25];
3139     coords[2] = PENTA18A_REF[26];
3140     break;
3141   case 9:
3142     coords[0] = PENTA18A_REF[27];
3143     coords[1] = PENTA18A_REF[28];
3144     coords[2] = PENTA18A_REF[29];
3145     break;
3146   case 10:
3147     coords[0] = PENTA18A_REF[30];
3148     coords[1] = PENTA18A_REF[31];
3149     coords[2] = PENTA18A_REF[32];
3150     break;
3151   case 11:
3152     coords[0] = PENTA18A_REF[33];
3153     coords[1] = PENTA18A_REF[34];
3154     coords[2] = PENTA18A_REF[35];
3155     break;
3156   case 12:
3157     coords[0] = PENTA18A_REF[36];
3158     coords[1] = PENTA18A_REF[37];
3159     coords[2] = PENTA18A_REF[38];
3160     break;
3161   case 13:
3162     coords[0] = PENTA18A_REF[39];
3163     coords[1] = PENTA18A_REF[40];
3164     coords[2] = PENTA18A_REF[41];
3165     break;
3166   case 14:
3167     coords[0] = PENTA18A_REF[42];
3168     coords[1] = PENTA18A_REF[43];
3169     coords[2] = PENTA18A_REF[44];
3170     break;
3171   case 15:
3172     coords[0] = PENTA18A_REF[45];
3173     coords[1] = PENTA18A_REF[46];
3174     coords[2] = PENTA18A_REF[47];
3175     break;
3176   case 16:
3177     coords[0] = PENTA18A_REF[48];
3178     coords[1] = PENTA18A_REF[49];
3179     coords[2] = PENTA18A_REF[50];
3180     break;
3181   case 17:
3182     coords[0] = PENTA18A_REF[51];
3183     coords[1] = PENTA18A_REF[52];
3184     coords[2] = PENTA18A_REF[53];
3185     break;
3186   LOCAL_COORD_MACRO_END;
3187   
3188   SHAPE_FUN_MACRO_BEGIN;
3189   funValue[0] = gc[0]*gc[1]*(gc[0]-1.0)*(2.0*gc[1]-1.0)/2.0;
3190   funValue[1] = gc[0]*gc[2]*(gc[0]-1.0)*(2.0*gc[2]-1.0)/2.0;
3191   funValue[2] = gc[0]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3192   
3193   funValue[3] = gc[0]*gc[1]*(gc[0]+1.0)*(2.0*gc[1]-1.0)/2.0;
3194   funValue[4] = gc[0]*gc[2]*(gc[0]+1.0)*(2.0*gc[2]-1.0)/2.0;
3195   funValue[5] = gc[0]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3196   
3197   funValue[6] = 2.0*gc[0]*gc[1]*gc[2]*(gc[0]-1.0);
3198   funValue[7] = -2.0*gc[0]*gc[2]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3199   funValue[8] = -2.0*gc[0]*gc[1]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3200   
3201   funValue[9] = -gc[1]*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[1]-1.0);
3202   funValue[10] = -gc[2]*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[2]-1.0);
3203   funValue[11] = -(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0);
3204   
3205   funValue[12] = 2.0*gc[0]*gc[1]*gc[2]*(gc[0]+1.0);
3206   funValue[13] = -2.0*gc[0]*gc[2]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3207   funValue[14] = -2.0*gc[0]*gc[1]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3208   
3209   funValue[15] = -4.0*gc[1]*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3210   funValue[16] = 4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3211   funValue[17] = 4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3212   SHAPE_FUN_MACRO_END;
3213   
3214   DEV_SHAPE_FUN_MACRO_BEGIN;
3215   devFunValue[0] = gc[1]*(2.0*gc[0]-1.0)*(2.0*gc[1]-1.0)/2.0;
3216   devFunValue[1] = gc[0]*(gc[0]-1.0)*(4.0*gc[1]-1.0)/2.0;
3217   devFunValue[2] = 0.0;
3218   
3219   devFunValue[3] = gc[2]*(2.0*gc[0]-1.0)*(2.0*gc[2]-1.0)/2.0;
3220   devFunValue[4] = 0.0;
3221   devFunValue[5] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]-1.0)/2.0;
3222   
3223   devFunValue[6] = (2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3224   devFunValue[7] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3225   devFunValue[8] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3226   
3227   devFunValue[9] = gc[1]*(2.0*gc[0]+1.0)*(2.0*gc[1]-1.0)/2.0;
3228   devFunValue[10] = gc[0]*(gc[0]+1.0)*(4.0*gc[1]-1.0)/2.0;
3229   devFunValue[11] = 0.0;
3230   
3231   devFunValue[12] = gc[2]*(2.0*gc[0]+1.0)*(2.0*gc[2]-1.0)/2.0;
3232   devFunValue[13] = 0.0;
3233   devFunValue[14] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]-1.0)/2.0;
3234   
3235   devFunValue[15] = (2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3236   devFunValue[16] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3237   devFunValue[17] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3238   
3239   devFunValue[18] = 2.0*gc[1]*gc[2]*(2.0*gc[0]-1.0);
3240   devFunValue[19] = 2.0*gc[0]*gc[2]*(gc[0]-1.0);
3241   devFunValue[20] = 2.0*gc[0]*gc[1]*(gc[0]-1.0);
3242   
3243   devFunValue[21] = -2.0*gc[2]*(2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3244   devFunValue[22] = -2.0*gc[0]*gc[2]*(gc[0]-1.0);
3245   devFunValue[23] = -2.0*gc[0]*(gc[0]-1.0)*(2.0*gc[2]+gc[1]-1.0);
3246   
3247   devFunValue[24] = -2.0*gc[1]*(2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3248   devFunValue[25] = -2.0*gc[0]*(gc[0]-1.0)*(2.0*gc[1]+gc[2]-1.0);
3249   devFunValue[26] = -2.0*gc[0]*gc[1]*(gc[0]-1.0);
3250   
3251   devFunValue[27] = -2.0*gc[0]*gc[1]*(2.0*gc[1]-1.0);
3252   devFunValue[28] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[1]-1.0);
3253   devFunValue[29] = 0.0;
3254   
3255   devFunValue[30] = -2.0*gc[0]*gc[2]*(2.0*gc[2]-1.0);
3256   devFunValue[31] = 0.0;
3257   devFunValue[32] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]-1.0);
3258   
3259   devFunValue[33] = -2.0*gc[0]*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0);
3260   devFunValue[34] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0);
3261   devFunValue[35] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0);
3262   
3263   devFunValue[36] = 2.0*gc[1]*gc[2]*(2.0*gc[0]+1.0);
3264   devFunValue[37] = 2.0*gc[0]*gc[2]*(gc[0]+1.0);
3265   devFunValue[38] = 2.0*gc[0]*gc[1]*(gc[0]+1.0);
3266   
3267   devFunValue[39] = -2.0*gc[2]*(2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3268   devFunValue[40] = -2.0*gc[0]*gc[2]*(gc[0]+1.0);
3269   devFunValue[41] = -2.0*gc[0]*(gc[0]+1.0)*(2.0*gc[2]+gc[1]-1.0);
3270   
3271   devFunValue[42] = -2.0*gc[1]*(2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3272   devFunValue[43] = -2.0*gc[0]*(gc[0]+1.0)*(2.0*gc[1]+gc[2]-1.0);
3273   devFunValue[44] = -2.0*gc[0]*gc[1]*(gc[0]+1.0);
3274   
3275   devFunValue[45] = -8.0*gc[0]*gc[1]*gc[2];
3276   devFunValue[46] = -4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3277   devFunValue[47] = -4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0);
3278   
3279   devFunValue[48] = 8.0*gc[0]*gc[2]*(gc[2]+gc[1]-1.0);
3280   devFunValue[49] = 4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3281   devFunValue[50] = 4.0*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[2]+gc[1]-1.0);
3282   
3283   devFunValue[51] = 8.0*gc[0]*gc[1]*(gc[2]+gc[1]-1.0);
3284   devFunValue[52] = 4.0*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[1]+gc[2]-1.0);
3285   devFunValue[53] = 4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0);
3286   DEV_SHAPE_FUN_MACRO_END;
3287 }
3288
3289 void GaussInfo::penta18bInit() 
3290 {
3291   LOCAL_COORD_MACRO_BEGIN;
3292   case 0:
3293     coords[0] = PENTA18B_REF[0];
3294     coords[1] = PENTA18B_REF[1];
3295     coords[2] = PENTA18B_REF[2];
3296     break;
3297   case 1:
3298     coords[0] = PENTA18B_REF[3];
3299     coords[1] = PENTA18B_REF[4];
3300     coords[2] = PENTA18B_REF[5];
3301     break;
3302   case 2:
3303     coords[0] = PENTA18B_REF[6];
3304     coords[1] = PENTA18B_REF[7];
3305     coords[2] = PENTA18B_REF[8];
3306     break;
3307   case 3:
3308     coords[0] = PENTA18B_REF[9];
3309     coords[1] = PENTA18B_REF[10];
3310     coords[2] = PENTA18B_REF[11];
3311     break;
3312   case 4:
3313     coords[0] = PENTA18B_REF[12];
3314     coords[1] = PENTA18B_REF[13];
3315     coords[2] = PENTA18B_REF[14];
3316     break;
3317   case 5:
3318     coords[0] = PENTA18B_REF[15];
3319     coords[1] = PENTA18B_REF[16];
3320     coords[2] = PENTA18B_REF[17];
3321     break;
3322   case 6:
3323     coords[0] = PENTA18B_REF[18];
3324     coords[1] = PENTA18B_REF[19];
3325     coords[2] = PENTA18B_REF[20];
3326     break;
3327   case 7:
3328     coords[0] = PENTA18B_REF[21];
3329     coords[1] = PENTA18B_REF[22];
3330     coords[2] = PENTA18B_REF[23];
3331     break;
3332   case 8:
3333     coords[0] = PENTA18B_REF[24];
3334     coords[1] = PENTA18B_REF[25];
3335     coords[2] = PENTA18B_REF[26];
3336     break;
3337   case 9:
3338     coords[0] = PENTA18B_REF[27];
3339     coords[1] = PENTA18B_REF[28];
3340     coords[2] = PENTA18B_REF[29];
3341     break;
3342   case 10:
3343     coords[0] = PENTA18B_REF[30];
3344     coords[1] = PENTA18B_REF[31];
3345     coords[2] = PENTA18B_REF[32];
3346     break;
3347   case 11:
3348     coords[0] = PENTA18B_REF[33];
3349     coords[1] = PENTA18B_REF[34];
3350     coords[2] = PENTA18B_REF[35];
3351     break;
3352   case 12:
3353     coords[0] = PENTA18B_REF[36];
3354     coords[1] = PENTA18B_REF[37];
3355     coords[2] = PENTA18B_REF[38];
3356     break;
3357   case 13:
3358     coords[0] = PENTA18B_REF[39];
3359     coords[1] = PENTA18B_REF[40];
3360     coords[2] = PENTA18B_REF[41];
3361     break;
3362   case 14:
3363     coords[0] = PENTA18B_REF[42];
3364     coords[1] = PENTA18B_REF[43];
3365     coords[2] = PENTA18B_REF[44];
3366     break;
3367   case 15:
3368     coords[0] = PENTA18B_REF[45];
3369     coords[1] = PENTA18B_REF[46];
3370     coords[2] = PENTA18B_REF[47];
3371     break;
3372   case 16:
3373     coords[0] = PENTA18B_REF[48];
3374     coords[1] = PENTA18B_REF[49];
3375     coords[2] = PENTA18B_REF[50];
3376     break;
3377   case 17:
3378     coords[0] = PENTA18B_REF[51];
3379     coords[1] = PENTA18B_REF[52];
3380     coords[2] = PENTA18B_REF[53];
3381     break;
3382   LOCAL_COORD_MACRO_END;
3383   
3384   SHAPE_FUN_MACRO_BEGIN;  
3385   funValue[0] = gc[0]*gc[1]*(gc[0]-1.0)*(2.0*gc[1]-1.0)/2.0;
3386   funValue[2] = gc[0]*gc[2]*(gc[0]-1.0)*(2.0*gc[2]-1.0)/2.0;
3387   funValue[1] = gc[0]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3388   
3389   funValue[3] = gc[0]*gc[1]*(gc[0]+1.0)*(2.0*gc[1]-1.0)/2.0;
3390   funValue[5] = gc[0]*gc[2]*(gc[0]+1.0)*(2.0*gc[2]-1.0)/2.0;
3391   funValue[4] = gc[0]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3392   
3393   funValue[8] = 2.0*gc[0]*gc[1]*gc[2]*(gc[0]-1.0);
3394   funValue[7] = -2.0*gc[0]*gc[2]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3395   funValue[6] = -2.0*gc[0]*gc[1]*(gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3396   
3397   funValue[12] = -gc[1]*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[1]-1.0);
3398   funValue[14] = -gc[2]*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[2]-1.0);
3399   funValue[13] = -(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0);
3400   
3401   funValue[11] = 2.0*gc[0]*gc[1]*gc[2]*(gc[0]+1.0);
3402   funValue[10] = -2.0*gc[0]*gc[2]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3403   funValue[9] = -2.0*gc[0]*gc[1]*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3404   
3405   funValue[17] = -4.0*gc[1]*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3406   funValue[16] = 4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3407   funValue[15] = 4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3408   SHAPE_FUN_MACRO_END;
3409   
3410   DEV_SHAPE_FUN_MACRO_BEGIN;
3411   devFunValue[0] = gc[1]*(2.0*gc[0]-1.0)*(2.0*gc[1]-1.0)/2.0;
3412   devFunValue[1] = gc[0]*(gc[0]-1.0)*(4.0*gc[1]-1.0)/2.0;
3413   devFunValue[2] = 0.0;
3414   
3415   devFunValue[6] = gc[2]*(2.0*gc[0]-1.0)*(2.0*gc[2]-1.0)/2.0;
3416   devFunValue[7] = 0.0;
3417   devFunValue[8] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]-1.0)/2.0;
3418   
3419   devFunValue[3] = (2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3420   devFunValue[4] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3421   devFunValue[5] = gc[0]*(gc[0]-1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3422   
3423   devFunValue[9] = gc[1]*(2.0*gc[0]+1.0)*(2.0*gc[1]-1.0)/2.0;
3424   devFunValue[10] = gc[0]*(gc[0]+1.0)*(4.0*gc[1]-1.0)/2.0;
3425   devFunValue[11] = 0.0;
3426   
3427   devFunValue[15] = gc[2]*(2.0*gc[0]+1.0)*(2.0*gc[2]-1.0)/2.0;
3428   devFunValue[16] = 0.0;
3429   devFunValue[17] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]-1.0)/2.0;
3430   
3431   devFunValue[12] = (2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0)/2.0;
3432   devFunValue[13] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3433   devFunValue[14] = gc[0]*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0)/2.0;
3434   
3435   devFunValue[24] = 2.0*gc[1]*gc[2]*(2.0*gc[0]-1.0);
3436   devFunValue[25] = 2.0*gc[0]*gc[2]*(gc[0]-1.0);
3437   devFunValue[26] = 2.0*gc[0]*gc[1]*(gc[0]-1.0);
3438   
3439   devFunValue[21] = -2.0*gc[2]*(2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3440   devFunValue[22] = -2.0*gc[0]*gc[2]*(gc[0]-1.0);
3441   devFunValue[23] = -2.0*gc[0]*(gc[0]-1.0)*(2.0*gc[2]+gc[1]-1.0);
3442   
3443   devFunValue[18] = -2.0*gc[1]*(2.0*gc[0]-1.0)*(gc[2]+gc[1]-1.0);
3444   devFunValue[19] = -2.0*gc[0]*(gc[0]-1.0)*(2.0*gc[1]+gc[2]-1.0);
3445   devFunValue[20] = -2.0*gc[0]*gc[1]*(gc[0]-1.0);
3446   
3447   devFunValue[36] = -2.0*gc[0]*gc[1]*(2.0*gc[1]-1.0);
3448   devFunValue[37] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[1]-1.0);
3449   devFunValue[38] = 0.0;
3450   
3451   devFunValue[42] = -2.0*gc[0]*gc[2]*(2.0*gc[2]-1.0);
3452   devFunValue[43] = 0.0;
3453   devFunValue[44] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]-1.0);
3454   
3455   devFunValue[39] = -2.0*gc[0]*(gc[2]+gc[1]-1.0)*(2.0*gc[2]+2.0*gc[1]-1.0);
3456   devFunValue[40] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0);
3457   devFunValue[41] = -(gc[0]-1.0)*(gc[0]+1.0)*(4.0*gc[2]+4.0*gc[1]-3.0);
3458   
3459   devFunValue[33] = 2.0*gc[1]*gc[2]*(2.0*gc[0]+1.0);
3460   devFunValue[34] = 2.0*gc[0]*gc[2]*(gc[0]+1.0);
3461   devFunValue[35] = 2.0*gc[0]*gc[1]*(gc[0]+1.0);
3462   
3463   devFunValue[30] = -2.0*gc[2]*(2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3464   devFunValue[31] = -2.0*gc[0]*gc[2]*(gc[0]+1.0);
3465   devFunValue[32] = -2.0*gc[0]*(gc[0]+1.0)*(2.0*gc[2]+gc[1]-1.0);
3466   
3467   devFunValue[27] = -2.0*gc[1]*(2.0*gc[0]+1.0)*(gc[2]+gc[1]-1.0);
3468   devFunValue[28] = -2.0*gc[0]*(gc[0]+1.0)*(2.0*gc[1]+gc[2]-1.0);
3469   devFunValue[29] = -2.0*gc[0]*gc[1]*(gc[0]+1.0);
3470   
3471   devFunValue[51] = -8.0*gc[0]*gc[1]*gc[2];
3472   devFunValue[52] = -4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3473   devFunValue[53] = -4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0);
3474   
3475   devFunValue[48] = 8.0*gc[0]*gc[2]*(gc[2]+gc[1]-1.0);
3476   devFunValue[49] = 4.0*gc[2]*(gc[0]-1.0)*(gc[0]+1.0);
3477   devFunValue[50] = 4.0*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[2]+gc[1]-1.0);
3478   
3479   devFunValue[45] = 8.0*gc[0]*gc[1]*(gc[2]+gc[1]-1.0);
3480   devFunValue[46] = 4.0*(gc[0]-1.0)*(gc[0]+1.0)*(2.0*gc[1]+gc[2]-1.0);
3481   devFunValue[47] = 4.0*gc[1]*(gc[0]-1.0)*(gc[0]+1.0);
3482   DEV_SHAPE_FUN_MACRO_END;
3483 }
3484
3485 /*!
3486  * Init Hehahedron Reference coordinates and Shape function.
3487  * Case A.
3488  */
3489 void GaussInfo::hexa8aInit() 
3490 {
3491   LOCAL_COORD_MACRO_BEGIN;
3492   case 0:
3493     coords[0] = HEXA8A_REF[0];
3494     coords[1] = HEXA8A_REF[1];
3495     coords[2] = HEXA8A_REF[2];
3496     break;
3497   case 1:
3498     coords[0] = HEXA8A_REF[3];
3499     coords[1] = HEXA8A_REF[4];
3500     coords[2] = HEXA8A_REF[5];
3501     break;
3502   case 2:
3503     coords[0] = HEXA8A_REF[6];
3504     coords[1] = HEXA8A_REF[7];
3505     coords[2] = HEXA8A_REF[8];
3506     break;
3507   case 3:
3508     coords[0] = HEXA8A_REF[9];
3509     coords[1] = HEXA8A_REF[10];
3510     coords[2] = HEXA8A_REF[11];
3511     break;
3512   case 4:
3513     coords[0] = HEXA8A_REF[12];
3514     coords[1] = HEXA8A_REF[13];
3515     coords[2] = HEXA8A_REF[14];
3516     break;
3517   case 5:
3518     coords[0] = HEXA8A_REF[15];
3519     coords[1] = HEXA8A_REF[16];
3520     coords[2] = HEXA8A_REF[17];
3521     break;
3522   case 6:
3523     coords[0] = HEXA8A_REF[18];
3524     coords[1] = HEXA8A_REF[19];
3525     coords[2] = HEXA8A_REF[20];
3526     break;
3527   case 7:
3528     coords[0] = HEXA8A_REF[21];
3529     coords[1] = HEXA8A_REF[22];
3530     coords[2] = HEXA8A_REF[23];
3531     break;
3532   LOCAL_COORD_MACRO_END;
3533   
3534   SHAPE_FUN_MACRO_BEGIN;
3535   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3536   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3537   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3538   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3539   
3540   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3541   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3542   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3543   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3544   SHAPE_FUN_MACRO_END;
3545
3546   DEV_SHAPE_FUN_MACRO_BEGIN;
3547
3548   devFunValue[0] = 0.125 * (-1.0) * (1.0 - gc[1])*(1.0 - gc[2]);
3549   devFunValue[1] = 0.125 * (1.0 - gc[0]) * (-1.0) * (1.0 - gc[2]);
3550   devFunValue[2] = 0.125 * (1.0 - gc[0]) * (1.0 - gc[1]) *(-1.0);
3551
3552   devFunValue[3] = 0.125*(1.0 - gc[1])*(1.0 - gc[2]);
3553   devFunValue[4] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 - gc[2]);
3554   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0);
3555
3556   devFunValue[6] = 0.125*(1.0 + gc[1])*(1.0 - gc[2]);
3557   devFunValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[2]);
3558   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0);
3559
3560   devFunValue[9 ] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 - gc[2]);
3561   devFunValue[10] = 0.125*(1.0 - gc[0])*(1.0 - gc[2]);
3562   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0);
3563
3564   devFunValue[12] = 0.125*(-1.0)*(1.0 - gc[1])*(1.0 + gc[2]);
3565   devFunValue[13] = 0.125*(1.0 - gc[0])*(-1.0)*(1.0 + gc[2]);
3566   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1]);
3567
3568   devFunValue[15] = 0.125*(1.0 - gc[1])*(1.0 + gc[2]);
3569   devFunValue[16] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 + gc[2]);
3570   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1]);
3571
3572   devFunValue[18] = 0.125*(1.0 + gc[1])*(1.0 + gc[2]);
3573   devFunValue[19] = 0.125*(1.0 + gc[0])*(1.0 + gc[2]);
3574   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1]);
3575
3576   devFunValue[21] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 + gc[2]);
3577   devFunValue[22] = 0.125*(1.0 - gc[0])*(1.0 + gc[2]);
3578   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1]);
3579
3580   DEV_SHAPE_FUN_MACRO_END;
3581 }
3582
3583 /*!
3584  * Init Hehahedron Reference coordinates and Shape function.
3585  * Case B.
3586  */
3587 void GaussInfo::hexa8bInit() 
3588 {
3589   LOCAL_COORD_MACRO_BEGIN;
3590   case 0:
3591     coords[0] = HEXA8B_REF[0];
3592     coords[1] = HEXA8B_REF[1];
3593     coords[2] = HEXA8B_REF[2];
3594     break;
3595   case 1:
3596     coords[0] = HEXA8B_REF[3];
3597     coords[1] = HEXA8B_REF[4];
3598     coords[2] = HEXA8B_REF[5];
3599     break;
3600   case 2:
3601     coords[0] = HEXA8B_REF[6];
3602     coords[1] = HEXA8B_REF[7];
3603     coords[2] = HEXA8B_REF[8];
3604     break;
3605   case 3:
3606     coords[0] = HEXA8B_REF[9];
3607     coords[1] = HEXA8B_REF[10];
3608     coords[2] = HEXA8B_REF[11];
3609     break;
3610   case 4:
3611     coords[0] = HEXA8B_REF[12];
3612     coords[1] = HEXA8B_REF[13];
3613     coords[2] = HEXA8B_REF[14];
3614     break;
3615   case 5:
3616     coords[0] = HEXA8B_REF[15];
3617     coords[1] = HEXA8B_REF[16];
3618     coords[2] = HEXA8B_REF[17];
3619     break;
3620   case 6:
3621     coords[0] = HEXA8B_REF[18];
3622     coords[1] = HEXA8B_REF[19];
3623     coords[2] = HEXA8B_REF[20];
3624     break;
3625   case 7:
3626     coords[0] = HEXA8B_REF[21];
3627     coords[1] = HEXA8B_REF[22];
3628     coords[2] = HEXA8B_REF[23];
3629     break;
3630   LOCAL_COORD_MACRO_END;
3631   
3632   SHAPE_FUN_MACRO_BEGIN;
3633   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3634   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
3635   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3636   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
3637   
3638   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3639   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
3640   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3641   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
3642   SHAPE_FUN_MACRO_END;
3643   
3644   DEV_SHAPE_FUN_MACRO_BEGIN;
3645   devFunValue[0] = 0.125 * (-1.0) * (1.0 - gc[1])*(1.0 - gc[2]);
3646   devFunValue[1] = 0.125 * (1.0 - gc[0]) * (-1.0) * (1.0 - gc[2]);
3647   devFunValue[2] = 0.125 * (1.0 - gc[0]) * (1.0 - gc[1]) *(-1.0);
3648   
3649   devFunValue[9] = 0.125*(1.0 - gc[1])*(1.0 - gc[2]);
3650   devFunValue[10] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 - gc[2]);
3651   devFunValue[11] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0);
3652
3653   devFunValue[6] = 0.125*(1.0 + gc[1])*(1.0 - gc[2]);
3654   devFunValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[2]);
3655   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0);
3656
3657   devFunValue[3] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 - gc[2]);
3658   devFunValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[2]);
3659   devFunValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0);
3660
3661   devFunValue[12] = 0.125*(-1.0)*(1.0 - gc[1])*(1.0 + gc[2]);
3662   devFunValue[13] = 0.125*(1.0 - gc[0])*(-1.0)*(1.0 + gc[2]);
3663   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1]);
3664
3665   devFunValue[21] = 0.125*(1.0 - gc[1])*(1.0 + gc[2]);
3666   devFunValue[22] = 0.125*(1.0 + gc[0])*(-1.0)*(1.0 + gc[2]);
3667   devFunValue[23] = 0.125*(1.0 + gc[0])*(1.0 - gc[1]);
3668
3669   devFunValue[18] = 0.125*(1.0 + gc[1])*(1.0 + gc[2]);
3670   devFunValue[19] = 0.125*(1.0 + gc[0])*(1.0 + gc[2]);
3671   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1]);
3672
3673   devFunValue[15] = 0.125*(-1.0)*(1.0 + gc[1])*(1.0 + gc[2]);
3674   devFunValue[16] = 0.125*(1.0 - gc[0])*(1.0 + gc[2]);
3675   devFunValue[17] = 0.125*(1.0 - gc[0])*(1.0 + gc[1]);
3676   DEV_SHAPE_FUN_MACRO_END;
3677 }
3678
3679 /*!
3680  * This shapefunc map is same as degenerated quad4bInit
3681  */
3682 void GaussInfo::hexa8DegQuad4aInit()
3683 {
3684   LOCAL_COORD_MACRO_BEGIN;
3685  case  0:
3686    coords[0] = -1.0;
3687    coords[1] =  1.0;
3688    coords[2] =  0.0;
3689    break;
3690  case  1:
3691    coords[0] = -1.0;
3692    coords[1] = -1.0;
3693    coords[2] =  0.0;
3694    break;
3695  case  2:
3696    coords[0] =  1.0;
3697    coords[1] = -1.0;
3698    coords[2] =  0.0;
3699    break;
3700  case  3:
3701    coords[0] =  1.0;
3702    coords[1] =  1.0;
3703    coords[2] =  0.0;
3704    break;
3705  case  4:
3706    coords[0] =  0.0;
3707    coords[1] =  0.0;
3708    coords[2] =  0.0;
3709    break;
3710  case  5:
3711    coords[0] =  0.0;
3712    coords[1] =  0.0;
3713    coords[2] =  0.0;
3714    break;
3715  case  6:
3716    coords[0] =  0.0;
3717    coords[1] =  0.0;
3718    coords[2] =  0.0;
3719    break;
3720  case  7:
3721    coords[0] =  0.0;
3722    coords[1] =  0.0;
3723    coords[2] =  0.0;
3724    break;
3725   LOCAL_COORD_MACRO_END;
3726   
3727   SHAPE_FUN_MACRO_BEGIN;
3728   funValue[0] = 0.25*(1.0 + gc[1])*(1.0 - gc[0]);
3729   funValue[1] = 0.25*(1.0 - gc[1])*(1.0 - gc[0]);
3730   funValue[2] = 0.25*(1.0 - gc[1])*(1.0 + gc[0]);
3731   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3732   funValue[4] = 0.;
3733   funValue[5] = 0.;
3734   funValue[6] = 0.;
3735   funValue[7] = 0.;
3736   SHAPE_FUN_MACRO_END;
3737   
3738   DEV_SHAPE_FUN_MACRO_BEGIN;
3739   devFunValue[0] = -0.25*(1.0+gc[1]);
3740   devFunValue[1] = 0.25*(1.0-gc[0]);
3741   devFunValue[2] = 0.0;
3742   
3743   devFunValue[3] = -0.25*(1-gc[1]);
3744   devFunValue[4] = -0.25*(1.0-gc[0]);
3745   devFunValue[5] = 0.0;
3746   
3747   devFunValue[6] = 0.25*(1.0-gc[1]);
3748   devFunValue[7] = -0.25*(1.0+gc[0]);
3749   devFunValue[8] = 0.0;
3750   
3751   devFunValue[9] = 0.25*(1.0+gc[1]);
3752   devFunValue[10] = 0.25*(1.0+gc[0]);
3753   devFunValue[11] = 0.0;
3754   
3755   devFunValue[12] = 0.0;
3756   devFunValue[13] = 0.0;
3757   devFunValue[14] = 0.0;
3758   
3759   devFunValue[15] = 0.0;
3760   devFunValue[16] = 0.0;
3761   devFunValue[17] = 0.0;
3762   
3763   devFunValue[18] = 0.0;
3764   devFunValue[19] = 0.0;
3765   devFunValue[20] = 0.0;
3766   
3767   devFunValue[21] = 0.0;
3768   devFunValue[22] = 0.0;
3769   devFunValue[23] = 0.0;
3770   DEV_SHAPE_FUN_MACRO_END;
3771 }
3772
3773 /*!
3774  * This shapefunc map is same as degenerated quad4bInit
3775  */
3776 void GaussInfo::hexa8DegQuad4bInit()
3777 {
3778   LOCAL_COORD_MACRO_BEGIN;
3779  case  0:
3780    coords[0] = -1.0;
3781    coords[1] = -1.0;
3782    coords[2] =  0.0;
3783    break;
3784  case  1:
3785    coords[0] =  1.0;
3786    coords[1] = -1.0;
3787    coords[2] =  0.0;
3788    break;
3789  case  2:
3790    coords[0] =  1.0;
3791    coords[1] =  1.0;
3792    coords[2] =  0.0;
3793    break;
3794  case  3:
3795    coords[0] = -1.0;
3796    coords[1] =  1.0;
3797    coords[2] =  0.0;
3798    break;
3799  case  4:
3800    coords[0] =  0.0;
3801    coords[1] =  0.0;
3802    coords[2] =  0.0;
3803    break;
3804  case  5:
3805    coords[0] =  0.0;
3806    coords[1] =  0.0;
3807    coords[2] =  0.0;
3808    break;
3809  case  6:
3810    coords[0] =  0.0;
3811    coords[1] =  0.0;
3812    coords[2] =  0.0;
3813    break;
3814  case  7:
3815    coords[0] =  0.0;
3816    coords[1] =  0.0;
3817    coords[2] =  0.0;
3818    break;
3819    LOCAL_COORD_MACRO_END;
3820
3821   SHAPE_FUN_MACRO_BEGIN;
3822   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
3823   funValue[1] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
3824   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3825   funValue[3] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
3826   funValue[4] = 0.;
3827   funValue[5] = 0.;
3828   funValue[6] = 0.;
3829   funValue[7] = 0.;
3830   SHAPE_FUN_MACRO_END;
3831      
3832   DEV_SHAPE_FUN_MACRO_BEGIN;
3833   devFunValue[0] = -0.25*(1.0-gc[1]);
3834   devFunValue[1] = -0.25*(1.0-gc[0]);
3835   devFunValue[2] = 0.0;
3836   
3837   devFunValue[3] = 0.25*(1-gc[1]);
3838   devFunValue[4] = -0.25*(1.0+gc[0]);
3839   devFunValue[5] = 0.0;
3840   
3841   devFunValue[6] = 0.25*(1.0+gc[1]);
3842   devFunValue[7] = 0.25*(1.0+gc[0]);
3843   devFunValue[8] = 0.0;
3844   
3845   devFunValue[9] = -0.25*(1.0+gc[1]);
3846   devFunValue[10] = 0.25*(1.0-gc[0]);
3847   devFunValue[11] = 0.0;
3848   
3849   devFunValue[12] = 0.0;
3850   devFunValue[13] = 0.0;
3851   devFunValue[14] = 0.0;
3852   
3853   devFunValue[15] = 0.0;
3854   devFunValue[16] = 0.0;
3855   devFunValue[17] = 0.0;
3856   
3857   devFunValue[18] = 0.0;
3858   devFunValue[19] = 0.0;
3859   devFunValue[20] = 0.0;
3860   
3861   devFunValue[21] = 0.0;
3862   devFunValue[22] = 0.0;
3863   devFunValue[23] = 0.0;
3864   DEV_SHAPE_FUN_MACRO_END;
3865 }
3866
3867 /*!
3868  * This shapefunc map is same as degenerated quad4cInit
3869  */
3870 void GaussInfo::hexa8DegQuad4cInit() 
3871 {
3872   LOCAL_COORD_MACRO_BEGIN;
3873  case  0:
3874    coords[0] = -1.0;
3875    coords[1] = -1.0;
3876    coords[2] =  0.0;
3877    break;
3878  case  1:
3879    coords[0] = -1.0;
3880    coords[1] =  1.0;
3881    coords[2] =  0.0;
3882    break;
3883  case  2:
3884    coords[0] =  1.0;
3885    coords[1] =  1.0;
3886    coords[2] =  0.0;
3887    break;
3888  case  3:
3889    coords[0] =  1.0;
3890    coords[1] = -1.0;
3891    coords[2] =  0.0;
3892    break;
3893  case  4:
3894    coords[0] =  0.0;
3895    coords[1] =  0.0;
3896    coords[2] =  0.0;
3897    break;
3898  case  5:
3899    coords[0] =  0.0;
3900    coords[1] =  0.0;
3901    coords[2] =  0.0;
3902    break;
3903  case  6:
3904    coords[0] =  0.0;
3905    coords[1] =  0.0;
3906    coords[2] =  0.0;
3907    break;
3908  case  7:
3909    coords[0] =  0.0;
3910    coords[1] =  0.0;
3911    coords[2] =  0.0;
3912    break;
3913
3914    LOCAL_COORD_MACRO_END;
3915
3916   SHAPE_FUN_MACRO_BEGIN;
3917   funValue[0] = 0.25*(1.0 - gc[0])*(1.0 - gc[1]);
3918   funValue[1] = 0.25*(1.0 - gc[0])*(1.0 + gc[1]);
3919   funValue[2] = 0.25*(1.0 + gc[0])*(1.0 + gc[1]);
3920   funValue[3] = 0.25*(1.0 + gc[0])*(1.0 - gc[1]);
3921   funValue[4] = 0. ;
3922   funValue[5] = 0. ;
3923   funValue[6] = 0. ;
3924   funValue[7] = 0. ;
3925   SHAPE_FUN_MACRO_END;
3926    
3927   DEV_SHAPE_FUN_MACRO_BEGIN;
3928   devFunValue[0] = -0.25*(1.0-gc[1]);
3929   devFunValue[1] = -0.25*(1.0-gc[0]);
3930   devFunValue[2] = 0.0;
3931   
3932   devFunValue[3] = -0.25*(1.0+gc[1]);
3933   devFunValue[4] = 0.25*(1.0-gc[0]);
3934   devFunValue[5] = 0.0;
3935   
3936   devFunValue[6] = 0.25*(1.0+gc[1]);
3937   devFunValue[7] = 0.25*(1.0+gc[0]);
3938   devFunValue[8] = 0.0;
3939   
3940   devFunValue[9] = 0.25*(1.0-gc[1]);
3941   devFunValue[10] = -0.25*(1.0+gc[0]);
3942   devFunValue[11] = 0.0;
3943   
3944   devFunValue[12] = 0.0;
3945   devFunValue[13] = 0.0;
3946   devFunValue[14] = 0.0;
3947   
3948   devFunValue[15] = 0.0;
3949   devFunValue[16] = 0.0;
3950   devFunValue[17] = 0.0;
3951   
3952   devFunValue[18] = 0.0;
3953   devFunValue[19] = 0.0;
3954   devFunValue[20] = 0.0;
3955   
3956   devFunValue[21] = 0.0;
3957   devFunValue[22] = 0.0;
3958   devFunValue[23] = 0.0;
3959   DEV_SHAPE_FUN_MACRO_END;
3960 }
3961
3962 /*!
3963  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
3964  * Case A.
3965  */
3966 void GaussInfo::hexa20aInit() 
3967 {
3968   LOCAL_COORD_MACRO_BEGIN;
3969   case 0:
3970     coords[0] = HEXA20A_REF[0];
3971     coords[1] = HEXA20A_REF[1];
3972     coords[2] = HEXA20A_REF[2];
3973     break;
3974   case 1:
3975     coords[0] = HEXA20A_REF[3];
3976     coords[1] = HEXA20A_REF[4];
3977     coords[2] = HEXA20A_REF[5];
3978     break;
3979   case 2:
3980     coords[0] = HEXA20A_REF[6];
3981     coords[1] = HEXA20A_REF[7];
3982     coords[2] = HEXA20A_REF[8];
3983     break;
3984   case 3:
3985     coords[0] = HEXA20A_REF[9];
3986     coords[1] = HEXA20A_REF[10];
3987     coords[2] = HEXA20A_REF[11];
3988     break;
3989   case 4:
3990     coords[0] = HEXA20A_REF[12];
3991     coords[1] = HEXA20A_REF[13];
3992     coords[2] = HEXA20A_REF[14];
3993     break;
3994   case 5:
3995     coords[0] = HEXA20A_REF[15];
3996     coords[1] = HEXA20A_REF[16];
3997     coords[2] = HEXA20A_REF[17];
3998     break;
3999   case 6:
4000     coords[0] = HEXA20A_REF[18];
4001     coords[1] = HEXA20A_REF[19];
4002     coords[2] = HEXA20A_REF[20];
4003     break;
4004   case 7:
4005     coords[0] = HEXA20A_REF[21];
4006     coords[1] = HEXA20A_REF[22];
4007     coords[2] = HEXA20A_REF[23];
4008     break;
4009   case 8:
4010     coords[0] = HEXA20A_REF[24];
4011     coords[1] = HEXA20A_REF[25];
4012     coords[2] = HEXA20A_REF[26];
4013     break;
4014   case 9:
4015     coords[0] = HEXA20A_REF[27];
4016     coords[1] = HEXA20A_REF[28];
4017     coords[2] = HEXA20A_REF[29];
4018     break;
4019   case 10:
4020     coords[0] = HEXA20A_REF[30];
4021     coords[1] = HEXA20A_REF[31];
4022     coords[2] = HEXA20A_REF[32];
4023     break;
4024   case 11:
4025     coords[0] = HEXA20A_REF[33];
4026     coords[1] = HEXA20A_REF[34];
4027     coords[2] = HEXA20A_REF[35];
4028     break;
4029   case 12:
4030     coords[0] = HEXA20A_REF[36];
4031     coords[1] = HEXA20A_REF[37];
4032     coords[2] = HEXA20A_REF[38];
4033     break;
4034   case 13:
4035     coords[0] = HEXA20A_REF[39];
4036     coords[1] = HEXA20A_REF[40];
4037     coords[2] = HEXA20A_REF[41];
4038     break;
4039   case 14:
4040     coords[0] = HEXA20A_REF[42];
4041     coords[1] = HEXA20A_REF[43];
4042     coords[2] = HEXA20A_REF[44];
4043     break;
4044   case 15:
4045     coords[0] = HEXA20A_REF[45];
4046     coords[1] = HEXA20A_REF[46];
4047     coords[2] = HEXA20A_REF[47];
4048     break;
4049   case 16:
4050     coords[0] = HEXA20A_REF[48];
4051     coords[1] = HEXA20A_REF[49];
4052     coords[2] = HEXA20A_REF[50];
4053     break;
4054   case 17:
4055     coords[0] = HEXA20A_REF[51];
4056     coords[1] = HEXA20A_REF[52];
4057     coords[2] = HEXA20A_REF[53];
4058     break;
4059   case 18:
4060     coords[0] = HEXA20A_REF[54];
4061     coords[1] = HEXA20A_REF[55];
4062     coords[2] = HEXA20A_REF[56];
4063     break;
4064   case 19:
4065     coords[0] = HEXA20A_REF[57];
4066     coords[1] = HEXA20A_REF[58];
4067     coords[2] = HEXA20A_REF[59];
4068     break;
4069   LOCAL_COORD_MACRO_END;
4070   
4071   SHAPE_FUN_MACRO_BEGIN;
4072   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
4073     (-2.0 - gc[0] - gc[1] - gc[2]);
4074   funValue[1] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
4075     (-2.0 + gc[0] - gc[1] - gc[2]);
4076   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
4077     (-2.0 + gc[0] + gc[1] - gc[2]);
4078   funValue[3] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
4079     (-2.0 - gc[0] + gc[1] - gc[2]);
4080   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
4081     (-2.0 - gc[0] - gc[1] + gc[2]);
4082   funValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
4083     (-2.0 + gc[0] - gc[1] + gc[2]);
4084   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
4085     (-2.0 + gc[0] + gc[1] + gc[2]);
4086   funValue[7] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
4087     (-2.0 - gc[0] + gc[1] + gc[2]);
4088   
4089   funValue[8] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4090   funValue[9] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
4091   funValue[10] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4092   funValue[11] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
4093   funValue[12] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
4094   funValue[13] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
4095   funValue[14] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
4096   funValue[15] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
4097   funValue[16] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4098   funValue[17] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
4099   funValue[18] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4100   funValue[19] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
4101   SHAPE_FUN_MACRO_END;
4102   
4103   DEV_SHAPE_FUN_MACRO_BEGIN;
4104
4105   devFunValue[0] = 0.125*(1. + gc[1] + gc[2] + 2* gc[0]) * (1.0 - gc[1])*(1.0 - gc[2]);
4106   devFunValue[1] = 0.125*(1.0 - gc[0])*(1. + gc[0] + gc[2] + 2 * gc[1])*(1.0 - gc[2]);
4107   devFunValue[2] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[0] + gc[1] + 2 *gc[2]);
4108
4109   devFunValue[3] = 0.125*(-1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4110   devFunValue[4] = 0.125*(1.0 + gc[0])*  (1 - gc[0] + gc[2] + 2*gc[1])  *(1.0 - gc[2]);
4111   devFunValue[5] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])* (1 - gc[0] + gc[1] + 2*gc[2]);
4112   
4113   devFunValue[6] = 0.125*(-1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4114   devFunValue[7] = 0.125*(1.0 + gc[0])* (-1 + gc[0] - gc[2] + 2*gc[1]) *(1.0 - gc[2]);
4115   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[0] - gc[1] + 2*gc[2]);
4116
4117   devFunValue[9] = 0.125*(1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4118   devFunValue[10] = 0.125*(1.0 - gc[0])*(-1 - gc[0] - gc[2] + 2*gc[1])*(1.0 - gc[2]);
4119   devFunValue[11] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[0] - gc[1] + 2*gc[2]);
4120   
4121   devFunValue[12] = 0.125*(1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4122   devFunValue[13] = 0.125*(1.0 - gc[0])*(1 + gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
4123   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1] + 2*gc[2]);
4124   
4125   devFunValue[15] = 0.125*(-1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4126   devFunValue[16] = 0.125*(1.0 + gc[0])*(1 - gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
4127   devFunValue[17] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1] + 2*gc[2]);
4128   
4129   devFunValue[18] = 0.125*(-1.0 + gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4130   devFunValue[19] = 0.125*(1.0 + gc[0])*(-1 + gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
4131   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1] + 2*gc[2]);
4132
4133   devFunValue[21] = 0.125*(1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4134   devFunValue[22] = 0.125*(1.0 - gc[0])*(-1 - gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
4135   devFunValue[23] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1] + 2*gc[2]);
4136
4137   devFunValue[24] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4138   devFunValue[25] = 0.25*(1.0 - gc[0]*gc[0])*(-1)*(1.0 - gc[2]);
4139   devFunValue[26] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(-1.0);
4140
4141   devFunValue[27] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[2]);
4142   devFunValue[28] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
4143   devFunValue[29] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(-1.0);
4144
4145   devFunValue[30] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4146   devFunValue[31] = 0.25*(1.0 - gc[0]*gc[0])*(1.0)*(1.0 - gc[2]);
4147   devFunValue[32] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(-1.0);
4148
4149   devFunValue[33] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 - gc[2]);
4150   devFunValue[34] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
4151   devFunValue[35] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(-1.0);
4152
4153   devFunValue[36] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 - gc[1]);
4154   devFunValue[37] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(-1.0);
4155   devFunValue[38] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
4156
4157   devFunValue[39] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[1]);
4158   devFunValue[40] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(-1.0);
4159   devFunValue[41] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
4160
4161   devFunValue[42] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[1]);
4162   devFunValue[43] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0]);
4163   devFunValue[44] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
4164
4165   devFunValue[45] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 + gc[1]);
4166   devFunValue[46] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0]);
4167   devFunValue[47] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
4168
4169   devFunValue[48] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4170   devFunValue[49] = 0.25*(1.0 - gc[0]*gc[0])*(-1.0)*(1.0 + gc[2]);
4171   devFunValue[50] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
4172
4173   devFunValue[51] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[2]);
4174   devFunValue[52] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
4175   devFunValue[53] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
4176
4177   devFunValue[54] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4178   devFunValue[55] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[2]);
4179   devFunValue[56] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
4180
4181   devFunValue[57] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 + gc[2]);
4182   devFunValue[58] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
4183   devFunValue[59] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
4184
4185   SHAPE_FUN_MACRO_END;
4186 }
4187
4188 /*!
4189  * Init Qaudratic Hehahedron Reference coordinates and Shape function.
4190  * Case B.
4191  */
4192 void GaussInfo::hexa20bInit()
4193 {
4194   LOCAL_COORD_MACRO_BEGIN;
4195   case 0:
4196     coords[0] = HEXA20B_REF[0];
4197     coords[1] = HEXA20B_REF[1];
4198     coords[2] = HEXA20B_REF[2];
4199     break;
4200   case 1:
4201     coords[0] = HEXA20B_REF[3];
4202     coords[1] = HEXA20B_REF[4];
4203     coords[2] = HEXA20B_REF[5];
4204     break;
4205   case 2:
4206     coords[0] = HEXA20B_REF[6];
4207     coords[1] = HEXA20B_REF[7];
4208     coords[2] = HEXA20B_REF[8];
4209     break;
4210   case 3:
4211     coords[0] = HEXA20B_REF[9];
4212     coords[1] = HEXA20B_REF[10];
4213     coords[2] = HEXA20B_REF[11];
4214     break;
4215   case 4:
4216     coords[0] = HEXA20B_REF[12];
4217     coords[1] = HEXA20B_REF[13];
4218     coords[2] = HEXA20B_REF[14];
4219     break;
4220   case 5:
4221     coords[0] = HEXA20B_REF[15];
4222     coords[1] = HEXA20B_REF[16];
4223     coords[2] = HEXA20B_REF[17];
4224     break;
4225   case 6:
4226     coords[0] = HEXA20B_REF[18];
4227     coords[1] = HEXA20B_REF[19];
4228     coords[2] = HEXA20B_REF[20];
4229     break;
4230   case 7:
4231     coords[0] = HEXA20B_REF[21];
4232     coords[1] = HEXA20B_REF[22];
4233     coords[2] = HEXA20B_REF[23];
4234     break;
4235   case 8:
4236     coords[0] = HEXA20B_REF[24];
4237     coords[1] = HEXA20B_REF[25];
4238     coords[2] = HEXA20B_REF[26];
4239     break;
4240   case 9:
4241     coords[0] = HEXA20B_REF[27];
4242     coords[1] = HEXA20B_REF[28];
4243     coords[2] = HEXA20B_REF[29];
4244     break;
4245   case 10:
4246     coords[0] = HEXA20B_REF[30];
4247     coords[1] = HEXA20B_REF[31];
4248     coords[2] = HEXA20B_REF[32];
4249     break;
4250   case 11:
4251     coords[0] = HEXA20B_REF[33];
4252     coords[1] = HEXA20B_REF[34];
4253     coords[2] = HEXA20B_REF[35];
4254     break;
4255   case 12:
4256     coords[0] = HEXA20B_REF[36];
4257     coords[1] = HEXA20B_REF[37];
4258     coords[2] = HEXA20B_REF[38];
4259     break;
4260   case 13:
4261     coords[0] = HEXA20B_REF[39];
4262     coords[1] = HEXA20B_REF[40];
4263     coords[2] = HEXA20B_REF[41];
4264     break;
4265   case 14:
4266     coords[0] = HEXA20B_REF[42];
4267     coords[1] = HEXA20B_REF[43];
4268     coords[2] = HEXA20B_REF[44];
4269     break;
4270   case 15:
4271     coords[0] = HEXA20B_REF[45];
4272     coords[1] = HEXA20B_REF[46];
4273     coords[2] = HEXA20B_REF[47];
4274     break;
4275   case 16:
4276     coords[0] = HEXA20B_REF[48];
4277     coords[1] = HEXA20B_REF[49];
4278     coords[2] = HEXA20B_REF[50];
4279     break;
4280   case 17:
4281     coords[0] = HEXA20B_REF[51];
4282     coords[1] = HEXA20B_REF[52];
4283     coords[2] = HEXA20B_REF[53];
4284     break;
4285   case 18:
4286     coords[0] = HEXA20B_REF[54];
4287     coords[1] = HEXA20B_REF[55];
4288     coords[2] = HEXA20B_REF[56];
4289     break;
4290   case 19:
4291     coords[0] = HEXA20B_REF[57];
4292     coords[1] = HEXA20B_REF[58];
4293     coords[2] = HEXA20B_REF[59];
4294     break;
4295   LOCAL_COORD_MACRO_END;
4296   
4297   SHAPE_FUN_MACRO_BEGIN;
4298   
4299   funValue[0] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
4300     (-2.0 - gc[0] - gc[1] - gc[2]);
4301   funValue[3] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 - gc[2])*
4302     (-2.0 + gc[0] - gc[1] - gc[2]);
4303   funValue[2] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
4304     (-2.0 + gc[0] + gc[1] - gc[2]);
4305   funValue[1] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 - gc[2])*
4306     (-2.0 - gc[0] + gc[1] - gc[2]);
4307   funValue[4] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
4308     (-2.0 - gc[0] - gc[1] + gc[2]);
4309   funValue[7] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(1.0 + gc[2])*
4310     (-2.0 + gc[0] - gc[1] + gc[2]);
4311   funValue[6] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
4312     (-2.0 + gc[0] + gc[1] + gc[2]);
4313   funValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[2])*
4314     (-2.0 - gc[0] + gc[1] + gc[2]);
4315   
4316   funValue[11] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4317   funValue[10] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
4318   funValue[9] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4319   funValue[8] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
4320   funValue[16] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
4321   funValue[19] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
4322   funValue[18] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
4323   funValue[17] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
4324   funValue[15] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4325   funValue[14] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
4326   funValue[13] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4327   funValue[12] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
4328   SHAPE_FUN_MACRO_END;
4329   
4330   DEV_SHAPE_FUN_MACRO_BEGIN;
4331   devFunValue[0] = 0.125*(1. + gc[1] + gc[2] + 2* gc[0]) * (1.0 - gc[1])*(1.0 - gc[2]);
4332   devFunValue[1] = 0.125*(1.0 - gc[0])*(1. + gc[0] + gc[2] + 2 * gc[1])*(1.0 - gc[2]);
4333   devFunValue[2] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(1.0 + gc[0] + gc[1] + 2 *gc[2]);
4334
4335   devFunValue[9] = 0.125*(-1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4336   devFunValue[10] = 0.125*(1.0 + gc[0])*  (1 - gc[0] + gc[2] + 2*gc[1])  *(1.0 - gc[2]);
4337   devFunValue[11] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])* (1 - gc[0] + gc[1] + 2*gc[2]);
4338   
4339   devFunValue[6] = 0.125*(-1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4340   devFunValue[7] = 0.125*(1.0 + gc[0])* (-1 + gc[0] - gc[2] + 2*gc[1]) *(1.0 - gc[2]);
4341   devFunValue[8] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(1.0 - gc[0] - gc[1] + 2*gc[2]);
4342
4343   devFunValue[3] = 0.125*(1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4344   devFunValue[4] = 0.125*(1.0 - gc[0])*(-1 - gc[0] - gc[2] + 2*gc[1])*(1.0 - gc[2]);
4345   devFunValue[5] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(1.0 + gc[0] - gc[1] + 2*gc[2]);
4346   
4347   devFunValue[12] = 0.125*(1.0 + gc[1] - gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4348   devFunValue[13] = 0.125*(1.0 - gc[0])*(1 + gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
4349   devFunValue[14] = 0.125*(1.0 - gc[0])*(1.0 - gc[1])*(-1.0 - gc[0] - gc[1] + 2*gc[2]);
4350   
4351   devFunValue[21] = 0.125*(-1.0 - gc[1] + gc[2] + 2*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4352   devFunValue[22] = 0.125*(1.0 + gc[0])*(1 - gc[0] - gc[2] + 2*gc[1])*(1.0 + gc[2]);
4353   devFunValue[23] = 0.125*(1.0 + gc[0])*(1.0 - gc[1])*(-1.0 + gc[0] - gc[1] + 2*gc[2]);
4354   
4355   devFunValue[18] = 0.125*(-1.0 + gc[1] + gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4356   devFunValue[19] = 0.125*(1.0 + gc[0])*(-1 + gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
4357   devFunValue[20] = 0.125*(1.0 + gc[0])*(1.0 + gc[1])*(-1.0 + gc[0] + gc[1] + 2*gc[2]);
4358
4359   devFunValue[15] = 0.125*(1.0 - gc[1] - gc[2] + 2*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4360   devFunValue[16] = 0.125*(1.0 - gc[0])*(-1 - gc[0] + gc[2] + 2*gc[1])*(1.0 + gc[2]);
4361   devFunValue[17] = 0.125*(1.0 - gc[0])*(1.0 + gc[1])*(-1.0 - gc[0] + gc[1] + 2*gc[2]);
4362
4363   devFunValue[33] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 - gc[2]);
4364   devFunValue[34] = 0.25*(1.0 - gc[0]*gc[0])*(-1)*(1.0 - gc[2]);
4365   devFunValue[35] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1])*(-1.0);
4366
4367   devFunValue[30] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[2]);
4368   devFunValue[31] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 - gc[2]);
4369   devFunValue[32] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0])*(-1.0);
4370
4371   devFunValue[27] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 - gc[2]);
4372   devFunValue[28] = 0.25*(1.0 - gc[0]*gc[0])*(1.0)*(1.0 - gc[2]);
4373   devFunValue[29] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1])*(-1.0);
4374
4375   devFunValue[24] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 - gc[2]);
4376   devFunValue[25] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 - gc[2]);
4377   devFunValue[26] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0])*(-1.0);
4378
4379   devFunValue[48] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 - gc[1]);
4380   devFunValue[49] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0])*(-1.0);
4381   devFunValue[50] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 - gc[1]);
4382
4383   devFunValue[57] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[1]);
4384   devFunValue[58] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0])*(-1.0);
4385   devFunValue[59] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 - gc[1]);
4386
4387   devFunValue[54] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[1]);
4388   devFunValue[55] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 + gc[0]);
4389   devFunValue[56] = 0.25*(-2.0*gc[2])*(1.0 + gc[0])*(1.0 + gc[1]);
4390
4391   devFunValue[51] = 0.25*(1.0 - gc[2]*gc[2])*(-1.0)*(1.0 + gc[1]);
4392   devFunValue[52] = 0.25*(1.0 - gc[2]*gc[2])*(1.0 - gc[0]);
4393   devFunValue[53] = 0.25*(-2.0*gc[2])*(1.0 - gc[0])*(1.0 + gc[1]);
4394
4395   devFunValue[45] = 0.25*(-2.0*gc[0])*(1.0 - gc[1])*(1.0 + gc[2]);
4396   devFunValue[46] = 0.25*(1.0 - gc[0]*gc[0])*(-1.0)*(1.0 + gc[2]);
4397   devFunValue[47] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 - gc[1]);
4398
4399   devFunValue[42] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[2]);
4400   devFunValue[43] = 0.25*(-2.0*gc[1])*(1.0 + gc[0])*(1.0 + gc[2]);
4401   devFunValue[44] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 + gc[0]);
4402
4403   devFunValue[39] = 0.25*(-2.0*gc[0])*(1.0 + gc[1])*(1.0 + gc[2]);
4404   devFunValue[40] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[2]);
4405   devFunValue[41] = 0.25*(1.0 - gc[0]*gc[0])*(1.0 + gc[1]);
4406
4407   devFunValue[36] = 0.25*(1.0 - gc[1]*gc[1])*(-1.0)*(1.0 + gc[2]);
4408   devFunValue[37] = 0.25*(-2.0*gc[1])*(1.0 - gc[0])*(1.0 + gc[2]);
4409   devFunValue[38] = 0.25*(1.0 - gc[1]*gc[1])*(1.0 - gc[0]);
4410   DEV_SHAPE_FUN_MACRO_END;
4411 }
4412
4413 void GaussInfo::hexa27aInit()
4414 {
4415   LOCAL_COORD_MACRO_BEGIN;
4416   case 0:
4417     coords[0] = HEXA27A_REF[0];
4418     coords[1] = HEXA27A_REF[1];
4419     coords[2] = HEXA27A_REF[2];
4420     break;
4421   case 1:
4422     coords[0] = HEXA27A_REF[3];
4423     coords[1] = HEXA27A_REF[4];
4424     coords[2] = HEXA27A_REF[5];
4425     break;
4426   case 2:
4427     coords[0] = HEXA27A_REF[6];
4428     coords[1] = HEXA27A_REF[7];
4429     coords[2] = HEXA27A_REF[8];
4430     break;
4431   case 3:
4432     coords[0] = HEXA27A_REF[9];
4433     coords[1] = HEXA27A_REF[10];
4434     coords[2] = HEXA27A_REF[11];
4435     break;
4436   case 4:
4437     coords[0] = HEXA27A_REF[12];
4438     coords[1] = HEXA27A_REF[13];
4439     coords[2] = HEXA27A_REF[14];
4440     break;
4441   case 5:
4442     coords[0] = HEXA27A_REF[15];
4443     coords[1] = HEXA27A_REF[16];
4444     coords[2] = HEXA27A_REF[17];
4445     break;
4446   case 6:
4447     coords[0] = HEXA27A_REF[18];
4448     coords[1] = HEXA27A_REF[19];
4449     coords[2] = HEXA27A_REF[20];
4450     break;
4451   case 7:
4452     coords[0] = HEXA27A_REF[21];
4453     coords[1] = HEXA27A_REF[22];
4454     coords[2] = HEXA27A_REF[23];
4455     break;
4456   case 8:
4457     coords[0] = HEXA27A_REF[24];
4458     coords[1] = HEXA27A_REF[25];
4459     coords[2] = HEXA27A_REF[26];
4460     break;
4461   case 9:
4462     coords[0] = HEXA27A_REF[27];
4463     coords[1] = HEXA27A_REF[28];
4464     coords[2] = HEXA27A_REF[29];
4465     break;
4466   case 10:
4467     coords[0] = HEXA27A_REF[30];
4468     coords[1] = HEXA27A_REF[31];
4469     coords[2] = HEXA27A_REF[32];
4470     break;
4471   case 11:
4472     coords[0] = HEXA27A_REF[33];
4473     coords[1] = HEXA27A_REF[34];
4474     coords[2] = HEXA27A_REF[35];
4475     break;
4476   case 12:
4477     coords[0] = HEXA27A_REF[36];
4478     coords[1] = HEXA27A_REF[37];
4479     coords[2] = HEXA27A_REF[38];
4480     break;
4481   case 13:
4482     coords[0] = HEXA27A_REF[39];
4483     coords[1] = HEXA27A_REF[40];
4484     coords[2] = HEXA27A_REF[41];
4485     break;
4486   case 14:
4487     coords[0] = HEXA27A_REF[42];
4488     coords[1] = HEXA27A_REF[43];
4489     coords[2] = HEXA27A_REF[44];
4490     break;
4491   case 15:
4492     coords[0] = HEXA27A_REF[45];
4493     coords[1] = HEXA27A_REF[46];
4494     coords[2] = HEXA27A_REF[47];
4495     break;
4496   case 16:
4497     coords[0] = HEXA27A_REF[48];
4498     coords[1] = HEXA27A_REF[49];
4499     coords[2] = HEXA27A_REF[50];
4500     break;
4501   case 17:
4502     coords[0] = HEXA27A_REF[51];
4503     coords[1] = HEXA27A_REF[52];
4504     coords[2] = HEXA27A_REF[53];
4505     break;
4506   case 18:
4507     coords[0] = HEXA27A_REF[54];
4508     coords[1] = HEXA27A_REF[55];
4509     coords[2] = HEXA27A_REF[56];
4510     break;
4511   case 19:
4512     coords[0] = HEXA27A_REF[57];
4513     coords[1] = HEXA27A_REF[58];
4514     coords[2] = HEXA27A_REF[59];
4515     break;
4516   case 20:
4517     coords[0] = HEXA27A_REF[60];
4518     coords[1] = HEXA27A_REF[61];
4519     coords[2] = HEXA27A_REF[62];
4520     break;
4521   case 21:
4522     coords[0] = HEXA27A_REF[63];
4523     coords[1] = HEXA27A_REF[64];
4524     coords[2] = HEXA27A_REF[65];
4525     break;
4526   case 22:
4527     coords[0] = HEXA27A_REF[66];
4528     coords[1] = HEXA27A_REF[67];
4529     coords[2] = HEXA27A_REF[68];
4530     break;
4531   case 23:
4532     coords[0] = HEXA27A_REF[69];
4533     coords[1] = HEXA27A_REF[70];
4534     coords[2] = HEXA27A_REF[71];
4535     break;
4536   case 24:
4537     coords[0] = HEXA27A_REF[72];
4538     coords[1] = HEXA27A_REF[73];
4539     coords[2] = HEXA27A_REF[74];
4540     break;
4541   case 25:
4542     coords[0] = HEXA27A_REF[75];
4543     coords[1] = HEXA27A_REF[76];
4544     coords[2] = HEXA27A_REF[77];
4545     break;
4546   case 26:
4547     coords[0] = HEXA27A_REF[78];
4548     coords[1] = HEXA27A_REF[79];
4549     coords[2] = HEXA27A_REF[80];
4550     break;
4551   LOCAL_COORD_MACRO_END;
4552   
4553   SHAPE_FUN_MACRO_BEGIN;
4554   
4555   funValue[0] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4556   funValue[1] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4557   funValue[2] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4558   funValue[3] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4559   funValue[4] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4560   funValue[5] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4561   funValue[6] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4562   funValue[7] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4563   funValue[8] =0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4564   funValue[9] =0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4565   funValue[10]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4566   funValue[11]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4567   funValue[12]=0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4568   funValue[13]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4569   funValue[14]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4570   funValue[15]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4571   funValue[16]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4572   funValue[17]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4573   funValue[18]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4574   funValue[19]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4575   funValue[20]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4576   funValue[21]=0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4577   funValue[22]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4578   funValue[23]=0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4579   funValue[24]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4580   funValue[25]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4581   funValue[26]=(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4582   
4583   SHAPE_FUN_MACRO_END;
4584   
4585   DEV_SHAPE_FUN_MACRO_BEGIN;
4586
4587   devFunValue[0] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4588   devFunValue[1] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
4589   devFunValue[2] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
4590
4591   devFunValue[3] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4592   devFunValue[4] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
4593   devFunValue[5] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
4594
4595   devFunValue[6] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4596   devFunValue[7] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
4597   devFunValue[8] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
4598
4599   devFunValue[9]  = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4600   devFunValue[10] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
4601   devFunValue[11] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
4602
4603   devFunValue[12] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4604   devFunValue[13] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
4605   devFunValue[14] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
4606
4607   devFunValue[15] = 0.125*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4608   devFunValue[16] = 0.125*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
4609   devFunValue[17] = 0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
4610
4611   devFunValue[18] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4612   devFunValue[19] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
4613   devFunValue[20] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
4614
4615   devFunValue[21] = 0.125*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4616   devFunValue[22] = 0.125*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
4617   devFunValue[23] = 0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
4618
4619   devFunValue[24] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4620   devFunValue[25] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
4621   devFunValue[26] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
4622
4623   devFunValue[27] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
4624   devFunValue[28] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]-1.);
4625   devFunValue[29] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]-1.);
4626   
4627   devFunValue[30] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4628   devFunValue[31] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]-1.);
4629   devFunValue[32] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
4630
4631   devFunValue[33] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
4632   devFunValue[34] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]-1.);
4633   devFunValue[35] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]-1.);
4634
4635   devFunValue[36] = 0.25*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4636   devFunValue[37] = 0.25*gc[0]*(gc[0]-1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
4637   devFunValue[38] = 0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
4638
4639   devFunValue[39] = 0.25*(-2*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
4640   devFunValue[40] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*gc[2]*(gc[2]+1.);
4641   devFunValue[41] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(2*gc[2]+1.);
4642
4643   devFunValue[42] = 0.25*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4644   devFunValue[43] = 0.25*gc[0]*(gc[0]+1.)*(-2*gc[1])*gc[2]*(gc[2]+1.);
4645   devFunValue[44] = 0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
4646
4647   devFunValue[45] = 0.25*(-2*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
4648   devFunValue[46] = 0.25*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*gc[2]*(gc[2]+1.);
4649   devFunValue[47] = 0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(2*gc[2]+1.);
4650
4651   devFunValue[48] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4652   devFunValue[49] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
4653   devFunValue[50] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
4654
4655   devFunValue[51] = 0.25*(2*gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4656   devFunValue[52] = 0.25*gc[0]*(gc[0]-1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
4657   devFunValue[53] = 0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
4658
4659   devFunValue[54] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4660   devFunValue[55] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
4661   devFunValue[56] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(-2*gc[2]);
4662
4663   devFunValue[57] = 0.25*(2*gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4664   devFunValue[58] = 0.25*gc[0]*(gc[0]+1.)*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
4665   devFunValue[59] = 0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(-2*gc[2]);
4666
4667   devFunValue[60] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
4668   devFunValue[61] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]-1.);
4669   devFunValue[62] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]-1.);
4670
4671   devFunValue[63] = 0.5*(2*gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4672   devFunValue[64] = 0.5*gc[0]*(gc[0]-1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
4673   devFunValue[65] = 0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
4674
4675   devFunValue[66] = 0.5*(-2*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
4676   devFunValue[67] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]+1.)*(1.-gc[2]*gc[2]);
4677   devFunValue[68] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(-2*gc[2]);
4678
4679   devFunValue[69] = 0.5*(2*gc[0]+1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4680   devFunValue[70] = 0.5*gc[0]*(gc[0]+1.)*(-2*gc[1])*(1.-gc[2]*gc[2]);
4681   devFunValue[71] = 0.5*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*(-2*gc[2]);
4682
4683   devFunValue[72] = 0.5*(-2*gc[0])*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
4684   devFunValue[73] = 0.5*(1.-gc[0]*gc[0])*(2*gc[1]-1.)*(1.-gc[2]*gc[2]);
4685   devFunValue[74] = 0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*(-2*gc[2]);
4686
4687   devFunValue[75] = 0.5*(-2*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
4688   devFunValue[76] = 0.5*(1.-gc[0]*gc[0])*(-2*gc[1])*gc[2]*(gc[2]+1.);
4689   devFunValue[77] = 0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(2*gc[2]+1.);
4690
4691   devFunValue[78] = (-2*gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
4692   devFunValue[79] = (1.-gc[0]*gc[0])*(-2*gc[1])*(1.-gc[2]*gc[2]);
4693   devFunValue[80] = (1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*(-2*gc[2]);
4694
4695   DEV_SHAPE_FUN_MACRO_END;
4696 }
4697
4698 ////////////////////////////////////////////////////////////////////////////////////////////////
4699 //                                GAUSS COORD CLASS                                           //
4700 ////////////////////////////////////////////////////////////////////////////////////////////////
4701 /*!
4702  * Constructor
4703  */
4704 GaussCoords::GaussCoords()
4705 {
4706 }
4707
4708 /*!
4709  * Destructor
4710  */
4711 GaussCoords::~GaussCoords()
4712 {
4713   GaussInfoVector::iterator it = _my_gauss_info.begin();
4714   for( ; it != _my_gauss_info.end(); it++ ) 
4715     {
4716       if((*it) != NULL)
4717         delete (*it);
4718     }
4719 }
4720
4721 /*!
4722  * Add Gauss localization info 
4723  */
4724 void GaussCoords::addGaussInfo( NormalizedCellType theGeometry,
4725                                 int coordDim,
4726                                 const double* theGaussCoord,
4727                                 mcIdType theNbGauss,
4728                                 const double* theReferenceCoord,
4729                                 mcIdType theNbRef)
4730 {
4731   GaussInfoVector::iterator it = _my_gauss_info.begin();
4732   for( ; it != _my_gauss_info.end(); it++ ) 
4733     {
4734       if( (*it)->getCellType() == theGeometry ) 
4735         {
4736           break;
4737         }
4738     }
4739
4740   DataVector aGaussCoord;
4741   for(int i = 0 ; i < theNbGauss*coordDim; i++ )
4742     aGaussCoord.push_back(theGaussCoord[i]);
4743
4744   DataVector aReferenceCoord;
4745   for(int i = 0 ; i < theNbRef*coordDim; i++ )
4746     aReferenceCoord.push_back(theReferenceCoord[i]);
4747
4748
4749   GaussInfo* info = new GaussInfo( theGeometry, aGaussCoord, FromIdType<int>(theNbGauss), aReferenceCoord, FromIdType<int>(theNbRef));
4750   info->initLocalInfo();
4751
4752   //If info with cell type doesn't exist add it
4753   if( it == _my_gauss_info.end() ) 
4754     {
4755       _my_gauss_info.push_back(info);
4756
4757       // If information exists, update it
4758     }
4759   else 
4760     {
4761       std::size_t index = std::distance(_my_gauss_info.begin(),it);
4762       delete (*it);
4763       _my_gauss_info[index] = info;
4764     }
4765 }
4766
4767
4768 /*!
4769  * Calculate gauss points coordinates
4770  */
4771 double* GaussCoords::calculateCoords( NormalizedCellType theGeometry, 
4772                                       const double *theNodeCoords, 
4773                                       const int theSpaceDim,
4774                                       const mcIdType *theIndex)
4775 {
4776   const GaussInfo *info = getInfoGivenCellType(theGeometry);
4777   int nbCoords = theSpaceDim * info->getNbGauss();
4778   double *aCoords = new double[nbCoords];
4779   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,aCoords);
4780   return aCoords;
4781 }
4782
4783
4784 void GaussCoords::calculateCoords( NormalizedCellType theGeometry, const double *theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
4785 {
4786   const GaussInfo *info = getInfoGivenCellType(theGeometry);
4787   calculateCoordsAlg(info,theNodeCoords,theSpaceDim,theIndex,result);
4788 }
4789
4790 void GaussCoords::calculateCoordsAlg(const GaussInfo *info, const double* theNodeCoords, const int theSpaceDim, const mcIdType *theIndex, double *result)
4791 {
4792   int aConn = info->getNbRef();
4793
4794   int nbCoords = theSpaceDim * info->getNbGauss();
4795   std::fill(result,result+nbCoords,0.);
4796
4797   for( int gaussId = 0; gaussId < info->getNbGauss(); gaussId++ ) 
4798     {
4799       double *coord=result+gaussId*theSpaceDim;
4800       const double *function=info->getFunctionValues(gaussId);
4801       for ( int connId = 0; connId < aConn ; connId++ ) 
4802         {
4803           const double* nodeCoord = theNodeCoords + (theIndex[connId]*theSpaceDim);
4804           for( int dimId = 0; dimId < theSpaceDim; dimId++ )
4805             coord[dimId] += nodeCoord[dimId]*function[connId];
4806         }
4807     }
4808 }
4809
4810 const GaussInfo *GaussCoords::getInfoGivenCellType(NormalizedCellType cellType)
4811 {
4812   GaussInfoVector::const_iterator it = _my_gauss_info.begin();
4813   //Try to find gauss localization info
4814   for( ; it != _my_gauss_info.end() ; it++ ) 
4815     if( (*it)->getCellType()==cellType) 
4816       return (*it);
4817   throw INTERP_KERNEL::Exception("Can't find gauss localization information !");
4818 }