diff --git a/Geo/MElement.h b/Geo/MElement.h
index ce3c8a8f41740083f4adc6e6af85bef165ea7fb2..5181230ae0546bc58ddf14d7ea4664555900347a 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -2550,25 +2550,32 @@ class MPyramid : public MElement {
   virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o) 
   {
     if(w == 1.) {
-      for(int i = 0; i < 4; i++){
-        s[i][0] =  0.25; 
-        s[i][1] =  0.25;
-        s[i][2] = -0.25; 
-      }
+        s[0][0] = -0.25 ; 
+        s[0][1] = -0.25 ;
+        s[0][2] = -0.25 ; 
+	s[1][0] =  0.25 ; 
+	s[1][1] = -0.25 ;
+	s[1][2] = -0.25 ; 
+	s[2][0] =  0.25 ; 
+	s[2][1] =  0.25 ;
+	s[2][2] = -0.25 ; 
+	s[3][0] = -0.25 ; 
+	s[3][1] =  0.25 ;
+	s[3][2] = -0.25 ; 
     }
     else{
-      s[0][0] = 0.25 * (-(1. - v) + v * w / (1. - w));
-      s[0][1] = 0.25 * (-(1. - u) + u * w / (1. - w));
-      s[0][2] = 0.25 * (-1.       + u * v / (1. - w) / (1. - w));
-      s[1][0] = 0.25 * ( (1. - v) + v * w / (1. - w));
-      s[1][1] = 0.25 * (-(1. + u) + u * w / (1. - w));
-      s[1][2] = 0.25 * (-1.       + u * v / (1. - w) / (1. - w));
-      s[2][0] = 0.25 * ( (1. + v) + v * w / (1. - w));
-      s[2][1] = 0.25 * ( (1. + u) + u * w / (1. - w));
-      s[2][2] = 0.25 * (-1.       + u * v / (1. - w) / (1. - w));
-      s[3][0] = 0.25 * (-(1. + v) + v * w / (1. - w));
-      s[3][1] = 0.25 * ( (1. - u) + u * w / (1. - w));
-      s[3][2] = 0.25 * (-1.       + u * v / (1. - w) / (1. - w));
+      s[0][0] = 0.25 * ( -(1. - v) + v * w / (1. - w)) ;
+      s[0][1] = 0.25 * ( -(1. - u) + u * w / (1. - w)) ;
+      s[0][2] = 0.25 * ( -1.     + u * v / (1. - w) + u * v * w / (1. - w) / (1. - w)) ; 
+      s[1][0] = 0.25 * (  (1. - v) - v * w / (1. - w)) ;
+      s[1][1] = 0.25 * ( -(1. + u) - u * w / (1. - w)) ;
+      s[1][2] = 0.25 * ( -1.     - u * v / (1. - w) - u * v * w / (1. - w) / (1. - w)) ; 
+      s[2][0] = 0.25 * (  (1. + v) + v * w / (1. - w)) ;
+      s[2][1] = 0.25 * (  (1. + u) + u * w / (1. - w)) ;
+      s[2][2] = 0.25 * ( -1.     + u * v / (1. - w) + u * v * w / (1. - w) / (1. - w)) ; 
+      s[3][0] = 0.25 * ( -(1. + v) - v * w / (1. - w)) ;
+      s[3][1] = 0.25 * (  (1. - u) - u * w / (1. - w)) ;
+      s[3][2] = 0.25 * ( -1.     - u * v / (1. - w) - u * v * w / (1. - w) / (1. - w)) ; 
     }
     s[4][0] = 0.; 
     s[4][1] = 0.;