diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 82910aaa5f328ac75ea292a7cb9066db8a30f111..f2ccb95d3563632f2835bf36c437a9ab08fe3322 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -31,6 +31,8 @@
 #include "linearSystem.h"
 #include "Options.h"
 #include "polynomialBasis.h"
+#include "Gauss.h"
+#include "PViewFactory.h"
 
 #if defined(HAVE_OPENGL)
 #include "drawContext.h"
@@ -414,6 +416,8 @@ binding::binding()
   Msg::registerBindings(this);
   linearSystem<double>::registerBindings(this);
   polynomialBasis::registerBindings(this);
+  gaussIntegration::registerBindings(this);
+  PViewFactory::registerBindings(this);
 #if defined(HAVE_SOLVER)
   function::registerBindings(this);
   linearSystemCSRGmm<double>::registerBindings(this);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 3d2cb90f24431058ab1960145c6993b461176904..424d24a0ed2bac36e4c6b569518a82c7588e0a11 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -258,6 +258,12 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
   return _computeDeterminantAndRegularize(this, jac);
 }
 
+//binded
+double MElement::getJacobianDeterminant(double u, double v, double w) {
+  double jac[3][3];
+  return getJacobian(u,v,w,jac);
+}
+
 
 double MElement::getJacobian(double gsf[][3], double jac[3][3])
 {
@@ -1030,4 +1036,7 @@ void MElement::registerBindings(binding *b)
   cm->setDescription("return the polynomial order the element");
   cm = cb->addMethod("getDim", &MElement::getDim);
   cm->setDescription("return the geometrical dimension of the element");
+  cm = cb->addMethod("getJacobianDeterminant", &MElement::getJacobianDeterminant);
+  cm->setDescription("return the jacobian of the determinant of the transformation");
+  cm->setArgNames("u","v","w",NULL);
 }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 1d392e94d8b8cf1113594912d0aa6a7aed5ec82a..405388a24e54300ba16d4c2057ecd586dd68ada4 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -220,6 +220,8 @@ class MElement
   double getJacobian(double gsf[][3], double jac[3][3]);
   double getJacobian(double u, double v, double w, double jac[3][3]);
   double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
+  //bindings : double[3][3] is not easy to bind, we could use fullMatrix instead
+  double getJacobianDeterminant(double u, double v, double w);
 
   // get the point in cartesian coordinates corresponding to the point
   // (u,v,w) in parametric coordinates
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index 68dc509f6b05d0cc327bf978141880e5535612e2..4d7f5c9233834347161143f8b3c36cea8b96ab51 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -14,6 +14,7 @@ set(SRC
     GaussQuadratureHex.cpp
     GaussQuadraturePri.cpp
     GaussLegendreSimplex.cpp
+    Gauss.cpp
   robustPredicates.cpp
   mathEvaluator.cpp
   Iso.cpp
diff --git a/Numeric/Gauss.cpp b/Numeric/Gauss.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f599c162984b3f7974f7a734ea19d196eb5ed026
--- /dev/null
+++ b/Numeric/Gauss.cpp
@@ -0,0 +1,54 @@
+#include "Gauss.h"
+#include "Bindings.h"
+static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix, fullMatrix<double> &wMatrix) {
+  pMatrix.resize(npts,3);
+  wMatrix.resize(npts,1);
+  for (int i=0;i<npts;i++) {
+    pMatrix(i,0) = pts[i].pt[0];
+    pMatrix(i,1) = pts[i].pt[1];
+    pMatrix(i,2) = pts[i].pt[2];
+    wMatrix(i,0) = pts[i].weight;
+  }
+}
+void gaussIntegration::registerBindings(binding *b)
+{
+  classBinding *cb = b->addClass<gaussIntegration>("gaussIntegration");
+  cb->setDescription("Gauss integration rules (points+weights)");
+  methodBinding *mb = cb->addMethod("getTriangle", &gaussIntegration::getTriangle);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a triangle");
+  mb->setArgNames("order","points","weights",NULL);
+  mb = cb->addMethod("getLine", &gaussIntegration::getLine);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a line");
+  mb->setArgNames("order","points","weights",NULL);
+  mb = cb->addMethod("getQuad", &gaussIntegration::getQuad);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a quad");
+  mb->setArgNames("order","points","weights",NULL);
+  mb = cb->addMethod("getTetrahedron", &gaussIntegration::getTetrahedron);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a tetrahedron");
+  mb->setArgNames("order","points","weights",NULL);
+  mb = cb->addMethod("getHexahedron", &gaussIntegration::getHexahedron);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a hexahedron");
+  mb->setArgNames("order","points","weights",NULL);
+  mb = cb->addMethod("getPrism", &gaussIntegration::getPrism);
+  mb->setDescription("get the integration rule that integrate exactly a polynom of given order on a prism");
+  mb->setArgNames("order","points","weights",NULL);
+}
+
+void gaussIntegration::getTriangle(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQTPts(order),getGQTPts(order),pts,weights);
+}
+void gaussIntegration::getLine(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQLPts(order),getGQLPts(order),pts,weights);
+}
+void gaussIntegration::getQuad(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQQPts(order),getGQQPts(order),pts,weights);
+}
+void gaussIntegration::getTetrahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQTetPts(order),getGQTetPts(order),pts,weights);
+}
+void gaussIntegration::getHexahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQHPts(order),getGQHPts(order),pts,weights);
+}
+void gaussIntegration::getPrism(int order, fullMatrix<double> &pts, fullMatrix<double> &weights) {
+  pts2fullMatrix(getNGQPriPts(order),getGQPriPts(order),pts,weights);
+}
diff --git a/Numeric/Gauss.h b/Numeric/Gauss.h
index d3f8791b59608c846a6054a7640887ed00008190..6be921404359def80fcfc4479cfc6ef8ebcefd2e 100644
--- a/Numeric/Gauss.h
+++ b/Numeric/Gauss.h
@@ -5,6 +5,8 @@
 
 #ifndef _GAUSS_H_
 #define _GAUSS_H_
+#include "fullMatrix.h"
+
 
 struct IntPt{
   double pt[3];
@@ -34,4 +36,16 @@ IntPt *getGQPriPts(int order);
 int getNGQHPts(int order);
 IntPt *getGQHPts(int order);
 
+//For now this class is only for bindings but maybe the interface is cleaner (it does not add new types) and it can replace the other interface
+class gaussIntegration {
+  public:
+  static void getTriangle(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void getLine(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void getQuad(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void getTetrahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void getHexahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void getPrism(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+  static void registerBindings(binding *b);
+};
+
 #endif
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index e0f38baa2161b0929dbb12243bcacb0f0feb2d8e..df6ee978196e83b76f6016e1cb01a71d626493e1 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -106,12 +106,12 @@ 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());
+    sf.resize (coord.size1(), coefficients.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];
+          sf(iPoint,i) += coefficients(i, j) * p[j];
     }
   }
   inline void df(double u, double v, double w, double grads[][3]) const
diff --git a/Post/CMakeLists.txt b/Post/CMakeLists.txt
index 6942bee19d6ba5500cfe15b4c5c5bff8de9bec43..cad2a94397a999262f178a52dd53495dfd380a39 100644
--- a/Post/CMakeLists.txt
+++ b/Post/CMakeLists.txt
@@ -9,6 +9,7 @@ set(SRC
       PViewDataList.cpp PViewDataListIO.cpp
       PViewDataGModel.cpp PViewDataGModelIO.cpp
     PViewOptions.cpp
+    PViewFactory.cpp
   adaptiveData.cpp shapeFunctions.cpp
   OctreePost.cpp
   ColorTable.cpp
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 5b0e24808f8934f053c60c3f2c80e320a2708769..125a00608964ef612307bfca1b7722d03ce1b96c 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -16,6 +16,7 @@
 
 class adaptiveData;
 class GModel;
+class GEntity;
 class nameData;
 class binding;
 
@@ -218,6 +219,8 @@ class PViewData {
                         bool append=false);
   virtual bool writeMSH(std::string fileName, bool binary=false);
   virtual bool writeMED(std::string fileName);
+  //
+  virtual GEntity *getEntity (int step, int entity) {return NULL;}
   static void registerBindings(binding *b);
 };
 
diff --git a/Post/PViewFactory.cpp b/Post/PViewFactory.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..df68730d59e1c58e02483f2362d655d1c9c38f3e
--- /dev/null
+++ b/Post/PViewFactory.cpp
@@ -0,0 +1,41 @@
+#include "PViewFactory.h"
+#include "GModel.h"
+#include "fullMatrix.h"
+#include "PView.h"
+#include <vector>
+#include "Bindings.h"
+
+PViewFactory::PViewFactory (std::string name, std::string type, GModel *model, int timeStep, int dim):_model(model),_name(name),_type(type),_timeStep(timeStep), _dim(dim)
+{
+}
+
+void PViewFactory::setEntry (int id, const fullMatrix<double> &val) 
+{
+  std::vector<double> &vv = _values[id];
+  vv.resize(val.size1()*val.size2());
+  int k=0;
+  for (int i=0;i<val.size1(); i++) {
+    for (int j=0;j<val.size2(); j++) {
+      vv[k++] = val(i,j);
+    }
+  }
+}
+
+PView *PViewFactory::createView ()
+{
+  return new PView(_name,_type,_model,_values,_timeStep,_dim);
+}
+
+void PViewFactory::registerBindings (class binding *b)
+{
+  classBinding *cb = b->addClass<PViewFactory>("PViewFactory");
+  cb->setDescription(" ");
+  methodBinding *mb = cb->addMethod("setEntry", &PViewFactory::setEntry);
+  mb->setDescription(" ");
+  mb->setArgNames("elementId", "values", NULL);
+  mb = cb->setConstructor<PViewFactory, std::string, std::string, GModel*,int,int>();
+  mb->setDescription(" ");
+  mb->setArgNames("name", "type", "model","timeStep", "dimension", NULL);
+  mb = cb->addMethod("createView",&PViewFactory::createView);
+  mb->setDescription(" ");
+}
diff --git a/Post/PViewFactory.h b/Post/PViewFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..3502d27e4ef1b778f4f0ba211bdaaf04187e115e
--- /dev/null
+++ b/Post/PViewFactory.h
@@ -0,0 +1,26 @@
+#ifndef _PVIEW_FACTORY_H_
+#define _PVIEW_FACTORY_H_
+//quick hack to have something that we can bind for the summer school...
+//this class has probably to be removed or rewritten
+
+#include<map>
+#include<vector>
+#include<string>
+class GModel;
+class PView;
+template <class t>
+class fullMatrix;
+class binding;
+
+class PViewFactory {
+  std::map<int,std::vector<double> > _values;
+  std::string _name, _type;
+  int _dim,_timeStep;
+  GModel *_model;
+  public:
+  PViewFactory(std::string name, std::string type, GModel *model, int timeStep, int dim);
+  void setEntry(int id, const fullMatrix<double> &val);
+  PView *createView();
+  static void registerBindings(binding *);
+};
+#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index bc99f8d333eb3188b14555e1d139c4a3c2aef844..04794331916c2b94115192146e5e85205ecb3272 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -24,7 +24,6 @@
 #include "PViewData.h"
 #endif
 
-/********************* deprecated api ? ****************************************/
 static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
 {
   int *startIndex;
@@ -97,6 +96,21 @@ void elasticitySolver::solve()
   printf("elastic energy=%f\n",energ);
 }
 
+void elasticitySolver::postSolve()
+{
+  GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+
+  double energ=0;
+  for (unsigned int i = 0; i < elasticFields.size(); i++)
+  {
+    SolverField<SVector3> Field(pAssembler, LagSpace);
+    IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+    BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+    Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+  }
+  printf("elastic energy=%f\n",energ);
+}
+
 void elasticitySolver::readInputFile(const std::string &fn)
 {
   FILE *f = fopen(fn.c_str(), "r");
@@ -247,7 +261,6 @@ void elasticitySolver::readInputFile(const std::string &fn)
   fclose(f);
 }
 
-/********************* end deprecated api ****************************************/
 
 #if defined (HAVE_LUA)
 
@@ -414,6 +427,30 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
 
 
 #if defined(HAVE_POST)
+static void deformation(dofManager<double> *a, MElement *e,
+                       double u, double v, double w, int _tag, double *eps){
+  double valx[256];
+  double valy[256];
+  double valz[256];
+  for (int k = 0; k < e->getNumVertices(); k++){
+    a->getDofValue(e->getVertex(k), 0, _tag, valx[k]);
+    a->getDofValue(e->getVertex(k), 1, _tag, valy[k]);
+    a->getDofValue(e->getVertex(k), 2, _tag, valz[k]);
+  }
+  double gradux[3];
+  double graduy[3];
+  double graduz[3];
+  e->interpolateGrad(valx, u, v, w, gradux);
+  e->interpolateGrad(valy, u, v, w, graduy);
+  e->interpolateGrad(valz, u, v, w, graduz);
+
+  eps[0] = gradux[0];
+  eps[1] = graduy[1];
+  eps[2] = graduz[2];
+  eps[3] = 0.5 * (gradux[1] + graduy[0]);
+  eps[4] = 0.5 * (gradux[2] + graduz[0]);
+  eps[5] = 0.5 * (graduy[2] + graduz[1]);
+}
 
 static double vonMises(dofManager<double> *a, MElement *e,
                        double u, double v, double w,
@@ -455,6 +492,10 @@ static double vonMises(dofManager<double> *a, MElement *e,
   return ComputeVonMises(s);
 }
 
+void elasticitySolver::getSolutionOnElement (MElement *el, fullMatrix<double> &sol) {
+  
+}
+
 PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
 {
   std::cout <<  "build Displacement View"<< std::endl;
@@ -624,6 +665,9 @@ void elasticitySolverRegisterBindings(binding *b)
   cm->setDescription ("reads an input file");
   cm->setArgNames("fileName",NULL);
 
+  cm = cb->addMethod("postSolve", &elasticitySolver::postSolve);
+  cm->setDescription (" ");
+
   cm = cb->addMethod("assemble", &elasticitySolver::assemble);
   cm->setDescription ("assembles the problem");
   cm->setArgNames ("linearSystem",NULL);
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 767ef8fcfa7c216ddef7470cef1d8cc485c85b44..58331aa5aaeb1465594b390f5ecec968e5c64551 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -100,6 +100,8 @@ class elasticitySolver
   void read(const std::string s) {readInputFile(s.c_str());}
   virtual void setMesh(const std::string &meshFileName);
   void solve();
+  void postSolve();
+  void getSolutionOnElement(MElement *el, fullMatrix<double> &sol);
   virtual PView *buildDisplacementView(const std::string postFileName);
   virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
   virtual PView *buildElasticEnergyView(const std::string postFileName);
diff --git a/benchmarks/elasticity/conge.lua b/benchmarks/elasticity/conge.lua
index 1319eb85f514103c7521d2a9fd251f115192d146..ee8860733889f1b9ae6c418719a50b7902707235 100644
--- a/benchmarks/elasticity/conge.lua
+++ b/benchmarks/elasticity/conge.lua
@@ -2,6 +2,37 @@ function neumannCondition (x,y,z)
   return {-100e6,0,0}
 end
 
+
+function errorOnView (view, dim)
+  basis = polynomialBasis.find(2)
+  gaussPoints = fullMatrix(1,1)
+  gaussWeights = fullMatrix(1,1)
+  gaussIntegration.getTriangle(3,gaussPoints,gaussWeights)
+
+  nnodes = 3 -- only triangles
+  viewData = view:getData()
+  values = fullMatrix(nnodes,dim)
+  nodes = fullMatrix(nnodes,3)
+  gaussXYZ = fullMatrix(gaussPoints:size1(),3)
+  gaussValues = fullMatrix(gaussPoints:size1(),dim)
+  nEntities = viewData:getNumEntities(-1)
+  psi = fullMatrix(gaussPoints:size1(),3)
+  basis:f(gaussPoints,psi)
+  for i=0,nEntities-1 do 
+    nElements = viewData:getNumElements(-1,i)
+    for j=0,nElements-1 do
+      if(viewData:getNumNodes(0,i,j)==nnodes) then
+        viewData:getAllValuesForElement(0,i,j,values)
+        viewData:getAllNodesForElement(0,i,j,nodes)
+      end
+      gaussXYZ:gemm(psi,nodes,1,0)
+      gaussValues:gemm(psi,values,1,0)
+    end
+  end
+end
+
+
+
 m = GModel()
 m:load("conge.geo")
 m:load("conge.msh")
@@ -19,33 +50,4 @@ view = e:buildVonMisesView("vonMises")
 view = e:buildDisplacementView("displacement")
 view:write("displacement.msh",5,false)
 
-viewData = view:getData()
-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")
+errorOnView(view,3)