diff --git a/FunctionSpace/BasisTest.cpp b/FunctionSpace/BasisTest.cpp
index a76ff1081c13d016fa37727094f51a5c495046d5..24cebdb53564e37d6a8fc016db7b75a6afb51bde 100644
--- a/FunctionSpace/BasisTest.cpp
+++ b/FunctionSpace/BasisTest.cpp
@@ -25,7 +25,8 @@ int ain(int argc, char** argv){
   Mesh msh(argv[1]);
 
   // Writer for .msh
-  WriterMsh writer(msh.getGroup(2).getAll());
+  WriterMsh writer;
+  writer.setDomain(msh.getGroup(2).getAll());
 
   // Plot Basis //
   HexNodeBasis b(1);
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 4f0800d708068646ff2937dbda26e444c57a6d9a..8a0101e828353f859e20efccc5330bad125e6bbb 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -22,6 +22,7 @@ set(SRC
   HexEdgeBasis.cpp
 
   FunctionSpace.cpp
+  FunctionSpaceScalar.cpp
   FunctionSpaceNode.cpp
 )
 
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index dce20c1b6df2a0763e3bae62c9a229d8d356a8b4..2fa2be3b4f748ade4b2780f1644ed7df45b97896 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/FunctionSpace/FunctionSpaceNode.cpp
@@ -1,12 +1,14 @@
 #include "FunctionSpaceNode.h"
 #include "Polynomial.h"
 #include "BasisScalar.h"
+#include "fullMatrix.h"
+#include "Mapper.h"
 #include "Exception.h"
 
 using namespace std;
 
 FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe, int order):
-FunctionSpace(goe, 0, order){
+FunctionSpaceScalar(goe, 0, order){
   // Just Call Super Constructor with basisType = 0
 }
 
@@ -14,9 +16,10 @@ FunctionSpace(goe, 0, order){
 FunctionSpaceNode::~FunctionSpaceNode(void){
 }
 
-vector<double> FunctionSpaceNode::
+void FunctionSpaceNode::
 interpolateAtNodes(const MElement& elem, 
-		   const vector<double>& coef) const{
+		   const vector<double>& coef,
+		   std::vector<double>& nodeValue) const{
   // Check
   unsigned int wS = coef.size();
   unsigned int bS = basis->getSize();
@@ -42,20 +45,33 @@ interpolateAtNodes(const MElement& elem,
   const vector<Polynomial>& fun = 
     static_cast<const BasisScalar*>(basis)->getBasis();
 
-  // Init Vector
-  vector<double> value(N);
+  // Init some stuff
+  fullMatrix<double> invJac(3, 3);
+  fullVector<double> xyz(3);
+  fullVector<double> origin(3);
+  
+  origin(0) = node[0]->x();
+  origin(1) = node[0]->y();
+  origin(2) = node[0]->z();
 
   // Interpolate
   for(unsigned int n = 0; n < N; n++){
-    const double x = node[n]->x();
-    const double y = node[n]->y();
-    const double z = node[n]->z();
+    // Map from physical to reference space 
+    xyz(0) = node[n]->x();
+    xyz(1) = node[n]->y();
+    xyz(2) = node[n]->z();
 
-    value[n] = 0;
+    element.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
+    invJac.invertInPlace();
+    
+    const fullVector<double> uvw = 
+      Mapper::invMap(xyz, origin, invJac);
+
+
+    // Interpolate
+    const int id = node[n]->getNum();
 
     for(unsigned int i = 0; i < bS; i++)
-      value[n] += fun[i].at(x, y, z) * coef[i];
+      nodeValue[id] += fun[i].at(uvw(0), uvw(1), uvw(2)) * coef[i];
   }
-
-  return value;
 }
diff --git a/FunctionSpace/FunctionSpaceNode.h b/FunctionSpace/FunctionSpaceNode.h
index d63400836da21e27d33a7dcffff40928fd729990..0af13f213e2786ae9101d848365bf4aebcf5ebbd 100644
--- a/FunctionSpace/FunctionSpaceNode.h
+++ b/FunctionSpace/FunctionSpaceNode.h
@@ -1,15 +1,17 @@
 #ifndef _FUNCTIONSPACENODE_H_
 #define _FUNCTIONSPACENODE_H_
 
-#include "FunctionSpace.h"
+#include "FunctionSpaceScalar.h"
 
-class FunctionSpaceNode: public FunctionSpace{
+class FunctionSpaceNode: public FunctionSpaceScalar{
+ public:
   FunctionSpaceNode(const GroupOfElement& goe, int order);
-
+  
   virtual ~FunctionSpaceNode(void);
-
-  std::vector<double> interpolateAtNodes(const MElement& element, 
-					 const std::vector<double>& coef) const;
+  
+  virtual void interpolateAtNodes(const MElement& element, 
+				  const std::vector<double>& coef,
+				  std::vector<double>& nodeValue) const;
 };
 
 #endif
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..43d737cde5657fff320d60f6469f2a07e968432e
--- /dev/null
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -0,0 +1,9 @@
+#include "FunctionSpaceScalar.h"
+
+FunctionSpaceScalar::FunctionSpaceScalar(const GroupOfElement& goe, 
+					  int basisType, int order):
+FunctionSpace(goe, basisType, order){
+}
+
+FunctionSpaceScalar::~FunctionSpaceScalar(void){
+}
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
new file mode 100644
index 0000000000000000000000000000000000000000..e8c09f65ff087a56b7c0d80caa7f894a9406f6c4
--- /dev/null
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -0,0 +1,18 @@
+#ifndef _FUNCTIONSPACESCALAR_H_
+#define _FUNCTIONSPACESCALAR_H_
+
+#include "FunctionSpace.h"
+
+class FunctionSpaceScalar: public FunctionSpace{
+ public:
+  FunctionSpaceScalar(const GroupOfElement& goe, 
+		      int basisType, int order);
+
+  virtual ~FunctionSpaceScalar(void);
+
+  virtual void interpolateAtNodes(const MElement& element, 
+				  const std::vector<double>& coef,
+				  std::vector<double>& nodeValue) const = 0;
+};
+
+#endif