diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h
index 53176cdd21d80165f30573bd8ae026c294e9ba95..678fa8bc1475a765e3b34b843f2d4944507eefac 100644
--- a/Post/shapeFunctions.h
+++ b/Post/shapeFunctions.h
@@ -830,22 +830,24 @@ public:
   inline int getNumGaussPoints(){ return 8; }
   void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
   {
-    double u8[8] = {0.3595161057791018,0.09633205020967324,
-                    0.3595161057791018,0.09633205020967324,
-                    0.6920507403468987,0.1854344369976602,
-                    0.6920507403468987,0.1854344369976602};
-    double v8[8] = {0.3595161057791018,0.3595161057791018,
-                    0.09633205020967324,0.09633205020967324,
-                    0.6920507403468987,0.6920507403468987,
-                    0.1854344369976602,0.1854344369976602};
+    //Integration points for old_pyramids
+    double u8[8] = {0.2631840555694285,-0.2631840555694285,
+		    0.2631840555694285,-0.2631840555694285,
+		    0.5066163033492386,-0.5066163033492386,
+		    0.5066163033492386,-0.5066163033492386};
+    double v8[8] = {0.2631840555694285,0.2631840555694285,
+		    -0.2631840555694285,-0.2631840555694285,
+		    0.5066163033492386,0.5066163033492386,
+		    -0.5066163033492386,-0.5066163033492386};
     double w8[8] = {0.544151844011225,0.544151844011225,
-                    0.544151844011225,0.544151844011225,
-                    0.122514822655441,0.122514822655441,
-                    0.122514822655441,0.122514822655441};
-    double p8[8] = {0.02519647051995625,0.02519647051995625,
-                    0.02519647051995625,0.02519647051995625,
-                    0.058136862813377,0.058136862813377,
-                    0.058136862813377,0.058136862813377};
+		    0.544151844011225,0.544151844011225,
+		    0.122514822655441,0.122514822655441,
+		    0.122514822655441,0.122514822655441};
+    double p8[8] = {0.100785882079825,0.100785882079825,
+		    0.100785882079825,0.100785882079825,
+		    0.232547451253508,0.232547451253508,
+		    0.232547451253508,0.232547451253508};
+
     if(num < 0 || num > 7) return;
     u = u8[num];
     v = v8[num];
@@ -855,6 +857,7 @@ public:
   void getShapeFunction(int num, double u, double v, double w, double &s) 
   {
     double r;
+
     if(w != 1. && num != 4) r = u*v*w / (1.-w);
     else                    r = 0.;
     switch(num) {
@@ -862,31 +865,46 @@ public:
     case 1  : s = 0.25 * ((1.+u) * (1.-v) - w - r); break;
     case 2  : s = 0.25 * ((1.+u) * (1.+v) - w + r); break;
     case 3  : s = 0.25 * ((1.-u) * (1.+v) - w - r); break;
-    case 4  : s = w; break;
-    default : s = 0.; break;
+    case 4  : s = w  ; break;
+    default : s = 0. ; break;
     }
   }
   void getGradShapeFunction(int num, double u, double v, double w, double s[3]) 
   {
     if(w == 1. && num != 4) {
-      s[0] =  0.25; 
-      s[1] =  0.25;
-      s[2] = -0.25; 
+      switch(num) {
+      case 0 : s[0] = -0.25 ; 
+	       s[1] = -0.25 ;
+	       s[2] = -0.25 ; break ;
+      case 1 : s[0] =  0.25 ; 
+	       s[1] = -0.25 ;
+	       s[2] = -0.25 ; break ;
+      case 2 : s[0] =  0.25 ; 
+	       s[1] =  0.25 ;
+               s[2] = -0.25 ; break ;
+      case 3 : s[0] = -0.25 ; 
+	       s[1] =  0.25 ;
+               s[2] = -0.25 ; break ;
+      case 4 : s[0] =  0. ; 
+	       s[1] =  0. ;
+	       s[2] =  1. ; break ;
+      default : s[0] = s[1] = s[2] = 0.; break;
+      }
     }
     else{
       switch(num) {
       case 0  : s[0] = 0.25 * ( -(1.-v) + v*w/(1.-w) ) ;
                 s[1] = 0.25 * ( -(1.-u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/SQU(1.-w) ) ; break ;
-      case 1  : s[0] = 0.25 * (  (1.-v) + v*w/(1.-w) ) ;
-                s[1] = 0.25 * ( -(1.+u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
+      case 1  : s[0] = 0.25 * (  (1.-v) - v*w/(1.-w) ) ;
+                s[1] = 0.25 * ( -(1.+u) - u*w/(1.-w) ) ;
+                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
       case 2  : s[0] = 0.25 * (  (1.+v) + v*w/(1.-w) ) ;
                 s[1] = 0.25 * (  (1.+u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/SQU(1.-w) ) ; break ;
-      case 3  : s[0] = 0.25 * ( -(1.+v) + v*w/(1.-w) ) ;
-                s[1] = 0.25 * (  (1.-u) + u*w/(1.-w) ) ;
-                s[2] = 0.25 * ( -1.     + u*v/SQU(1.-w) ) ; break ;
+                s[2] = 0.25 * ( -1.     + u*v/(1.-w) + u*v*w/SQU(1.-w) ) ; break ;
+      case 3  : s[0] = 0.25 * ( -(1.+v) - v*w/(1.-w) ) ;
+                s[1] = 0.25 * (  (1.-u) - u*w/(1.-w) ) ;
+                s[2] = 0.25 * ( -1.     - u*v/(1.-w) - u*v*w/SQU(1.-w) ) ; break ;
       case 4  : s[0] = 0. ; 
                 s[1] = 0. ;
                 s[2] = 1. ; break ;