diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index f487391c1216e1ab3bac141869a3a6dc70111b8f..e1e515091f34950e9f7e8e1de80799fa579c94d8 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -81,12 +81,18 @@ class Basis{
   virtual unsigned int getNOrientation(void) const = 0;
   virtual unsigned int getOrientation(const MElement& element) const = 0;
 
+  // Functions Permutation //
+  virtual void getFunctionPermutation(const MElement& element,
+                                      unsigned int* indexPermutation) const = 0;
+
   // Direct Access to Evaluated Functions //
-  virtual fullMatrix<double>* getFunctions(const MElement& element,
-                                           double u, double v, double w) const = 0;
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            const MElement& element,
+                            double u, double v, double w) const = 0;
 
-  virtual fullMatrix<double>* getFunctions(unsigned int orientation,
-                                           double u, double v, double w) const = 0;
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            unsigned int orientation,
+                            double u, double v, double w) const = 0;
 
   // Precompute Functions //
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const = 0;
diff --git a/FunctionSpace/BasisHierarchical0From.cpp b/FunctionSpace/BasisHierarchical0From.cpp
index 27fd1a3cadbdf6600974cfb66d5a84cd38b6ef1b..36778207b728b765f595b9a21e5931742732c55d 100644
--- a/FunctionSpace/BasisHierarchical0From.cpp
+++ b/FunctionSpace/BasisHierarchical0From.cpp
@@ -59,37 +59,32 @@ getOrientation(const MElement& element) const{
   return refSpace->getPermutation(element);
 }
 
-fullMatrix<double>* BasisHierarchical0From::
-getFunctions(const MElement& element,
-	     double u, double v, double w) const{
+void BasisHierarchical0From::
+getFunctionPermutation(const MElement& element,
+                       unsigned int* indexPermutation) const{
+}
 
-  // Alloc Memory //
-  fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
+void BasisHierarchical0From::
+getFunctions(fullMatrix<double>& retValues,
+             const MElement& element,
+             double u, double v, double w) const{
 
   // Define Orientation //
   unsigned int orientation = refSpace->getPermutation(element);
 
   // Fill Matrix //
   for(unsigned int i = 0; i < nFunction; i++)
-    (*values)(i, 0) = basis[orientation][i]->at(u, v, w);
-
-  // Return //
-  return values;
+    retValues(i, 0) = basis[orientation][i]->at(u, v, w);
 }
 
-fullMatrix<double>* BasisHierarchical0From::
-getFunctions(unsigned int orientation,
-	     double u, double v, double w) const{
-
-  // Alloc Memory //
-  fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
+void BasisHierarchical0From::
+getFunctions(fullMatrix<double>& retValues,
+             unsigned int orientation,
+             double u, double v, double w) const{
 
   // Fill Matrix //
   for(unsigned int i = 0; i < nFunction; i++)
-    (*values)(i, 0) = basis[orientation][i]->at(u, v, w);
-
-  // Return //
-  return values;
+    retValues(i, 0) = basis[orientation][i]->at(u, v, w);
 }
 
 void BasisHierarchical0From::
@@ -114,10 +109,10 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
   for(unsigned int i = 0; i < nRefSpace; i++)
     for(unsigned int j = 0; j < nFunction; j++)
       for(unsigned int k = 0; k < nPoint; k++)
-	(*preEvaluatedFunction[i])(j, k) =
-	  basis[i][j]->at(point(k, 0),
-			  point(k, 1),
-			  point(k, 2));
+        (*preEvaluatedFunction[i])(j, k) =
+          basis[i][j]->at(point(k, 0),
+                          point(k, 1),
+                          point(k, 2));
 
   // PreEvaluated //
   preEvaluated = true;
@@ -152,14 +147,14 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
   for(unsigned int i = 0; i < nRefSpace; i++){
     for(unsigned int j = 0; j < nFunction; j++){
       for(unsigned int k = 0; k < nPoint; k++){
-	tmp = Polynomial::at(*grad[i][j],
-			     point(k, 0),
-			     point(k, 1),
-			     point(k, 2));
-
-	(*preEvaluatedGradFunction[i])(j, 3 * k)     = tmp(0);
-	(*preEvaluatedGradFunction[i])(j, 3 * k + 1) = tmp(1);
-	(*preEvaluatedGradFunction[i])(j, 3 * k + 2) = tmp(2);
+        tmp = Polynomial::at(*grad[i][j],
+                             point(k, 0),
+                             point(k, 1),
+                             point(k, 2));
+
+        (*preEvaluatedGradFunction[i])(j, 3 * k)     = tmp(0);
+        (*preEvaluatedGradFunction[i])(j, 3 * k + 1) = tmp(1);
+        (*preEvaluatedGradFunction[i])(j, 3 * k + 2) = tmp(2);
       }
     }
   }
@@ -205,7 +200,7 @@ void BasisHierarchical0From::getGrad(void) const{
   for(unsigned int s = 0; s < nRefSpace; s++)
     for(unsigned int f = 0 ; f < nFunction; f++)
       grad[s][f] =
-	new vector<Polynomial>(basis[s][f]->gradient());
+        new vector<Polynomial>(basis[s][f]->gradient());
 
   // Has Grad //
   hasGrad = true;
@@ -219,22 +214,22 @@ string BasisHierarchical0From::toString(void) const{
   stream << "Vertex Based:" << endl;
   for(; i < nVertex; i++)
     stream << "f("  << i + 1                 << ") = "
-	   << basis[refSpace][i]->toString() << endl;
+           << basis[refSpace][i]->toString() << endl;
 
   stream << "Edge Based:"   << endl;
   for(; i < nVertex + nEdge; i++)
     stream << "f(" << i + 1                  << ") = "
-	   << basis[refSpace][i]->toString() << endl;
+           << basis[refSpace][i]->toString() << endl;
 
   stream << "Face Based:"   << endl;
   for(; i < nVertex + nEdge + nFace; i++)
     stream << "f(" << i + 1                  << ") = "
-	   << basis[refSpace][i]->toString() << endl;
+           << basis[refSpace][i]->toString() << endl;
 
   stream << "Cell Based:"   << endl;
   for(; i < nVertex + nEdge + nFace + nCell; i++)
     stream << "f("  << i + 1                 << ") = "
-	   << basis[refSpace][i]->toString() << endl;
+           << basis[refSpace][i]->toString() << endl;
 
   return stream.str();
 }
diff --git a/FunctionSpace/BasisHierarchical0From.h b/FunctionSpace/BasisHierarchical0From.h
index 6acc773198e2a3aa3a2f8951c305cfc64965a9cf..15c180f28b0690dd5ec0ba58b1de0a5c78e84ced 100644
--- a/FunctionSpace/BasisHierarchical0From.h
+++ b/FunctionSpace/BasisHierarchical0From.h
@@ -39,11 +39,16 @@ class BasisHierarchical0From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual fullMatrix<double>* getFunctions(const MElement& element,
-					   double u, double v, double w) const;
+  virtual void getFunctionPermutation(const MElement& element,
+                                      unsigned int* indexPermutation) const;
 
-  virtual fullMatrix<double>* getFunctions(unsigned int orientation,
-					   double u, double v, double w) const;
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            const MElement& element,
+                            double u, double v, double w) const;
+
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            unsigned int orientation,
+                            double u, double v, double w) const;
 
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
diff --git a/FunctionSpace/BasisHierarchical1From.cpp b/FunctionSpace/BasisHierarchical1From.cpp
index 89e303310794ed81bf48a40e14365240ad6e135c..826d4b64c9dac718f4f1e0425ffb4b59559cfccb 100644
--- a/FunctionSpace/BasisHierarchical1From.cpp
+++ b/FunctionSpace/BasisHierarchical1From.cpp
@@ -25,7 +25,7 @@ BasisHierarchical1From::~BasisHierarchical1From(void){
   if(hasCurl){
     for(unsigned int i = 0; i < nRefSpace; i++){
       for(unsigned int j = 0; j < nFunction; j++)
-	delete curl[i][j];
+        delete curl[i][j];
 
       delete[] curl[i];
     }
@@ -59,12 +59,15 @@ getOrientation(const MElement& element) const{
   return refSpace->getPermutation(element);
 }
 
-fullMatrix<double>* BasisHierarchical1From::
-getFunctions(const MElement& element,
-	     double u, double v, double w) const{
+void BasisHierarchical1From::
+getFunctionPermutation(const MElement& element,
+                       unsigned int* indexPermutation) const{
+}
 
-  // Alloc Memory //
-  fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
+void BasisHierarchical1From::
+getFunctions(fullMatrix<double>& retValues,
+             const MElement& element,
+             double u, double v, double w) const{
 
   // Define Orientation //
   unsigned int orientation = refSpace->getPermutation(element);
@@ -74,34 +77,26 @@ getFunctions(const MElement& element,
     fullVector<double> eval =
       Polynomial::at(*basis[orientation][i], u, v, w);
 
-    (*values)(i, 0) = eval(0);
-    (*values)(i, 1) = eval(1);
-    (*values)(i, 2) = eval(2);
+    retValues(i, 0) = eval(0);
+    retValues(i, 1) = eval(1);
+    retValues(i, 2) = eval(2);
   }
-
-  // Return //
-  return values;
 }
 
-fullMatrix<double>* BasisHierarchical1From::
-getFunctions(unsigned int orientation,
-	     double u, double v, double w) const{
-
-  // Alloc Memory //
-  fullMatrix<double>* values = new fullMatrix<double>(nFunction, 3);
+void BasisHierarchical1From::
+getFunctions(fullMatrix<double>& retValues,
+             unsigned int orientation,
+             double u, double v, double w) const{
 
   // Fill Vector //
   for(unsigned int i = 0; i < nFunction; i++){
     fullVector<double> eval =
       Polynomial::at(*basis[orientation][i], u, v, w);
 
-    (*values)(i, 0) = eval(0);
-    (*values)(i, 1) = eval(1);
-    (*values)(i, 2) = eval(2);
+    retValues(i, 0) = eval(0);
+    retValues(i, 1) = eval(1);
+    retValues(i, 2) = eval(2);
   }
-
-  // Return //
-  return values;
 }
 
 void BasisHierarchical1From::
@@ -129,14 +124,14 @@ preEvaluateFunctions(const fullMatrix<double>& point) const{
   for(unsigned int i = 0; i < nRefSpace; i++){
     for(unsigned int j = 0; j < nFunction; j++){
       for(unsigned int k = 0; k < nPoint; k++){
-	tmp = Polynomial::at(*basis[i][j],
-			     point(k, 0),
-			     point(k, 1),
-			     point(k, 2));
-
-	(*preEvaluatedFunction[i])(j, 3 * k)     = tmp(0);
-	(*preEvaluatedFunction[i])(j, 3 * k + 1) = tmp(1);
-	(*preEvaluatedFunction[i])(j, 3 * k + 2) = tmp(2);
+        tmp = Polynomial::at(*basis[i][j],
+                             point(k, 0),
+                             point(k, 1),
+                             point(k, 2));
+
+        (*preEvaluatedFunction[i])(j, 3 * k)     = tmp(0);
+        (*preEvaluatedFunction[i])(j, 3 * k + 1) = tmp(1);
+        (*preEvaluatedFunction[i])(j, 3 * k + 2) = tmp(2);
       }
     }
   }
@@ -174,14 +169,14 @@ preEvaluateDerivatives(const fullMatrix<double>& point) const{
   for(unsigned int i = 0; i < nRefSpace; i++){
     for(unsigned int j = 0; j < nFunction; j++){
       for(unsigned int k = 0; k < nPoint; k++){
-	tmp = Polynomial::at(*curl[i][j],
-			     point(k, 0),
-			     point(k, 1),
-			     point(k, 2));
-
-	(*preEvaluatedCurlFunction[i])(j, 3 * k)     = tmp(0);
-	(*preEvaluatedCurlFunction[i])(j, 3 * k + 1) = tmp(1);
-	(*preEvaluatedCurlFunction[i])(j, 3 * k + 2) = tmp(2);
+        tmp = Polynomial::at(*curl[i][j],
+                             point(k, 0),
+                             point(k, 1),
+                             point(k, 2));
+
+        (*preEvaluatedCurlFunction[i])(j, 3 * k)     = tmp(0);
+        (*preEvaluatedCurlFunction[i])(j, 3 * k + 1) = tmp(1);
+        (*preEvaluatedCurlFunction[i])(j, 3 * k + 2) = tmp(2);
       }
     }
   }
@@ -227,7 +222,7 @@ void BasisHierarchical1From::getCurl(void) const{
   for(unsigned int s = 0; s < nRefSpace; s++)
     for(unsigned int f = 0 ; f < nFunction; f++)
       curl[s][f] =
-	new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
+        new vector<Polynomial>(Polynomial::curl(*basis[s][f]));
 
   // Has Curl //
   hasCurl = true;
@@ -241,30 +236,30 @@ string BasisHierarchical1From::toString(void) const{
   stream << "Vertex Based:" << endl;
   for(; i < nVertex; i++)
     stream << "f("   << i + 1                               << ") = " << endl
-	   << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
+           << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
 
   stream << "Edge Based:"   << endl;
   for(; i < nVertex + nEdge; i++)
     stream << " f("  << i + 1                               << ") = " << endl
-	   << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
+           << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
 
   stream << "Face Based:"   << endl;
   for(; i < nVertex + nEdge + nFace; i++)
     stream << " f("  << i + 1                               << ") = " << endl
-	   << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
+           << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
 
   stream << "Cell Based:"   << endl;
   for(; i < nVertex + nEdge + nFace + nCell; i++)
     stream << " f("  << i + 1                               << ") = " << endl
-	   << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
-	   << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
+           << "\t[ " << (*basis[refSpace][i])[0].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[1].toString() << " ]"   << endl
+           << "\t[ " << (*basis[refSpace][i])[2].toString() << " ]"   << endl;
 
   return stream.str();
 }
diff --git a/FunctionSpace/BasisHierarchical1From.h b/FunctionSpace/BasisHierarchical1From.h
index 82faf8005b682e845c4a1ec8c062e849257b06ab..ad0f99e352086b8f614323cafcedd306d9de0c15 100644
--- a/FunctionSpace/BasisHierarchical1From.h
+++ b/FunctionSpace/BasisHierarchical1From.h
@@ -39,11 +39,16 @@ class BasisHierarchical1From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual fullMatrix<double>* getFunctions(const MElement& element,
-					   double u, double v, double w) const;
+  virtual void getFunctionPermutation(const MElement& element,
+                                      unsigned int* indexPermutation) const;
 
-  virtual fullMatrix<double>* getFunctions(unsigned int orientation,
-					   double u, double v, double w) const;
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            const MElement& element,
+                            double u, double v, double w) const;
+
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            unsigned int orientation,
+                            double u, double v, double w) const;
 
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
@@ -88,4 +93,3 @@ class BasisHierarchical1From: public BasisLocal{
  */
 
 #endif
-
diff --git a/FunctionSpace/BasisLagrange.cpp b/FunctionSpace/BasisLagrange.cpp
index 12656d71dc9cc3ec01e1c7990e8cd80fb378547a..36a0d8732972fc0b7dcf49adc12115838c8dcb78 100644
--- a/FunctionSpace/BasisLagrange.cpp
+++ b/FunctionSpace/BasisLagrange.cpp
@@ -12,16 +12,22 @@ BasisLagrange::~BasisLagrange(void){
 
 unsigned int BasisLagrange::
 getNOrientation(void) const{
-  return refSpace->getNPermutation();
+  return 1;
 }
 
 unsigned int BasisLagrange::
 getOrientation(const MElement& element) const{
-  return refSpace->getPermutation(element);
+  return 0;
 }
 
-fullMatrix<double>* BasisLagrange::
-getFunctions(const MElement& element,
+void BasisLagrange::
+getFunctionPermutation(const MElement& element,
+                       unsigned int* indexPermutation) const{
+}
+
+void BasisLagrange::
+getFunctions(fullMatrix<double>& retValues,
+             const MElement& element,
              double u, double v, double w) const{
 
   // Fill Matrix //
@@ -34,14 +40,12 @@ getFunctions(const MElement& element,
   lBasis->f(point, tmp);
 
   // Transpose 'tmp': otherwise not coherent with df !!
-  fullMatrix<double> values = tmp.transpose();
-
-  // Get Inorder & Return //
-  return inorder(refSpace->getPermutation(element), values);
+  retValues = tmp.transpose();
 }
 
-fullMatrix<double>* BasisLagrange::
-getFunctions(unsigned int orientation,
+void BasisLagrange::
+getFunctions(fullMatrix<double>& retValues,
+             unsigned int orientation,
              double u, double v, double w) const{
 
   // Fill Matrix //
@@ -54,10 +58,7 @@ getFunctions(unsigned int orientation,
   lBasis->f(point, tmp);
 
   // Transpose 'tmp': otherwise not coherent with df !!
-  fullMatrix<double> values = tmp.transpose();
-
-  // Get Inorder & Return //
-  return inorder(orientation, values);
+  retValues = tmp.transpose();
 }
 
 void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
@@ -164,25 +165,3 @@ project(const MElement& element,
   // Return ;
   return newCoef;
 }
-
-fullMatrix<double>* BasisLagrange::inorder(unsigned int orientation,
-                                           fullMatrix<double>& mat) const{
-  // Matrix Size //
-  const int nRow = mat.size1();
-  const int nCol = mat.size2();
-
-  // Order //
-  const vector<unsigned int>* order =
-    refSpace->getAllLagrangeNode()[orientation];
-
-  // Alloc new matrix //
-  fullMatrix<double>* ret = new fullMatrix<double>(nRow, nCol);
-
-  // Populate //
-  for(int i = 0; i < nRow; i++)
-    for(int j = 0; j < nCol; j++)
-      (*ret)(i, j) = mat((*order)[i], j);
-
-  // Return //
-  return ret;
-}
diff --git a/FunctionSpace/BasisLagrange.h b/FunctionSpace/BasisLagrange.h
index 37d85d31e34d7df1569df87f830fd580055fa32d..105802265c1436a7f55fe4c922135257922c7adf 100644
--- a/FunctionSpace/BasisLagrange.h
+++ b/FunctionSpace/BasisLagrange.h
@@ -36,11 +36,16 @@ class BasisLagrange: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
-  virtual fullMatrix<double>* getFunctions(const MElement& element,
-                                           double u, double v, double w) const;
+  virtual void getFunctionPermutation(const MElement& element,
+                                      unsigned int* indexPermutation) const;
 
-  virtual fullMatrix<double>* getFunctions(unsigned int orientation,
-                                           double u, double v, double w) const;
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            const MElement& element,
+                            double u, double v, double w) const;
+
+  virtual void getFunctions(fullMatrix<double>& retValues,
+                            unsigned int orientation,
+                            double u, double v, double w) const;
 
   virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
   virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
@@ -72,9 +77,6 @@ class BasisLagrange: public BasisLocal{
 
  protected:
   BasisLagrange(void);
-
-  fullMatrix<double>* inorder(unsigned int orientation,
-                              fullMatrix<double>& mat) const;
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 42eb4218cb8e1afa0e1b789247df94a49a032c04..fa675e8c97eeec68d0cab8c717576123fdc6c609 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -26,19 +26,18 @@ interpolate(const MElement& element,
   eelement.xyz2uvw(phys, uvw);
 
   // Get Basis Functions //
-  fullMatrix<double>* fun =
-    (*basis)[0]->getFunctions(element, uvw[0], uvw[1], uvw[2]);
+  const unsigned int nFun = (*basis)[0]->getNFunction();
+  fullMatrix<double>  fun(nFun, 1);
 
-  const unsigned int nFun = fun->size1();
+  (*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
 
   // Interpolate (in Reference Place) //
   double val = 0;
 
   for(unsigned int i = 0; i < nFun; i++)
-    val += (*fun)(i, 0) * coef[i];
+    val += fun(i, 0) * coef[i];
 
   // Return Interpolated Value //
-  delete fun;
   return val;
 }
 
@@ -48,18 +47,17 @@ interpolateInRefSpace(const MElement& element,
                       const fullVector<double>& uvw) const{
 
   // Get Basis Functions //
-  fullMatrix<double>* fun =
-    (*basis)[0]->getFunctions(element, uvw(0), uvw(1), uvw(2));
+  const unsigned int nFun = (*basis)[0]->getNFunction();
+  fullMatrix<double>  fun(nFun, 1);
 
-  const unsigned int nFun = fun->size1();
+  (*basis)[0]->getFunctions(fun, element, uvw(0), uvw(1), uvw(2));
 
   // Interpolate (in Reference Place) //
   double val = 0;
 
   for(unsigned int i = 0; i < nFun; i++)
-    val += (*fun)(i, 0) * coef[i];
+    val += fun(i, 0) * coef[i];
 
   // Return Interpolated Value //
-  delete fun;
   return val;
 }
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index ce4b89207bf9a2a392174957c4ee71ba3bc7cbaa..cc7fab61a605713e48dab5087e2067fdb242c739 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -32,10 +32,10 @@ interpolate(const MElement& element,
   invJac.invertInPlace();
 
   // Get Basis Functions //
-  fullMatrix<double>* fun =
-    (*basis)[0]->getFunctions(element, uvw[0], uvw[1], uvw[2]);
+  const unsigned int nFun = (*basis)[0]->getNFunction();
+  fullMatrix<double>  fun(nFun, 3);
 
-  const unsigned int nFun = fun->size1();
+  (*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
 
   // Interpolate (in Reference Place) //
   fullMatrix<double> val(1, 3);
@@ -44,14 +44,12 @@ interpolate(const MElement& element,
   val(0, 2) = 0;
 
   for(unsigned int i = 0; i < nFun; i++){
-    val(0, 0) += (*fun)(i, 0) * coef[i];
-    val(0, 1) += (*fun)(i, 1) * coef[i];
-    val(0, 2) += (*fun)(i, 2) * coef[i];
+    val(0, 0) += fun(i, 0) * coef[i];
+    val(0, 1) += fun(i, 1) * coef[i];
+    val(0, 2) += fun(i, 2) * coef[i];
   }
 
   // Return Interpolated Value //
-  delete fun;
-
   fullVector<double> map(3);
   Mapper::hCurl(val, 0, 0, invJac, map);
   return map;
@@ -72,11 +70,10 @@ interpolateInRefSpace(const MElement& element,
   invJac.invertInPlace();
 
   // Get Basis Functions //
-  fullMatrix<double>* fun =
-    (*basis)[0]->getFunctions(element, uvw(0), uvw(1), uvw(2));
-
-  const unsigned int nFun = fun->size1();
+  const unsigned int nFun = (*basis)[0]->getNFunction();
+  fullMatrix<double>  fun(nFun, 3);
 
+  (*basis)[0]->getFunctions(fun, element, uvw(0), uvw(1), uvw(2));
 
   // Interpolate (in Reference Place) //
   fullMatrix<double> val(1, 3);
@@ -85,14 +82,12 @@ interpolateInRefSpace(const MElement& element,
   val(0, 2) = 0;
 
   for(unsigned int i = 0; i < nFun; i++){
-    val(0, 0) += (*fun)(i, 0) * coef[i];
-    val(0, 1) += (*fun)(i, 1) * coef[i];
-    val(0, 2) += (*fun)(i, 2) * coef[i];
+    val(0, 0) += fun(i, 0) * coef[i];
+    val(0, 1) += fun(i, 1) * coef[i];
+    val(0, 2) += fun(i, 2) * coef[i];
   }
 
   // Return Interpolated Value //
-  delete fun;
-
   fullVector<double> map(3);
   Mapper::hCurl(val, 0, 0, invJac, map);
   return map;