diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 868209ecad7f281037bdabf15bd84c20f3413c1f..fc7ce91029320816d46aa851b476b7c421841ffe 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -295,6 +295,143 @@ vector<int> FunctionSpace::getFaceClosure(const MElement& element){
   return faceClosure;
 }
 
+const vector<const Polynomial*>
+FunctionSpace::locBasis(const MElement& element, 
+			const BasisScalar& basis){
+  // Get Basis //
+  const unsigned int nFNode = basis.getNVertexBased();
+  const unsigned int nFEdge = basis.getNEdgeBased();
+  const unsigned int nFFace = basis.getNFaceBased();
+  const unsigned int nFCell = basis.getNCellBased();
+
+  // Get Closure //
+  vector<int> edgeClosure = getEdgeClosure(element);
+  vector<int> faceClosure = getFaceClosure(element); 
+
+  // Get Functions //
+  vector<const Polynomial*> fun(basis.getSize());
+  unsigned int i = 0;
+  
+  // Vertex Based
+  for(unsigned int j = 0; j < nFNode; j++){
+    fun[i] = &basis.getNodeFunction(j);
+    i++;
+  }
+
+  // Edge Based
+  // Number of basis function *per* edge
+  //  --> should always be an integer !
+  const unsigned int nEdge     = edgeClosure.size();
+  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] = 
+	&basis.getEdgeFunction(edgeClosure[k], fEdge);
+      
+      fEdge++;
+      i++;
+    }
+  }
+
+  // Face Based
+  // Number of basis function *per* face
+  //  --> should always be an integer !
+  const unsigned int nFace     = faceClosure.size();
+  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] = 
+	&basis.getFaceFunction(faceClosure[k], fFace);
+      
+      fFace++;
+      i++;
+    }
+  }
+
+  // Cell Based
+  for(unsigned int j = 0; j < nFCell; j++){
+    fun[i] = &basis.getCellFunction(j);
+    i++;
+  }
+
+
+  // Return //
+  return fun;
+}
+
+const vector<const vector<Polynomial>*>
+FunctionSpace::locBasis(const MElement& element, 
+			const BasisVector& basis){
+  // Get Basis //
+  const unsigned int nFNode = basis.getNVertexBased();
+  const unsigned int nFEdge = basis.getNEdgeBased();
+  const unsigned int nFFace = basis.getNFaceBased();
+  const unsigned int nFCell = basis.getNCellBased();
+
+  // Get Closure //
+  vector<int> edgeClosure = getEdgeClosure(element);
+  vector<int> faceClosure = getFaceClosure(element); 
+
+  // Get Functions //
+  vector<const vector<Polynomial>*> fun(basis.getSize());
+  unsigned int i = 0;
+  
+  // Vertex Based
+  for(unsigned int j = 0; j < nFNode; j++){
+    fun[i] = &basis.getNodeFunction(j);
+    i++;
+  }
+
+  // Edge Based
+  // Number of basis function *per* edge
+  //  --> should always be an integer !
+  const unsigned int nEdge     = edgeClosure.size();
+  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] = 
+	&basis.getEdgeFunction(edgeClosure[k], fEdge);
+      
+      fEdge++;
+      i++;
+    }
+  }
+
+  // Face Based
+  // Number of basis function *per* face
+  //  --> should always be an integer !
+  const unsigned int nFace     = faceClosure.size();
+  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] = 
+	&basis.getFaceFunction(faceClosure[k], fFace);
+      
+      fFace++;
+      i++;
+    }
+  }
+
+  // Cell Based
+  for(unsigned int j = 0; j < nFCell; j++){
+    fun[i] = &basis.getCellFunction(j);
+    i++;
+  }
+
+
+  // Return //
+  return fun;  
+}
+
+
 string FunctionSpace::toString(void) const{
   return basis->toString();
 }
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 18f0cf387354a07a783b9aea50f2c53d16564da0..76cba19839739302d6323a3256790ec9bea2165a 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -6,6 +6,9 @@
 #include <string>
 
 #include "Basis.h"
+#include "BasisScalar.h"
+#include "BasisVector.h"
+
 #include "Comparators.h"
 #include "Dof.h"
 #include "GroupOfDof.h"
@@ -99,6 +102,16 @@ class FunctionSpace{
   // Closure
   static std::vector<int> getEdgeClosure(const MElement& element);
   static std::vector<int> getFaceClosure(const MElement& element);
+
+  // Local Basis
+  static 
+    const std::vector<const Polynomial*>
+    locBasis(const MElement& element, 
+	     const BasisScalar& basis);
+  static 
+    const std::vector<const std::vector<Polynomial>*>
+    locBasis(const MElement& element, 
+	     const BasisVector& basis);
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index 42effb1d876ab00ad7a15613692c1a66cdedd76f..1907ca14144a20ced7677445ae4898c6f1ed7be4 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -10,6 +10,10 @@ FunctionSpaceEdge::FunctionSpaceEdge(const GroupOfElement& goe,
 				     int order){
   // Build 1Form Basis //
   build(goe, 1, order); 
+
+  // Init BasisVector //
+  basisVector = 
+    static_cast<const BasisVector*>(basis);
 }
     
 FunctionSpaceEdge::~FunctionSpaceEdge(void){
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 8480ee7543e27dae8b6d5e1fd61d8005b2d45ddc..e3f3447793cc9f2c7d28b4ce878533f0d0c15704 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -1,5 +1,4 @@
 #include "FunctionSpaceScalar.h"
-#include "Exception.h"
 
 using namespace std;
 
@@ -12,147 +11,3 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
   if(hasGrad)
     delete gradBasis;
 }
-
-const vector<const Polynomial*> FunctionSpaceScalar::
-getLocalFunctions(const MElement& element) const{
-  
-  // Get Basis //
-  const unsigned int nFNode = basisScalar->getNVertexBased();
-  const unsigned int nFEdge = basisScalar->getNEdgeBased();
-  const unsigned int nFFace = basisScalar->getNFaceBased();
-  const unsigned int nFCell = basisScalar->getNCellBased();
-
-  // Get Closure //
-  vector<int> edgeClosure = getEdgeClosure(element);
-  vector<int> faceClosure = getFaceClosure(element); 
-
-  // Get Functions //
-  vector<const Polynomial*> fun(basisScalar->getSize());
-  unsigned int i = 0;
-  
-  // Vertex Based
-  for(unsigned int j = 0; j < nFNode; j++){
-    fun[i] = &basisScalar->getNodeFunction(j);
-    i++;
-  }
-
-  // Edge Based
-  // Number of basis function *per* edge
-  //  --> should always be an integer !
-  const unsigned int nEdge     = edgeClosure.size();
-  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] = 
-	&basisScalar->getEdgeFunction(edgeClosure[k], fEdge);
-      
-      fEdge++;
-      i++;
-    }
-  }
-
-  // Face Based
-  // Number of basis function *per* face
-  //  --> should always be an integer !
-  const unsigned int nFace     = faceClosure.size();
-  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] = 
-	&basisScalar->getFaceFunction(faceClosure[k], fFace);
-      
-      fFace++;
-      i++;
-    }
-  }
-
-  // Cell Based
-  for(unsigned int j = 0; j < nFCell; j++){
-    fun[i] = &basisScalar->getCellFunction(j);
-    i++;
-  }
-
-
-  // Return //
-  return fun;
-}
-
-const vector<const vector<Polynomial>*> FunctionSpaceScalar::
-getGradLocalFunctions(const MElement& element) const{
-
-  // Got Grad Basis ? //
-  // --> mutable data 
-  //  --> Just a 'cache memory' 
-  if(!hasGrad){
-    gradBasis = new GradBasis(*basisScalar);
-    hasGrad   = true;
-  }
-
-  // Get Basis //
-  const unsigned int nFNode = gradBasis->getNVertexBased();
-  const unsigned int nFEdge = gradBasis->getNEdgeBased();
-  const unsigned int nFFace = gradBasis->getNFaceBased();
-  const unsigned int nFCell = gradBasis->getNCellBased();
-
-  // Get Closure //
-  vector<int> edgeClosure = getEdgeClosure(element);
-  vector<int> faceClosure = getFaceClosure(element); 
-
-  // Get Functions //
-  vector<const vector<Polynomial>*> fun(gradBasis->getSize());
-  unsigned int i = 0;
-  
-  // Vertex Based
-  for(unsigned int j = 0; j < nFNode; j++){
-    fun[i] = &gradBasis->getNodeFunction(j);
-    i++;
-  }
-
-  // Edge Based
-  // Number of basis function *per* edge
-  //  --> should always be an integer !
-  const unsigned int nEdge     = edgeClosure.size();
-  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] = 
-	&gradBasis->getEdgeFunction(edgeClosure[k], fEdge);
-      
-      fEdge++;
-      i++;
-    }
-  }
-
-  // Face Based
-  // Number of basis function *per* face
-  //  --> should always be an integer !
-  const unsigned int nFace     = faceClosure.size();
-  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] = 
-	&gradBasis->getFaceFunction(faceClosure[k], fFace);
-      
-      fFace++;
-      i++;
-    }
-  }
-
-  // Cell Based
-  for(unsigned int j = 0; j < nFCell; j++){
-    fun[i] = &gradBasis->getCellFunction(j);
-    i++;
-  }
-
-
-  // Return //
-  return fun;  
-}
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 8b0999e0b70746fa8924bde330f76829de4aab2a..18a86c9773eb28265dcd3e581b750f3d633d5062 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -41,8 +41,6 @@ class FunctionSpaceScalar : public FunctionSpace{
 
   const std::vector<const std::vector<Polynomial>*>
     getGradLocalFunctions(const MElement& element) const;
-  
-  const BasisScalar& getBasis(const MElement& element) const;
 
  protected:
   FunctionSpaceScalar(void);
@@ -77,23 +75,29 @@ class FunctionSpaceScalar : public FunctionSpace{
    @param element A MElement
    @return Returns the basis functions associated
    to the given element (with correct @em closure)
-   **
-
-   @fn FunctionSpaceScalar::getBasis
-   @param element A MElement of the support 
-   of this FunctionSpace
-   @return Returns the Basis (BasisScalar) associated
-   to the given MElement
  */
 
-
 //////////////////////
 // Inline Functions //
 //////////////////////
-/*
-inline const BasisScalar& FunctionSpaceScalar::
-getBasis(const MElement& element) const{
-  return *basisScalar;
+
+inline const std::vector<const Polynomial*> 
+FunctionSpaceScalar::getLocalFunctions(const MElement& element) const{
+  return locBasis(element, *basisScalar);
 }
-*/
+
+inline const std::vector<const std::vector<Polynomial>*>
+FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
+
+  // Got Grad Basis ? //
+  // --> mutable data 
+  //  --> Just a 'cache memory' 
+  if(!hasGrad){
+    gradBasis = new GradBasis(*basisScalar);
+    hasGrad   = true;
+  }
+
+  return locBasis(element, *gradBasis);
+}
+
 #endif
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index c32694405427fd031b5454602efa11b57a1e5051..fc9b5337be9b7d7519263c7ec7676d549a22be25 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -1,77 +1,19 @@
 #include "FunctionSpaceVector.h"
-#include "Exception.h"
 
 using namespace std;
 
-FunctionSpaceVector::~FunctionSpaceVector(void){
-}
-
-const vector<const vector<Polynomial>*> FunctionSpaceVector::
-getLocalFunctions(const MElement& element) const{
-
-  // Get Basis //
-  const BasisVector&  basis = getBasis(element);
-
-  const unsigned int nFNode = basis.getNVertexBased();
-  const unsigned int nFEdge = basis.getNEdgeBased();
-  const unsigned int nFFace = basis.getNFaceBased();
-  const unsigned int nFCell = basis.getNCellBased();
-
-  // Get Closure //
-  vector<int> edgeClosure = getEdgeClosure(element);
-  vector<int> faceClosure = getFaceClosure(element); 
-
-  // Get Functions //
-  vector<const vector<Polynomial>*> fun(basis.getSize());
-  unsigned int i = 0;
-  
-  // Vertex Based
-  for(unsigned int j = 0; j < nFNode; j++){
-    fun[i] = &basis.getNodeFunction(j);
-    i++;
-  }
+FunctionSpaceVector::FunctionSpaceVector(void){
+  hasCurl   = false;
+  curlBasis = NULL;
 
-  // Edge Based
-  // Number of basis function *per* edge
-  //  --> should always be an integer !
-  const unsigned int nEdge     = edgeClosure.size();
-  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] = 
-	&basis.getEdgeFunction(edgeClosure[k], fEdge);
-      
-      fEdge++;
-      i++;
-    }
-  }
-
-  // Face Based
-  // Number of basis function *per* face
-  //  --> should always be an integer !
-  const unsigned int nFace     = faceClosure.size();
-  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] = 
-	&basis.getFaceFunction(faceClosure[k], fFace);
-      
-      fFace++;
-      i++;
-    }
-  }
-
-  // Cell Based
-  for(unsigned int j = 0; j < nFCell; j++){
-    fun[i] = &basis.getCellFunction(j);
-    i++;
-  }
+  hasDiv   = false;
+  divBasis = NULL;
+}
 
+FunctionSpaceVector::~FunctionSpaceVector(void){
+  if(hasCurl)
+    delete curlBasis;
 
-  // Return //
-  return fun;
+  if(hasDiv)
+    delete divBasis;
 }
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 30d2ddb9140a9b6ce11cebdbe7b80e1bbfdd8dce..3b4ff01119df8033b64cbc57a0a8ef1108b28d55 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -3,6 +3,8 @@
 
 #include "fullMatrix.h"
 #include "BasisVector.h"
+#include "CurlBasis.h"
+#include "DivBasis.h"
 #include "FunctionSpace.h"
 
 /**
@@ -21,6 +23,15 @@
 
 
 class FunctionSpaceVector : public FunctionSpace{
+ protected:
+  const BasisVector* basisVector;
+
+  mutable bool       hasCurl;
+  mutable CurlBasis* curlBasis;
+
+  mutable bool       hasDiv;
+  mutable DivBasis*  divBasis;
+
  public:
   virtual ~FunctionSpaceVector(void);
 
@@ -32,7 +43,14 @@ class FunctionSpaceVector : public FunctionSpace{
   const std::vector<const std::vector<Polynomial>*>
     getLocalFunctions(const MElement& element) const;
 
-  const BasisVector& getBasis(const MElement& element) const;
+  const std::vector<const std::vector<Polynomial>*>
+    getCurlLocalFunctions(const MElement& element) const;
+
+  const std::vector<const Polynomial*> 
+    getDivLocalFunctions(const MElement& element) const;
+
+ protected:
+  FunctionSpaceVector(void);
 };
 
 
@@ -64,23 +82,43 @@ class FunctionSpaceVector : public FunctionSpace{
    @param element A MElement
    @return Returns the basis functions associated
    to the given element (with correct @em closure)
-   **
-
-   @fn FunctionSpaceVector::getBasis
-   @param element A MElement of the support 
-   of this FunctionSpace
-   @return Returns the Basis (BasisVector) associated
-   to the given MElement
- */
-
+*/
 
 //////////////////////
 // Inline Functions //
 //////////////////////
 
-inline const BasisVector& FunctionSpaceVector::
-getBasis(const MElement& element) const{
-  return static_cast<const BasisVector&>(*basis);
+inline const std::vector<const std::vector<Polynomial>*> 
+FunctionSpaceVector::getLocalFunctions(const MElement& element) const{
+  return locBasis(element, *basisVector);
+}
+
+inline const std::vector<const std::vector<Polynomial>*>
+FunctionSpaceVector::getCurlLocalFunctions(const MElement& element) const{
+
+  // Got Curl Basis ? //
+  // --> mutable data 
+  //  --> Just a 'cache memory' 
+  if(!hasCurl){
+    curlBasis = new CurlBasis(*basisVector);
+    hasCurl   = true;
+  }
+
+  return locBasis(element, *curlBasis);
+}
+
+inline const std::vector<const Polynomial*> 
+FunctionSpaceVector::getDivLocalFunctions(const MElement& element) const{
+
+  // Got Div Basis ? //
+  // --> mutable data 
+  //  --> Just a 'cache memory' 
+  if(!hasDiv){
+    divBasis = new DivBasis(*basisVector);
+    hasDiv   = true;
+  }
+
+  return locBasis(element, *divBasis);
 }
 
 #endif