diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 89c4f644be0325102e358e0813651e5549613318..82910aaa5f328ac75ea292a7cb9066db8a30f111 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -30,6 +30,7 @@
 #include "GmshMessage.h"
 #include "linearSystem.h"
 #include "Options.h"
+#include "polynomialBasis.h"
 
 #if defined(HAVE_OPENGL)
 #include "drawContext.h"
@@ -412,6 +413,7 @@ binding::binding()
   gmshOptions::registerBindings(this);
   Msg::registerBindings(this);
   linearSystem<double>::registerBindings(this);
+  polynomialBasis::registerBindings(this);
 #if defined(HAVE_SOLVER)
   function::registerBindings(this);
   linearSystemCSRGmm<double>::registerBindings(this);
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index efdd93e4b387aed3d1a934f46dd9ec9ce34ef301..36e1bb5cdf2152a688f916603e96c65ad7e36263 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -14,16 +14,16 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
   int order = (o == -1) ? getPolynomialOrder() : o;
   
   switch (order) {
-  case 1: return &polynomialBases::find(MSH_LIN_2);
-  case 2: return &polynomialBases::find(MSH_LIN_3);
-  case 3: return &polynomialBases::find(MSH_LIN_4);
-  case 4: return &polynomialBases::find(MSH_LIN_5);
-  case 5: return &polynomialBases::find(MSH_LIN_6);
-  case 6: return &polynomialBases::find(MSH_LIN_7);
-  case 7: return &polynomialBases::find(MSH_LIN_8);
-  case 8: return &polynomialBases::find(MSH_LIN_9);
-  case 9: return &polynomialBases::find(MSH_LIN_10);
-  case 10: return &polynomialBases::find(MSH_LIN_11);
+  case 1: return polynomialBases::find(MSH_LIN_2);
+  case 2: return polynomialBases::find(MSH_LIN_3);
+  case 3: return polynomialBases::find(MSH_LIN_4);
+  case 4: return polynomialBases::find(MSH_LIN_5);
+  case 5: return polynomialBases::find(MSH_LIN_6);
+  case 6: return polynomialBases::find(MSH_LIN_7);
+  case 7: return polynomialBases::find(MSH_LIN_8);
+  case 8: return polynomialBases::find(MSH_LIN_9);
+  case 9: return polynomialBases::find(MSH_LIN_10);
+  case 10: return polynomialBases::find(MSH_LIN_11);
   default: Msg::Error("Order %d line function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 75c6183938ead8a66bb8c972ccf0650073c850d3..81b6f2cb3dd0ee3033b712a06a93ac8e3a9da678 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -51,7 +51,7 @@ class MPoint : public MElement {
     s[0][0] = s[0][1] = s[0][2] = 0.;
   }
 
-  virtual const polynomialBasis* getFunctionSpace(int o) const { return &polynomialBases::find(MSH_PNT); }
+  virtual const polynomialBasis* getFunctionSpace(int o) const { return polynomialBases::find(MSH_PNT); }
   virtual bool isInside(double u, double v, double w)
   {
     return true;
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index e245eaeefc6b5819325a869ae2f006eb0d6f23c1..d6e99b70637c3b896acfbb0e3bb41378a9721bbf 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -38,15 +38,15 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
   
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_PRI_6);
-    case 2: return &polynomialBases::find(MSH_PRI_18);
+    case 1: return polynomialBases::find(MSH_PRI_6);
+    case 2: return polynomialBases::find(MSH_PRI_18);
     default: Msg::Error("Order %d prism function space not implemented", order);
     }
   }
   else { 
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_PRI_6);
-    case 2: return &polynomialBases::find(MSH_PRI_18);
+    case 1: return polynomialBases::find(MSH_PRI_6);
+    case 2: return polynomialBases::find(MSH_PRI_18);
     default: Msg::Error("Order %d prism function space not implemented", order);
     }
   }
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index e5d4afe433637afaf6653b827f0ff711610ae2cf..92eacb3a1bfd9464fad8520f21af9b0c420a51c9 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -24,29 +24,29 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-      case 1: return &polynomialBases::find(MSH_QUA_4);
-      case 2: return &polynomialBases::find(MSH_QUA_8);
-      case 3: return &polynomialBases::find(MSH_QUA_12);
-      case 4: return &polynomialBases::find(MSH_QUA_16I);
-      case 5: return &polynomialBases::find(MSH_QUA_20);
-      case 6: return &polynomialBases::find(MSH_QUA_24);
-      case 7: return &polynomialBases::find(MSH_QUA_28);
-      case 8: return &polynomialBases::find(MSH_QUA_32);
-      case 9: return &polynomialBases::find(MSH_QUA_36I);
-      case 10: return &polynomialBases::find(MSH_QUA_40);
+      case 1: return polynomialBases::find(MSH_QUA_4);
+      case 2: return polynomialBases::find(MSH_QUA_8);
+      case 3: return polynomialBases::find(MSH_QUA_12);
+      case 4: return polynomialBases::find(MSH_QUA_16I);
+      case 5: return polynomialBases::find(MSH_QUA_20);
+      case 6: return polynomialBases::find(MSH_QUA_24);
+      case 7: return polynomialBases::find(MSH_QUA_28);
+      case 8: return polynomialBases::find(MSH_QUA_32);
+      case 9: return polynomialBases::find(MSH_QUA_36I);
+      case 10: return polynomialBases::find(MSH_QUA_40);
     }
   }
   switch (order) {
-    case 1: return &polynomialBases::find(MSH_QUA_4);
-    case 2: return &polynomialBases::find(MSH_QUA_9);
-    case 3: return &polynomialBases::find(MSH_QUA_16);
-    case 4: return &polynomialBases::find(MSH_QUA_25);
-    case 5: return &polynomialBases::find(MSH_QUA_36);
-    case 6: return &polynomialBases::find(MSH_QUA_49);
-    case 7: return &polynomialBases::find(MSH_QUA_64);
-    case 8: return &polynomialBases::find(MSH_QUA_81);
-    case 9: return &polynomialBases::find(MSH_QUA_100);
-    case 10: return &polynomialBases::find(MSH_QUA_121);
+    case 1: return polynomialBases::find(MSH_QUA_4);
+    case 2: return polynomialBases::find(MSH_QUA_9);
+    case 3: return polynomialBases::find(MSH_QUA_16);
+    case 4: return polynomialBases::find(MSH_QUA_25);
+    case 5: return polynomialBases::find(MSH_QUA_36);
+    case 6: return polynomialBases::find(MSH_QUA_49);
+    case 7: return polynomialBases::find(MSH_QUA_64);
+    case 8: return polynomialBases::find(MSH_QUA_81);
+    case 9: return polynomialBases::find(MSH_QUA_100);
+    case 10: return polynomialBases::find(MSH_QUA_121);
     default: Msg::Error("Order %d triangle function space not implemented", order);
   }
   return 0;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index ca1cc08c550005fe23176be1e96f516697133403..2cdeb6704e7a1db921f88b915f6208f88b100fff 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -111,26 +111,26 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
   
   if ((nv == 0) && (o == -1)) {
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_TET_4);
-    case 2: return &polynomialBases::find(MSH_TET_10);
-    case 3: return &polynomialBases::find(MSH_TET_20);
-    case 4: return &polynomialBases::find(MSH_TET_34);
-    case 5: return &polynomialBases::find(MSH_TET_52);
+    case 1: return polynomialBases::find(MSH_TET_4);
+    case 2: return polynomialBases::find(MSH_TET_10);
+    case 3: return polynomialBases::find(MSH_TET_20);
+    case 4: return polynomialBases::find(MSH_TET_34);
+    case 5: return polynomialBases::find(MSH_TET_52);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
   else { 
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_TET_4);
-    case 2: return &polynomialBases::find(MSH_TET_10);
-    case 3: return &polynomialBases::find(MSH_TET_20);
-    case 4: return &polynomialBases::find(MSH_TET_35);
-    case 5: return &polynomialBases::find(MSH_TET_56);
-    case 6: return &polynomialBases::find(MSH_TET_84);
-    case 7: return &polynomialBases::find(MSH_TET_120);
-    case 8: return &polynomialBases::find(MSH_TET_165);
-    case 9: return &polynomialBases::find(MSH_TET_220);
-    case 10: return &polynomialBases::find(MSH_TET_286);
+    case 1: return polynomialBases::find(MSH_TET_4);
+    case 2: return polynomialBases::find(MSH_TET_10);
+    case 3: return polynomialBases::find(MSH_TET_20);
+    case 4: return polynomialBases::find(MSH_TET_35);
+    case 5: return polynomialBases::find(MSH_TET_56);
+    case 6: return polynomialBases::find(MSH_TET_84);
+    case 7: return polynomialBases::find(MSH_TET_120);
+    case 8: return polynomialBases::find(MSH_TET_165);
+    case 9: return polynomialBases::find(MSH_TET_220);
+    case 10: return polynomialBases::find(MSH_TET_286);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index ce58c06cec6cf440da51af506e9ab527fc8077e8..78fc812fa18cec6fcd7c15f3ddfde45f3c86f301 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -73,31 +73,31 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
 
   if ((nf == 0) && (o == -1)) {
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_TRI_3);
-    case 2: return &polynomialBases::find(MSH_TRI_6);
-    case 3: return &polynomialBases::find(MSH_TRI_9);
-    case 4: return &polynomialBases::find(MSH_TRI_12);
-    case 5: return &polynomialBases::find(MSH_TRI_15I);
-    case 6: return &polynomialBases::find(MSH_TRI_18);
-    case 7: return &polynomialBases::find(MSH_TRI_21I);
-    case 8: return &polynomialBases::find(MSH_TRI_24);
-    case 9: return &polynomialBases::find(MSH_TRI_27);
-    case 10: return &polynomialBases::find(MSH_TRI_30);
+    case 1: return polynomialBases::find(MSH_TRI_3);
+    case 2: return polynomialBases::find(MSH_TRI_6);
+    case 3: return polynomialBases::find(MSH_TRI_9);
+    case 4: return polynomialBases::find(MSH_TRI_12);
+    case 5: return polynomialBases::find(MSH_TRI_15I);
+    case 6: return polynomialBases::find(MSH_TRI_18);
+    case 7: return polynomialBases::find(MSH_TRI_21I);
+    case 8: return polynomialBases::find(MSH_TRI_24);
+    case 9: return polynomialBases::find(MSH_TRI_27);
+    case 10: return polynomialBases::find(MSH_TRI_30);
     default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
     }
   }
   else { 
     switch (order) {
-    case 1: return &polynomialBases::find(MSH_TRI_3);
-    case 2: return &polynomialBases::find(MSH_TRI_6);
-    case 3: return &polynomialBases::find(MSH_TRI_10);
-    case 4: return &polynomialBases::find(MSH_TRI_15);
-    case 5: return &polynomialBases::find(MSH_TRI_21);
-    case 6: return &polynomialBases::find(MSH_TRI_28);
-    case 7: return &polynomialBases::find(MSH_TRI_36);
-    case 8: return &polynomialBases::find(MSH_TRI_45);
-    case 9: return &polynomialBases::find(MSH_TRI_55);
-    case 10: return &polynomialBases::find(MSH_TRI_66);
+    case 1: return polynomialBases::find(MSH_TRI_3);
+    case 2: return polynomialBases::find(MSH_TRI_6);
+    case 3: return polynomialBases::find(MSH_TRI_10);
+    case 4: return polynomialBases::find(MSH_TRI_15);
+    case 5: return polynomialBases::find(MSH_TRI_21);
+    case 6: return polynomialBases::find(MSH_TRI_28);
+    case 7: return polynomialBases::find(MSH_TRI_36);
+    case 8: return polynomialBases::find(MSH_TRI_45);
+    case 9: return polynomialBases::find(MSH_TRI_55);
+    case 10: return polynomialBases::find(MSH_TRI_66);
     default: Msg::Error("Order %d triangle function space not implemented", order);
     }
   }
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index bdae7c5046cb87de2cc9a7750625545a499ad1a8..b6c5c828520cffa8dd9bca9d1133ab887604c22a 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -515,15 +515,15 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
   
   switch (nPts){
   case 2 :
-    points = polynomialBases::find(MSH_TRI_10).points;
+    points = polynomialBases::find(MSH_TRI_10)->points;
     start = 9;
     break;
   case 3 :
-    points = polynomialBases::find(MSH_TRI_15).points;
+    points = polynomialBases::find(MSH_TRI_15)->points;
     start = 12;
     break;
   case 4 :
-    points = polynomialBases::find(MSH_TRI_21).points;
+    points = polynomialBases::find(MSH_TRI_21)->points;
     start = 15;
     break;
   default :  
@@ -608,11 +608,11 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
 
   switch (nPts){
   case 3 :
-    points = polynomialBases::find(MSH_TET_35).points;
+    points = polynomialBases::find(MSH_TET_35)->points;
     start = 34;
     break;
   case 4 :
-    points = polynomialBases::find(MSH_TET_56).points;
+    points = polynomialBases::find(MSH_TET_56)->points;
     start = 52;
     break;
   default :  
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 79a94761f0024b9ef9f51fedd8b2699a52b7d21b..474f4f025098e6db7fc1d0dbe20fbddeb56f2d51 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -906,10 +906,10 @@ static void generate1dVertexClosure(polynomialBasis::clCont &closure)
 }
 std::map<int, polynomialBasis> polynomialBases::fs;
 
-const polynomialBasis &polynomialBases::find(int tag)
+const polynomialBasis *polynomialBases::find(int tag)
 {
   std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
-  if (it != fs.end())     return it->second;
+  if (it != fs.end())     return &it->second;
   polynomialBasis F;
   F.numFaces = -1;
 
@@ -1305,7 +1305,7 @@ const polynomialBasis &polynomialBases::find(int tag)
 //   }
 
   fs.insert(std::make_pair(tag, F));
-  return fs[tag];
+  return &fs[tag];
 }
 
 
@@ -1317,8 +1317,8 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
   std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
   if (it != injector.end()) return it->second;
 
-  const polynomialBasis& fs1 = find(tag1);
-  const polynomialBasis& fs2 = find(tag2);
+  const polynomialBasis& fs1 = *find(tag1);
+  const polynomialBasis& fs2 = *find(tag2);
 
   fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
 
@@ -1332,3 +1332,15 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
   injector.insert(std::make_pair(key, inj));
   return injector[key];
 }
+
+#include "Bindings.h"
+void polynomialBasis::registerBindings(binding *b) {
+  classBinding *cb = b->addClass<polynomialBasis>("polynomialBasis");
+  cb->setDescription("polynomial shape functions for elements");
+  methodBinding *mb = cb->addMethod("f",(void (polynomialBasis::*)(fullMatrix<double>&, fullMatrix<double>&))&polynomialBasis::f);
+  mb->setDescription("evaluate the shape functions");
+  mb->setArgNames("nodes","values",NULL);
+  mb = cb->addMethod("find",&polynomialBases::find);
+  mb->setDescription("return the polynomial basis corresponding to an element type");
+  mb->setArgNames("elementType",NULL);
+}
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 80cece67066ae08efb2630a70f1d26d208fbcc21..e0f38baa2161b0929dbb12243bcacb0f0feb2d8e 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -62,6 +62,7 @@ inline double pow_int (const double &a , const int &n) {
   }
 } 
 
+class binding;
 
 // presently those function spaces are only for simplices and quads;
 // should be extended to other elements like hexes
@@ -102,6 +103,17 @@ class polynomialBasis
       }
     }
   }
+  // I would favour this interface that allows optimizations (not implemented) and is easier to bind
+  inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) {
+    double p[256];
+    sf.resize (coefficients.size1(), coord.size1());
+    for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
+      evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
+      for (int i = 0; i < coefficients.size1(); i++)
+        for (int j = 0; j < coefficients.size2(); j++)
+          sf(i,iPoint) += coefficients(i, j) * p[j];
+    }
+  }
   inline void df(double u, double v, double w, double grads[][3]) const
   {
     switch (monomials.size2()) {
@@ -242,6 +254,7 @@ class polynomialBasis
       break;
     }
   }
+  static void registerBindings(binding *b);
 };
 
 class polynomialBases
@@ -250,7 +263,7 @@ class polynomialBases
   static std::map<int, polynomialBasis> fs;
   static std::map<std::pair<int, int>, fullMatrix<double> > injector;
  public :
-  static const polynomialBasis &find(int);
+  static const polynomialBasis *find(int);
   static const fullMatrix<double> &findInjector(int, int);
 };
 
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 21da15f8795c07d10f6c329e7dc10e7817a4c40e..2a1dc767dbd46638527ad0d6a8d55e4325cf54ac 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -120,6 +120,27 @@ bool PViewData::combineSpace(nameData &nd)
   return false; 
 }
 
+double PViewData::getValueBinding(int step, int ent, int ele, int nod, int comp)
+{
+  double val;
+  getValue(step,ent,ele,nod,comp,val);
+  return val;
+}
+
+void PViewData::getAllValuesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m) {
+  int nNodes = getNumNodes(step,ent,ele);
+  int nComponents = getNumComponents(step,ent,ele);
+  for (int i=0; i<nNodes; i++)
+    for (int j=0; j<nComponents; j++)
+       getValue(step,ent,ele,i,j,m(i,j));
+}
+
+void PViewData::getAllNodesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m) {
+  int nNodes = getNumNodes(step,ent,ele);
+  for (int i=0; i<nNodes; i++)
+    getNode(step,ent,ele,i,m(i,0),m(i,1),m(i,2));
+}
+
 #include "Bindings.h"
 void PViewData::registerBindings(binding *b) {
   classBinding *cb = b->addClass<PViewData>("PViewData");
@@ -127,5 +148,35 @@ void PViewData::registerBindings(binding *b) {
   methodBinding *cm;
   cm = cb->addMethod("getNumEntities",&PViewData::getNumEntities);
   cm->setArgNames("step",NULL);
-  cm->setDescription("return the number of entities for a given time step");
+  cm->setDescription("return the number of entities for a given time step (-1 for default)");
+  cm = cb->addMethod("getNumElements",&PViewData::getNumElements);
+  cm->setArgNames("step","entity",NULL);
+  cm->setDescription("return the number of entities for a given time step (-1 for default) for a given entity (-1 for all)");
+  cm = cb->addMethod("getNumTriangles",&PViewData::getNumTriangles);
+  cm->setArgNames("step",NULL);
+  cm->setDescription("return the number of triangles for a given time step (-1 for default)");
+
+  cm = cb->addMethod("getNumNodes",&PViewData::getNumNodes);
+  cm->setArgNames("step","entity","element",NULL);
+  cm->setDescription("return the number of nodes of one element of an entity of a time step (-1 for default time step)");
+
+  cm = cb->addMethod("getNumValues",&PViewData::getNumValues);
+  cm->setArgNames("step","entity","element",NULL);
+  cm->setDescription("return the number of values of one element of an entity of a time step (-1 for default time step)");
+
+  cm = cb->addMethod("getNumComponents",&PViewData::getNumComponents);
+  cm->setArgNames("step","entity","element",NULL);
+  cm->setDescription("return the number of components of one element of an entity of a time step (-1 for default time step)");
+
+  cm = cb->addMethod("getValue",&PViewData::getValueBinding);
+  cm->setArgNames("step","entity","element","node","component",NULL);
+  cm->setDescription("return one of the values");
+
+  cm = cb->addMethod("getAllValuesForElement",&PViewData::getAllValuesForElementBinding);
+  cm->setArgNames("step","entity","element","values",NULL);
+  cm->setDescription("fill a fullMatrix with all values from the elements. The fullMatrix should have the size nbNodes x nbComponents.");
+
+  cm = cb->addMethod("getAllNodesForElement",&PViewData::getAllNodesForElementBinding);
+  cm->setArgNames("step","entity","element","coordinates",NULL);
+  cm->setDescription("fill a fullMatrix with all coordinates from the nodes of the elements. The fullMatrix should have the size nbNodes x 3");
 }
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 24879b7e10529c6c611d812c540f3deb85709f22..5b0e24808f8934f053c60c3f2c80e320a2708769 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -132,6 +132,10 @@ class PViewData {
   virtual void getValue(int step, int ent, int ele, int nod, int comp, double &val){}
   virtual void setValue(int step, int ent, int ele, int nod, int comp, double val);
 
+  double getValueBinding(int step, int ent, int ele, int nod, int comp);
+  void getAllValuesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m);
+  void getAllNodesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m);
+
   // return a scalar value (same as value for scalars, norm for
   // vectors, etc.) associated with the node-th node from the ele-th
   // element in the ent-th entity
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 8c39b8d7f922382385455341a8c2a3d49315c1ea..53a68f5ca2f3a7fef2cfd7e94f052c7e42bd39dd 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -784,7 +784,7 @@ void PViewDataList::setOrder2(int type)
   case TYPE_PRI: typeMSH = MSH_PRI_18; break;
   case TYPE_PYR: typeMSH = MSH_PYR_14; break;
   }
-  const polynomialBasis *fs = &polynomialBases::find(typeMSH);
+  const polynomialBasis *fs = polynomialBases::find(typeMSH);
   if(!fs){
     Msg::Error("Could not find function space for element type %d", typeMSH);
     return;
diff --git a/benchmarks/elasticity/conge.lua b/benchmarks/elasticity/conge.lua
index 041ae640b0ef5f31d30428780caa6852cf41e19b..1319eb85f514103c7521d2a9fd251f115192d146 100644
--- a/benchmarks/elasticity/conge.lua
+++ b/benchmarks/elasticity/conge.lua
@@ -15,13 +15,37 @@ sys = linearSystemCSRTaucs()
 e:assemble(sys)
 sys:systemSolve()
 view = e:buildVonMisesView("vonMises")
-print(view)
-view:write("vonMises.msh",5,false)
+
 view = e:buildDisplacementView("displacement")
 view:write("displacement.msh",5,false)
-view = e:buildLagrangeMultiplierView("lagrangeMultiplier")
-view:write("lagrangeMultiplier.msh",5,false)
-view = e:buildElasticEnergyView("elasticEnergy")
-view:write("elasticEnergy.msh",5,false)
+
 viewData = view:getData()
-print(viewData:getNumEntities(-1))
+nEntities = viewData:getNumEntities(-1)
+v = fullMatrix(3,3)
+nodes = fullMatrix(3,3)
+for i=0,nEntities-1 do 
+  nElements = viewData:getNumElements(-1,i)
+  for j=0,nElements-1 do
+    if(viewData:getNumNodes(0,i,j)==3) then
+      viewData:getAllValuesForElement(0,i,j,v)
+      viewData:getAllNodesForElement(0,i,j,nodes)
+    end
+  end
+end
+
+basis = polynomialBasis.find(2)
+points = fullMatrix(3,3)
+points:set(0,0,0);
+points:set(0,1,0);
+points:set(0,2,0);
+
+points:set(1,0,1);
+points:set(1,1,0);
+points:set(1,2,0);
+
+points:set(2,0,0);
+points:set(2,1,1);
+points:set(2,2,0);
+
+basis:f(points,v)
+v:print("f")