diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index a6c7d63d7ca09e71dc71041154e6efff87ab883e..4de8f2ce404986639338b9031713d860fab6f869 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -17,6 +17,9 @@ set(SRC
   CurlBasis.cpp
   DivBasis.cpp
 
+  EvaluatedBasis.cpp
+  EvaluatedBasisScalar.cpp
+  EvaluatedBasisVector.cpp
 
   LineNodeBasis.cpp
   LineEdgeBasis.cpp
diff --git a/FunctionSpace/EvaluatedBasis.cpp b/FunctionSpace/EvaluatedBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b06376ebb6af29eb34bb6df8c47c7db0a43c1d2
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasis.cpp
@@ -0,0 +1,7 @@
+#include "EvaluatedBasis.h"
+
+EvaluatedBasis::EvaluatedBasis(void){
+}
+
+EvaluatedBasis::~EvaluatedBasis(void){
+}
diff --git a/FunctionSpace/EvaluatedBasis.h b/FunctionSpace/EvaluatedBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..0ef09504102d050344ef6df0c811cab49db1fbc8
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasis.h
@@ -0,0 +1,101 @@
+#ifndef _EVALUATEDBASIS_H_
+#define _EVALUATEDBASIS_H_
+
+class EvaluatedBasis{
+ protected:
+  bool scalar;
+
+  unsigned int nVertex;
+  unsigned int nEdge;
+  unsigned int nFace;
+  unsigned int nCell;
+
+  unsigned int nEdgeClosure;
+  unsigned int nFaceClosure;
+
+  unsigned int nPoint;
+
+ public:
+  //! Deletes this EvaluatedBasis
+  //!
+  virtual ~EvaluatedBasis(void);
+
+  //! @return Returns:
+  //! @li @c true, if the evaluated
+  //! values are @em scalar  
+  //! @li @c false, otherwise  
+  bool isScalar(void) const;
+
+  //! @return Returns the number of @em Vertex
+  //! @em Based functions of this EvaluatedBasis
+  unsigned int getNVertexBased(void) const;
+
+  //! @return Returns the number of @em Edge
+  //! @em Based functions of this EvaluatedBasis
+  unsigned int getNEdgeBased(void) const;
+
+  //! @return Returns the number of @em Face
+  //! @em Based functions of this EvaluatedBasis
+  unsigned int getNFaceBased(void) const;
+
+  //! @return Returns the number of @em Cell
+  //! @em Based functions of this EvaluatedBasis
+  unsigned int getNCellBased(void) const;
+
+  //! @return Returns the number of @em edge
+  //! @em closures for this EvaluatedBasis
+  unsigned int getNEdgeClosure(void) const;
+
+  //! @return Returns the number of @em face
+  //! @em closures for this EvaluatedBasis
+  unsigned int getNFaceClosure(void) const;
+
+  //! @return Returns the number of 
+  //! evaluation Points
+  unsigned int getNPoint(void) const;
+
+ protected:
+  //! @internal
+  //! Instantiate a new EvaluatedBasis
+  //!
+  //! @endinternal
+  EvaluatedBasis(void);
+};
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline bool EvaluatedBasis::isScalar(void) const{
+  return scalar;
+}
+
+inline unsigned int EvaluatedBasis::getNVertexBased(void) const{
+  return nVertex;
+}
+
+inline unsigned int EvaluatedBasis::getNEdgeBased(void) const{
+  return nEdge;
+}
+
+inline unsigned int EvaluatedBasis::getNFaceBased(void) const{
+  return nFace;
+}
+
+inline unsigned int EvaluatedBasis::getNCellBased(void) const{
+  return nCell;
+}
+
+inline unsigned int EvaluatedBasis::getNEdgeClosure(void) const{
+  return nEdgeClosure;
+}
+
+inline unsigned int EvaluatedBasis::getNFaceClosure(void) const{
+  return nFaceClosure;
+}
+
+inline unsigned int EvaluatedBasis::getNPoint(void) const{
+  return nPoint;
+}
+
+#endif
diff --git a/FunctionSpace/EvaluatedBasisScalar.cpp b/FunctionSpace/EvaluatedBasisScalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..036e819e079198c31068b8a995a65b906c5ddb79
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasisScalar.cpp
@@ -0,0 +1,124 @@
+#include "EvaluatedBasisScalar.h"
+
+using namespace std;
+
+EvaluatedBasisScalar::
+EvaluatedBasisScalar(const BasisScalar& basis, 
+		     const fullMatrix<double>& point){
+  // Data //
+  scalar = true;
+  
+  nVertex = basis.getNVertexBased();
+  nEdge   = basis.getNEdgeBased();
+  nFace   = basis.getNFaceBased();
+  nCell   = basis.getNCellBased();
+
+  nEdgeClosure = basis.getNEdgeClosure();
+  nFaceClosure = basis.getNFaceClosure();
+
+  nPoint = point.size1();
+
+
+  // Alloc //
+  // Node
+  node = new vector<vector<double>*>(nVertex);
+
+  for(unsigned int j = 0; j < nVertex; j++)
+    (*node)[j] = new vector<double>(nPoint);
+
+  // Edge
+  edge = new vector<vector<vector<double>*>*>(nEdgeClosure);
+
+  for(unsigned int i = 0; i < nEdgeClosure; i++){
+    (*edge)[i] = new vector<vector<double>*>(nEdge);
+    
+    for(unsigned int j = 0; j < nEdge; j++)
+      (*(*edge)[i])[j] = new vector<double>(nPoint);
+  }
+
+  // Face
+  face = new vector<vector<vector<double>*>*>(nFaceClosure);
+
+  for(unsigned int i = 0; i < nFaceClosure; i++){
+    (*face)[i] = new vector<vector<double>*>(nFace);
+
+    for(unsigned int j = 0; j < nFace; j++)
+      (*(*face)[i])[j] = new vector<double>(nPoint);
+  }  
+
+  // Cell
+  cell = new vector<vector<double>*>(nCell);
+
+  for(unsigned int j = 0; j < nCell; j++)
+    (*cell)[j] = new vector<double>(nPoint);
+
+
+  // Evaluate //
+  // Node
+  for(unsigned int j = 0; j < nVertex; j++)
+    for(unsigned int k = 0; k < nPoint; k++)
+      (*(*node)[j])[k] = 
+	basis.getNodeFunction(j).at(point(k, 0),
+				    point(k, 1),
+				    point(k, 2));
+
+  // Edge
+  for(unsigned int i = 0; i < nEdgeClosure; i++)
+    for(unsigned int j = 0; j < nEdge; j++)
+      for(unsigned int k = 0; k < nPoint; k++)
+	(*(*(*edge)[i])[j])[k] = 
+	  basis.getEdgeFunction(i, j).at(point(k, 0),
+					 point(k, 1),
+					 point(k, 2));
+  
+  // Face
+  for(unsigned int i = 0; i < nFaceClosure; i++)
+    for(unsigned int j = 0; j < nFace; j++)
+      for(unsigned int k = 0; k < nPoint; k++)
+	(*(*(*face)[i])[j])[k] = 
+	  basis.getFaceFunction(i, j).at(point(k, 0),
+					 point(k, 1),
+					 point(k, 2));
+
+  // Cell
+  for(unsigned int j = 0; j < nCell; j++)
+    for(unsigned int k = 0; k < nPoint; k++)
+      (*(*cell)[j])[k] = 
+	basis.getCellFunction(j).at(point(k, 0),
+				    point(k, 1),
+				    point(k, 2));
+}
+
+EvaluatedBasisScalar::~EvaluatedBasisScalar(void){
+  // Node //
+  for(unsigned int j = 0; j < nVertex; j++)
+    delete (*node)[j];
+
+  delete node;
+
+  // Edge //
+  for(unsigned int i = 0; i < nEdgeClosure; i++){
+    for(unsigned int j = 0; j < nEdge; j++)
+      delete (*(*edge)[i])[j];
+
+    delete (*edge)[i];
+  }
+
+  delete edge;
+
+  // Face //
+  for(unsigned int i = 0; i < nFaceClosure; i++){
+    for(unsigned int j = 0; j < nFace; j++)
+      delete (*(*face)[i])[j];
+
+    delete (*face)[i];
+  }  
+
+  delete face;
+
+  // Cell //
+  for(unsigned int j = 0; j < nCell; j++)
+    delete (*cell)[j];
+
+  delete cell;
+}
diff --git a/FunctionSpace/EvaluatedBasisScalar.h b/FunctionSpace/EvaluatedBasisScalar.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ecdeb0f4feece35c6c99d4c96afa258040f84eb
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasisScalar.h
@@ -0,0 +1,84 @@
+#ifndef _EVALUATEDBASISSCALAR_H_
+#define _EVALUATEDBASISSCALAR_H_
+
+#include <vector>
+#include "fullMatrix.h"
+#include "BasisScalar.h"
+#include "EvaluatedBasis.h"
+
+class EvaluatedBasisScalar: public EvaluatedBasis{
+ protected:
+  std::vector            <std::vector<double>*>*   node;
+  std::vector<std::vector<std::vector<double>*>*>* edge;
+  std::vector<std::vector<std::vector<double>*>*>* face;
+  std::vector            <std::vector<double>*>*   cell;
+
+ public:
+  //! @param basis the Basis to Evaluate
+  //! @param point the Evaluation Points
+  //!
+  //! Instanciates the requested EvaluatedBasisScalar
+  //!
+  EvaluatedBasisScalar(const BasisScalar& basis, 
+		       const fullMatrix<double>& point);
+
+  //! Deletes this EvaluatedBasisScalar
+  //!
+  virtual ~EvaluatedBasisScalar(void);
+
+  //! @param i A natural number
+  //! @return Returns the evaluation of the @c i%th @em Vertex Based 
+  //! Basis Function (at every evaluation Points)
+  const std::vector<double>&
+    getNodeEvaluation(unsigned int i) const;
+  
+  //! @param i A natural number
+  //! @param closure A natural number
+  //! @return Returns the evaluation of the @c i%th @em Edge Based 
+  //! Basis Function (at every evaluation Points), with the @c closure%th Closure
+  const std::vector<double>&
+    getEdgeEvaluation(unsigned int closure, unsigned int i) const;
+  
+  //! @param i A natural number
+  //! @param closure A natural number
+  //! @return Returns the evaluation of the @c i%th @em Face Based 
+  //! Basis Function (at every evaluation Points), with the @c closure%th Closure
+  const std::vector<double>&
+    getFaceEvaluation(unsigned int closure, unsigned int i) const;
+ 
+  //! @param i A natural number
+  //! @return Returns the evaluation of the @c i%th @em Cell Based 
+  //! Basis Function (at every evaluation Points)
+  const std::vector<double>&
+    getCellEvaluation(unsigned int i) const;
+};
+
+//////////////////////
+// Inline Function //
+//////////////////////
+
+inline
+const std::vector<double>& 
+EvaluatedBasisScalar::getNodeEvaluation(unsigned int i) const{
+  return *(*node)[i];
+}
+
+inline  
+const std::vector<double>& 
+EvaluatedBasisScalar::getEdgeEvaluation(unsigned int closure, unsigned int i) const{
+  return *(*(*edge)[closure])[i];
+}
+
+inline
+const std::vector<double>& 
+EvaluatedBasisScalar::getFaceEvaluation(unsigned int closure, unsigned int i) const{
+  return *(*(*face)[closure])[i];
+}
+
+inline
+const std::vector<double>&
+EvaluatedBasisScalar::getCellEvaluation(unsigned int i) const{
+  return *(*cell)[i];
+}
+
+#endif
diff --git a/FunctionSpace/EvaluatedBasisVector.cpp b/FunctionSpace/EvaluatedBasisVector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac437d38aab985dadd62a2c70c571ea280cb8450
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasisVector.cpp
@@ -0,0 +1,129 @@
+#include "Polynomial.h"
+#include "EvaluatedBasisVector.h"
+
+using namespace std;
+
+EvaluatedBasisVector::
+EvaluatedBasisVector(const BasisVector& basis, 
+		     const fullMatrix<double>& point){
+  // Data //
+  scalar = false;
+  
+  nVertex = basis.getNVertexBased();
+  nEdge   = basis.getNEdgeBased();
+  nFace   = basis.getNFaceBased();
+  nCell   = basis.getNCellBased();
+
+  nEdgeClosure = basis.getNEdgeClosure();
+  nFaceClosure = basis.getNFaceClosure();
+
+  nPoint = point.size1();
+
+
+  // Alloc //
+  // Node
+  node = new vector<vector<fullVector<double> >*>(nVertex);
+
+  for(unsigned int j = 0; j < nVertex; j++)
+    (*node)[j] = new vector<fullVector<double> >(nPoint);
+
+  // Edge
+  edge = new vector<vector<vector<fullVector<double> >*>*>(nEdgeClosure);
+
+  for(unsigned int i = 0; i < nEdgeClosure; i++){
+    (*edge)[i] = new vector<vector<fullVector<double> >*>(nEdge);
+    
+    for(unsigned int j = 0; j < nEdge; j++)
+      (*(*edge)[i])[j] = new vector<fullVector<double> >(nPoint);
+  }
+
+  // Face
+  face = new vector<vector<vector<fullVector<double> >*>*>(nFaceClosure);
+
+  for(unsigned int i = 0; i < nFaceClosure; i++){
+    (*face)[i] = new vector<vector<fullVector<double> >*>(nFace);
+
+    for(unsigned int j = 0; j < nFace; j++)
+      (*(*face)[i])[j] = new vector<fullVector<double> >(nPoint);
+  }  
+
+  // Cell
+  cell = new vector<vector<fullVector<double> >*>(nCell);
+
+  for(unsigned int j = 0; j < nCell; j++)
+    (*cell)[j] = new vector<fullVector<double> >(nPoint);
+
+
+  // Evaluate //
+  // Node
+  for(unsigned int j = 0; j < nVertex; j++)
+    for(unsigned int k = 0; k < nPoint; k++)
+      (*(*node)[j])[k] = 
+	Polynomial::at(basis.getNodeFunction(j),
+		       point(k, 0),
+		       point(k, 1),
+		       point(k, 2));
+
+  // Edge
+  for(unsigned int i = 0; i < nEdgeClosure; i++)
+    for(unsigned int j = 0; j < nEdge; j++)
+      for(unsigned int k = 0; k < nPoint; k++)
+	(*(*(*edge)[i])[j])[k] = 
+	  Polynomial::at(basis.getEdgeFunction(i, j),
+			 point(k, 0),
+			 point(k, 1),
+			 point(k, 2));
+  
+  // Face
+  for(unsigned int i = 0; i < nFaceClosure; i++)
+    for(unsigned int j = 0; j < nFace; j++)
+      for(unsigned int k = 0; k < nPoint; k++)
+	(*(*(*face)[i])[j])[k] = 
+	  Polynomial::at(basis.getFaceFunction(i, j),
+			 point(k, 0),
+			 point(k, 1),
+			 point(k, 2));
+
+  // Cell
+  for(unsigned int j = 0; j < nCell; j++)
+    for(unsigned int k = 0; k < nPoint; k++)
+      (*(*cell)[j])[k] = 
+	Polynomial::at(basis.getCellFunction(j),
+		       point(k, 0),
+		       point(k, 1),
+		       point(k, 2));
+}
+
+EvaluatedBasisVector::~EvaluatedBasisVector(void){
+  // Node //
+  for(unsigned int j = 0; j < nVertex; j++)
+    delete (*node)[j];
+
+  delete node;
+
+  // Edge //
+  for(unsigned int i = 0; i < nEdgeClosure; i++){
+    for(unsigned int j = 0; j < nEdge; j++)
+      delete (*(*edge)[i])[j];
+
+    delete (*edge)[i];
+  }
+
+  delete edge;
+
+  // Face //
+  for(unsigned int i = 0; i < nFaceClosure; i++){
+    for(unsigned int j = 0; j < nFace; j++)
+      delete (*(*face)[i])[j];
+
+    delete (*face)[i];
+  }  
+
+  delete face;
+
+  // Cell //
+  for(unsigned int j = 0; j < nCell; j++)
+    delete (*cell)[j];
+
+  delete cell;
+}
diff --git a/FunctionSpace/EvaluatedBasisVector.h b/FunctionSpace/EvaluatedBasisVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..b9b686df74b6a6a8d8fc40965b5daa394e5c256c
--- /dev/null
+++ b/FunctionSpace/EvaluatedBasisVector.h
@@ -0,0 +1,84 @@
+#ifndef _EVALUATEDBASISVECTOR_H_
+#define _EVALUATEDBASISVECTOR_H_
+
+#include <vector>
+#include "fullMatrix.h"
+#include "BasisVector.h"
+#include "EvaluatedBasis.h"
+
+class EvaluatedBasisVector: public EvaluatedBasis{
+ protected:
+  std::vector            <std::vector<fullVector<double> >*>*   node;
+  std::vector<std::vector<std::vector<fullVector<double> >*>*>* edge;
+  std::vector<std::vector<std::vector<fullVector<double> >*>*>* face;
+  std::vector            <std::vector<fullVector<double> >*>*   cell;
+
+ public:
+  //! @param basis the Basis to Evaluate
+  //! @param point the Evaluation Points
+  //!
+  //! Instanciates the requested EvaluatedBasisVector
+  //!
+  EvaluatedBasisVector(const BasisVector& basis, 
+		       const fullMatrix<double>& point);
+
+  //! Deletes this EvaluatedBasisVector
+  //!
+  virtual ~EvaluatedBasisVector(void);
+
+  //! @param i A natural number
+  //! @return Returns the evaluation of the @c i%th @em Vertex Based 
+  //! Basis Function (at every evaluation Points)
+  const std::vector<fullVector<double> >&
+    getNodeEvaluation(unsigned int i) const;
+  
+  //! @param i A natural number
+  //! @param closure A natural number
+  //! @return Returns the evaluation of the @c i%th @em Edge Based 
+  //! Basis Function (at every evaluation Points), with the @c closure%th Closure
+  const std::vector<fullVector<double> >&
+    getEdgeEvaluation(unsigned int closure, unsigned int i) const;
+  
+  //! @param i A natural number
+  //! @param closure A natural number
+  //! @return Returns the evaluation of the @c i%th @em Face Based 
+  //! Basis Function (at every evaluation Points), with the @c closure%th Closure
+  const std::vector<fullVector<double> >&
+    getFaceEvaluation(unsigned int closure, unsigned int i) const;
+ 
+  //! @param i A natural number
+  //! @return Returns the evaluation of the @c i%th @em Cell Based 
+  //! Basis Function (at every evaluation Points)
+  const std::vector<fullVector<double> >&
+    getCellEvaluation(unsigned int i) const;
+};
+
+//////////////////////
+// Inline Function //
+//////////////////////
+
+inline
+const std::vector<fullVector<double> >& 
+EvaluatedBasisVector::getNodeEvaluation(unsigned int i) const{
+  return *(*node)[i];
+}
+
+inline  
+const std::vector<fullVector<double> >& 
+EvaluatedBasisVector::getEdgeEvaluation(unsigned int closure, unsigned int i) const{
+  return *(*(*edge)[closure])[i];
+}
+
+inline
+const std::vector<fullVector<double> >& 
+EvaluatedBasisVector::getFaceEvaluation(unsigned int closure, unsigned int i) const{
+  return *(*(*face)[closure])[i];
+}
+
+inline
+const std::vector<fullVector<double> >&
+EvaluatedBasisVector::getCellEvaluation(unsigned int i) const{
+  return *(*cell)[i];
+}
+
+#endif
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index ffbce1bc6cf39fc7d6dc20b3fd79110b87c7bd64..a0c6d9848f06d85061a54ff09b04ede32c4774fb 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -395,7 +395,7 @@ FunctionSpace::locBasis(const MElement& element,
   // Edge Based
   // Number of basis function *per* edge
   //  --> should always be an integer !
-  const unsigned int nEdge     = edgeClosure.size();
+  const unsigned int nEdge = edgeClosure.size();
 
   if(nEdge){
     const unsigned int nFPerEdge = nFEdge / nEdge;
@@ -415,7 +415,7 @@ FunctionSpace::locBasis(const MElement& element,
   // Face Based
   // Number of basis function *per* face
   //  --> should always be an integer !
-  const unsigned int nFace     = faceClosure.size();
+  const unsigned int nFace = faceClosure.size();
 
   if(nFace){
     const unsigned int nFPerFace = nFFace / nFace;
@@ -443,6 +443,155 @@ FunctionSpace::locBasis(const MElement& element,
   return fun;  
 }
 
+const vector<const vector<double>*>
+FunctionSpace::locEvalBasis(const MElement& element, 
+			    const EvaluatedBasisScalar& evalBasis){
+  // Get Basis //
+  const unsigned int nFNode = evalBasis.getNVertexBased();
+  const unsigned int nFEdge = evalBasis.getNEdgeBased();
+  const unsigned int nFFace = evalBasis.getNFaceBased();
+  const unsigned int nFCell = evalBasis.getNCellBased();
+
+  // Get Closure //
+  vector<int> edgeClosure = getEdgeClosure(element);
+  vector<int> faceClosure = getFaceClosure(element); 
+
+  // Get Functions //
+  vector<const vector<double>*> fun(nFNode + nFEdge + nFFace + nFCell);
+  unsigned int i = 0;
+  
+  // Vertex Based
+  for(unsigned int j = 0; j < nFNode; j++){
+    fun[i] = &evalBasis.getNodeEvaluation(j);
+    i++;
+  }
+
+  // Edge Based
+  // Number of basis function *per* edge
+  //  --> should always be an integer !
+  const unsigned int nEdge = edgeClosure.size();
+
+  if(nEdge){
+    const unsigned int nFPerEdge = nFEdge / nEdge;
+    unsigned int fEdge = 0;
+    
+    for(unsigned int j = 0; j < nFPerEdge; j++){
+      for(unsigned int k = 0; k < nEdge; k++){
+	fun[i] = 
+	  &evalBasis.getEdgeEvaluation(edgeClosure[k], fEdge);
+	
+	fEdge++;
+	i++;
+      }
+    }
+  }
+
+  // Face Based
+  // Number of basis function *per* face
+  //  --> should always be an integer !
+  const unsigned int nFace = faceClosure.size();
+
+  if(nFace){
+    const unsigned int nFPerFace = nFFace / nFace;
+    unsigned int fFace = 0;
+    
+    for(unsigned int j = 0; j < nFPerFace; j++){
+      for(unsigned int k = 0; k < nFace; k++){
+	fun[i] = 
+	  &evalBasis.getFaceEvaluation(faceClosure[k], fFace);
+	
+	fFace++;
+	i++;
+      }
+    }
+  }
+
+  // Cell Based
+  for(unsigned int j = 0; j < nFCell; j++){
+    fun[i] = &evalBasis.getCellEvaluation(j);
+    i++;
+  }
+
+
+  // Return //
+  return fun;  
+}
+
+const vector<const vector<fullVector<double> >*>
+FunctionSpace::locEvalBasis(const MElement& element, 
+			    const EvaluatedBasisVector& evalBasis){
+  // Get Basis //
+  const unsigned int nFNode = evalBasis.getNVertexBased();
+  const unsigned int nFEdge = evalBasis.getNEdgeBased();
+  const unsigned int nFFace = evalBasis.getNFaceBased();
+  const unsigned int nFCell = evalBasis.getNCellBased();
+
+  // Get Closure //
+  vector<int> edgeClosure = getEdgeClosure(element);
+  vector<int> faceClosure = getFaceClosure(element); 
+
+  // Get Functions //
+  vector<const vector<fullVector<double> >*> 
+    fun(nFNode + nFEdge + nFFace + nFCell);
+  
+  unsigned int i = 0;
+  
+  // Vertex Based
+  for(unsigned int j = 0; j < nFNode; j++){
+    fun[i] = &evalBasis.getNodeEvaluation(j);
+    i++;
+  }
+
+  // Edge Based
+  // Number of basis function *per* edge
+  //  --> should always be an integer !
+  const unsigned int nEdge = edgeClosure.size();
+
+  if(nEdge){
+    const unsigned int nFPerEdge = nFEdge / nEdge;
+    unsigned int fEdge = 0;
+    
+    for(unsigned int j = 0; j < nFPerEdge; j++){
+      for(unsigned int k = 0; k < nEdge; k++){
+	fun[i] = 
+	  &evalBasis.getEdgeEvaluation(edgeClosure[k], fEdge);
+	
+	fEdge++;
+	i++;
+      }
+    }
+  }
+
+  // Face Based
+  // Number of basis function *per* face
+  //  --> should always be an integer !
+  const unsigned int nFace = faceClosure.size();
+
+  if(nFace){
+    const unsigned int nFPerFace = nFFace / nFace;
+    unsigned int fFace = 0;
+    
+    for(unsigned int j = 0; j < nFPerFace; j++){
+      for(unsigned int k = 0; k < nFace; k++){
+	fun[i] = 
+	  &evalBasis.getFaceEvaluation(faceClosure[k], fFace);
+	
+	fFace++;
+	i++;
+      }
+    }
+  }
+
+  // Cell Based
+  for(unsigned int j = 0; j < nFCell; j++){
+    fun[i] = &evalBasis.getCellEvaluation(j);
+    i++;
+  }
+
+
+  // Return //
+  return fun;  
+}
 
 string FunctionSpace::toString(void) const{
   return basis->toString();
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 74ce8e90f663e6a8ae65edcc24bf978d3e8dc6a0..82d7dd64ad41af6d6a28d15e01da329896883f7a 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -9,6 +9,9 @@
 #include "BasisScalar.h"
 #include "BasisVector.h"
 
+#include "EvaluatedBasisScalar.h"
+#include "EvaluatedBasisVector.h"
+
 #include "Comparators.h"
 #include "Dof.h"
 #include "GroupOfDof.h"
@@ -117,6 +120,16 @@ class FunctionSpace{
     const std::vector<const std::vector<Polynomial>*>
     locBasis(const MElement& element, 
 	     const BasisVector& basis);
+
+  static 
+    const std::vector<const std::vector<double>*>
+    locEvalBasis(const MElement& element, 
+		 const EvaluatedBasisScalar& evalBasis);
+
+  static 
+    const std::vector<const std::vector<fullVector<double> >*>
+    locEvalBasis(const MElement& element,
+		 const EvaluatedBasisVector& evalBasis);
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index e3f3447793cc9f2c7d28b4ce878533f0d0c15704..8936e00335e42aca34caf1e841fddf1ce5e55648 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -5,9 +5,52 @@ using namespace std;
 FunctionSpaceScalar::FunctionSpaceScalar(void){
   hasGrad   = false;
   gradBasis = NULL;
+  
+  locPreEvaluated  = false;
+  gradPreEvaluated = false;  
 }
 
 FunctionSpaceScalar::~FunctionSpaceScalar(void){
   if(hasGrad)
     delete gradBasis;
+  
+  if(locPreEvaluated)
+    delete evalLoc;
+
+  if(gradPreEvaluated)
+    delete evalGrad;
+}
+
+void FunctionSpaceScalar::
+preEvaluateLocalFunctions(fullMatrix<double>& points){
+  // Delete Old Struct (if any) //
+  if(locPreEvaluated)
+    delete evalLoc;
+
+  // New Struct //
+  evalLoc = new EvaluatedBasisScalar(*basisScalar, points);
+
+  // PreEvaluated //
+  locPreEvaluated = true;
+}
+
+void FunctionSpaceScalar::
+preEvaluateGradLocalFunctions(fullMatrix<double>& points){
+  // Got Grad Basis ? //
+  // --> mutable data 
+  //  --> Just a 'cache memory' 
+  if(!hasGrad){
+    gradBasis = new GradBasis(*basisScalar);
+    hasGrad   = true;
+  }
+
+  // Delete Old Struct (if any) //
+  if(gradPreEvaluated)
+    delete evalGrad;
+
+  // New Struct //
+  evalGrad = new EvaluatedBasisVector(*gradBasis, points);
+
+  // PreEvaluated //
+  gradPreEvaluated = true;
 }
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index b4323f50832846ef5be06b7805ec27b1f2e32d10..de7c3ababfbaa78dbacafaefcdee6afcd9028a7f 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -2,8 +2,11 @@
 #define _FUNCTIONSPACESCALAR_H_
 
 #include "fullMatrix.h"
+#include "Exception.h"
 #include "BasisScalar.h"
 #include "GradBasis.h"
+#include "EvaluatedBasisScalar.h"
+#include "EvaluatedBasisVector.h"
 #include "FunctionSpace.h"
 
 /**
@@ -30,6 +33,12 @@ class FunctionSpaceScalar : public FunctionSpace{
   mutable bool       hasGrad;
   mutable GradBasis* gradBasis;
 
+  bool locPreEvaluated;
+  bool gradPreEvaluated;
+  
+  EvaluatedBasisScalar* evalLoc;
+  EvaluatedBasisVector* evalGrad;
+
  public:
   virtual ~FunctionSpaceScalar(void);
 
@@ -49,6 +58,15 @@ class FunctionSpaceScalar : public FunctionSpace{
   const std::vector<const std::vector<Polynomial>*>
     getGradLocalFunctions(const MElement& element) const;
 
+  void preEvaluateLocalFunctions(fullMatrix<double>& points);
+  void preEvaluateGradLocalFunctions(fullMatrix<double>& points);
+
+  const std::vector<const std::vector<double>*>
+    getEvaluatedLocalFunctions(const MElement& element) const;
+
+  const std::vector<const std::vector<fullVector<double> >*>
+    getEvaluatedGradLocalFunctions(const MElement& element) const;
+  
  protected:
   FunctionSpaceScalar(void);
 };
@@ -121,7 +139,7 @@ FunctionSpaceScalar::getLocalFunctions(const MElement& element) const{
 
 inline const std::vector<const std::vector<Polynomial>*>
 FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
-
+  
   // Got Grad Basis ? //
   // --> mutable data 
   //  --> Just a 'cache memory' 
@@ -129,8 +147,24 @@ FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
     gradBasis = new GradBasis(*basisScalar);
     hasGrad   = true;
   }
-
+  
   return locBasis(element, *gradBasis);
 }
 
+inline const std::vector<const std::vector<double>*>
+FunctionSpaceScalar::getEvaluatedLocalFunctions(const MElement& element) const{
+  if(!locPreEvaluated)
+    throw Exception("Local Basis Functions not PreEvaluated");
+
+  return locEvalBasis(element, *evalLoc);
+}
+
+inline const std::vector<const std::vector<fullVector<double> >*>
+FunctionSpaceScalar::getEvaluatedGradLocalFunctions(const MElement& element) const{
+  if(!gradPreEvaluated)
+    throw Exception("Gradients of Local Basis Functions not PreEvaluated");
+
+  return locEvalBasis(element, *evalGrad);
+}
+
 #endif