diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index 1907ca14144a20ced7677445ae4898c6f1ed7be4..f7dc3b54ac08488cb48918b8e963870ba15dd669 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -61,3 +61,40 @@ interpolate(const MElement& element,
   // Return Interpolated Value //
   return val;
 }
+
+fullVector<double> FunctionSpaceEdge::
+interpolateInRefSpace(const MElement& element, 
+		      const std::vector<double>& coef,
+		      const fullVector<double>& uvw) const{
+ 
+  // Const Cast For MElement //
+  MElement& eelement = 
+    const_cast<MElement&>(element);
+  
+  // Get Basis Functions //
+  const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
+  const unsigned int nFun                      = fun.size();
+
+  // Get Jacobian //
+  fullMatrix<double>  invJac(3, 3);        
+  eelement.getJacobian(uvw(0), uvw(1), uvw(2), invJac);
+  invJac.invertInPlace();
+ 
+  // Interpolate (in Reference Place) //
+  fullVector<double> val(3); 
+  val(0) = 0; 
+  val(1) = 0; 
+  val(2) = 0;
+
+  for(unsigned int i = 0; i < nFun; i++){
+    fullVector<double> vi = 
+      Mapper::grad(Polynomial::at(*fun[i], uvw(0), uvw(1), uvw(2)),
+		   invJac);
+    
+    vi.scale(coef[i]);
+    val.axpy(vi, 1);
+  }
+
+  // Return Interpolated Value //
+  return val;
+}
diff --git a/FunctionSpace/FunctionSpaceEdge.h b/FunctionSpace/FunctionSpaceEdge.h
index bc03737d2c7987355ab1953011f57db575788829..0fbbaf9b11e7efec5a8a10cd1abe7b8b18ffd50c 100644
--- a/FunctionSpace/FunctionSpaceEdge.h
+++ b/FunctionSpace/FunctionSpaceEdge.h
@@ -21,6 +21,11 @@ class FunctionSpaceEdge : public FunctionSpaceVector{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const;
+
+  virtual fullVector<double> 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const;
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index 10928c769c632edda0f450eaa29eaa0fa5e1769a..0adb1198b6f35ce4a2ae5de059d9d4b1cda83dcc 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/FunctionSpace/FunctionSpaceNode.cpp
@@ -48,3 +48,23 @@ interpolate(const MElement& element,
   // Return Interpolated Value //
   return val;
 }
+
+double FunctionSpaceNode::
+interpolateInRefSpace(const MElement& element, 
+		      const std::vector<double>& coef,
+		      const fullVector<double>& uvw) const{
+
+  // Get Basis Functions //
+  const vector<const Polynomial*> fun = getLocalFunctions(element);
+  const unsigned int nFun             = fun.size();
+  
+  // Interpolate (in Reference Place) //
+  double val = 0; 
+
+  for(unsigned int i = 0; i < nFun; i++)
+    val += 
+      fun[i]->at(uvw(0), uvw(1), uvw(2)) * coef[i];
+
+  // Return Interpolated Value //
+  return val;
+}
diff --git a/FunctionSpace/FunctionSpaceNode.h b/FunctionSpace/FunctionSpaceNode.h
index 59655bcd4839e8a287b1461b2b725d2db5a649c1..d682c7da9e63da2c7f9f443c197f1d29a3ed49e0 100644
--- a/FunctionSpace/FunctionSpaceNode.h
+++ b/FunctionSpace/FunctionSpaceNode.h
@@ -21,6 +21,11 @@ class FunctionSpaceNode : public FunctionSpaceScalar{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const;
+
+  virtual double 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const;
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 18a86c9773eb28265dcd3e581b750f3d633d5062..c51c27458b2ea5f45285978f7418dd425624fed8 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -35,6 +35,11 @@ class FunctionSpaceScalar : public FunctionSpace{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const = 0;
+
+  virtual double 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const = 0;
   
   const std::vector<const Polynomial*> 
     getLocalFunctions(const MElement& element) const;
@@ -57,7 +62,26 @@ class FunctionSpaceScalar : public FunctionSpace{
    @param coef The coefficients of the interpolation
    @param xyz The coordinate 
    (of point @em inside the given @c element)
-   of the interpolation
+   of the interpolation in the @em Physical Space
+
+   @return Returns the (scalar) interpolated value
+
+   @warning
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   
+   @todo
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   ---> check
+   **
+
+   @fn FunctionSpaceScalar::interpolateInRefSpace
+   @param element The MElement to interpolate on
+   @param coef The coefficients of the interpolation
+   @param uvw The coordinate 
+   (of point @em inside the given @c element)
+   of the interpolation in the @em Reference Space
 
    @return Returns the (scalar) interpolated value
 
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 3b4ff01119df8033b64cbc57a0a8ef1108b28d55..7227f306b283d9c0b7b2cbc4b718191153f14379 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -40,6 +40,11 @@ class FunctionSpaceVector : public FunctionSpace{
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const = 0;
 
+  virtual fullVector<double> 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const = 0;
+
   const std::vector<const std::vector<Polynomial>*>
     getLocalFunctions(const MElement& element) const;
 
@@ -64,7 +69,26 @@ class FunctionSpaceVector : public FunctionSpace{
    @param coef The coefficients of the interpolation
    @param xyz The coordinate 
    (of point @em inside the given @c element)
-   of the interpolation
+   of the interpolation in the @em Physical Space
+
+   @return Returns the (vectorial) interpolated value
+
+   @warning
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   
+   @todo
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   ---> check
+   **
+
+   @fn FunctionSpaceVector::interpolateInRefSpace
+   @param element The MElement to interpolate on
+   @param coef The coefficients of the interpolation
+   @param xyz The coordinate 
+   (of point @em inside the given @c element)
+   of the interpolation in the @em Reference Space
 
    @return Returns the (vectorial) interpolated value
 
diff --git a/FunctionSpace/LagrangeBasis.cpp b/FunctionSpace/LagrangeBasis.cpp
index 63b51e1efe49a18193412f0740bb508973948332..9402df5f4cb408d94e3520fb73e4bb0e729ba6ed 100644
--- a/FunctionSpace/LagrangeBasis.cpp
+++ b/FunctionSpace/LagrangeBasis.cpp
@@ -9,24 +9,50 @@ LagrangeBasis::LagrangeBasis(void){
 LagrangeBasis::~LagrangeBasis(void){
 }
 
-fullVector<double> LagrangeBasis::
-project(const fullVector<double>& coef,
-	const std::vector<const Polynomial*>& basis){
+vector<double> LagrangeBasis::
+project(const MElement& element,
+	const vector<double>& coef,
+	const FunctionSpaceScalar& fSpace){
 
   // Init New Coefs //
   const unsigned int size = point->size1();
-  fullVector<double> newCoef(size);
+  vector<double> newCoef(size);
 
   // Interpolation at Lagrange Points //
-  const unsigned int nCoef = coef.size();
- 
   for(unsigned int i = 0; i < size; i++){
-    newCoef(i) = 0;
+    fullVector<double> uvw(3);
+    uvw(0) = (*point)(i, 0);
+    uvw(1) = (*point)(i, 1);
+    uvw(2) = (*point)(i, 2);
+
+    newCoef[i] = fSpace.interpolateInRefSpace(element, 
+					      coef, 
+					      uvw);
+  }
+    
+  // Return ;
+  return newCoef;
+}
+
+vector<fullVector<double> > LagrangeBasis::
+project(const MElement& element,
+	const vector<double>& coef,
+	const FunctionSpaceVector& fSpace){
+  
+  // Init New Coefs //
+  const unsigned int size = point->size1();
+  vector<fullVector<double> > newCoef(size);
+
+  // Interpolation at Lagrange Points //
+  for(unsigned int i = 0; i < size; i++){
+    fullVector<double> uvw(3);
+    uvw(0) = (*point)(i, 0);
+    uvw(1) = (*point)(i, 1);
+    uvw(2) = (*point)(i, 2);
     
-    for(unsigned int j = 0; j < nCoef; j++)
-      newCoef(i) += coef(j) * basis[j]->at((*point)(i, 0),
-					   (*point)(i, 1),
-					   (*point)(i, 2));
+    newCoef[i] = fSpace.interpolateInRefSpace(element, 
+					      coef, 
+					      uvw);
   }
     
   // Return ;
diff --git a/FunctionSpace/LagrangeBasis.h b/FunctionSpace/LagrangeBasis.h
index 0d29cb6c782322d9792899c9e38d9c9755719bc2..7daa7f49fe30bfb1bc20d444899117dd4d90bf63 100644
--- a/FunctionSpace/LagrangeBasis.h
+++ b/FunctionSpace/LagrangeBasis.h
@@ -6,6 +6,10 @@
 #include "polynomialBasis.h"
 #include "BasisScalar.h"
 
+#include "MElement.h"
+#include "FunctionSpaceScalar.h"
+#include "FunctionSpaceVector.h"
+
 /**
    @interface LagrangeBasis
    @brief Interoface for Lagrange Basis
@@ -39,13 +43,28 @@ class LagrangeBasis: public BasisScalar{
   //! @return Returns the Monomial Matrix
   const fullMatrix<double>& getMonomial(void) const;
 
-  //! @param coef A vector of Real numbers
-  //! @param basis A vector of Polynomials 
-  //! in a @em Reference Space
+  //! @param element A MElement
+  //! @param coef A vector of coefficient associated 
+  //! to the given Element
+  //! @param fSpace The (scalar) Function Space 
+  //! of the given Coefficients
   //! @return Projects the given Coefficients in this LagrangeBasis@n
-  fullVector<double> project(const fullVector<double>& coef,
-			     const std::vector<const Polynomial*>& basis);
-
+  std::vector<double> project(const MElement& element,
+			      const std::vector<double>& coef,
+			      const FunctionSpaceScalar& fSpace);
+
+  //! @param element A MElement
+  //! @param coef A vector of coefficient associated 
+  //! to the given Element
+  //! @param fSpace The (vectorial) Function Space 
+  //! of the given Coefficients
+  //! @return Projects the given Coefficients in this LagrangeBasis@n
+  //! @note Each Coefficients will be projected into a vector
+  //! with the same dimesion as the vectorial Polynomials
+  std::vector<fullVector<double> > 
+    project(const MElement& element,
+	    const std::vector<double>& coef,
+	    const FunctionSpaceVector& fSpace);
  protected:
   //! Returns a new LagrangeBasis
   //!