diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e4b8571cf4a3291bbbedaf7e03d447ed72c69da8..40ca858f718aee2698bb0591896b59441ad5d9a0 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1207,6 +1207,26 @@ const fullMatrix<double> &MElement::getGradShapeFunctionsAtIntegrationPoints (in
   }
   return mat;
 }
+
+const fullMatrix<double> &MElement::getShapeFunctionsAtIntegrationPoints (int integrationOrder, int functionSpaceOrder) {
+  static std::map <std::pair<int,int>, fullMatrix<double> > F;
+  const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
+  fullMatrix<double> &mat = F [std::make_pair(fs->type, integrationOrder)];
+  if (mat.size1()!=0) return mat;
+  int npts;
+  IntPt *pts;
+  getIntegrationPoints (integrationOrder, &npts, &pts);
+  mat.resize (fs->points.size1(), npts*3);
+  double f[512];
+  for (int i = 0; i < npts; i++) {
+    fs->f (pts[i].pt[0], pts[i].pt[1], pts[i].pt[2], f);
+    for (int j = 0; j < fs->points.size1(); j++) {
+      mat (j, i) = f[j];
+    }
+  }
+  return mat;
+}
+
 const fullMatrix<double> &MElement::getGradShapeFunctionsAtNodes (int functionSpaceOrder) {
   static std::map <int, fullMatrix<double> > DF;
   const polynomialBasis *fs = getFunctionSpace (functionSpaceOrder);
diff --git a/Geo/MElement.h b/Geo/MElement.h
index cc10d2cc8c7df21483a1d35c1f6ce7098b22eb11..22133a09f4126df8fc54324d0941f2484d1114e4 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -239,6 +239,8 @@ class MElement
                                      int order=-1);
   const fullMatrix<double> &getGradShapeFunctionsAtIntegrationPoints
     (int integrationOrder, int functionSpaceOrder=-1);
+  const fullMatrix<double> &getShapeFunctionsAtIntegrationPoints
+    (int integrationOrder, int functionSpaceOrder=-1);
   const fullMatrix<double> &getGradShapeFunctionsAtNodes (int functionSpaceOrder=-1);
 
   // return the Jacobian of the element evaluated at point (u,v,w) in
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 85a250c70e1f5a81c317175694e25c742a493345..2c9e7a61222c35d50ef5c29dd62f6518af28a9b5 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -420,7 +420,7 @@ class fullMatrix
   }
 #endif
   ;
-  inline fullMatrix<scalar> transpose()
+  inline fullMatrix<scalar> transpose() const
   {
     fullMatrix<scalar> T(size2(), size1());
     for(int i = 0; i < size1(); i++)