diff --git a/FunctionSpace/HexEdgeBasis.cpp b/FunctionSpace/HexEdgeBasis.cpp
index 7ebbbdee4101e8be34135ecf1fff54d4a65fb1b4..34f607aab48f4e01ba67df137da5507220448ab6 100644
--- a/FunctionSpace/HexEdgeBasis.cpp
+++ b/FunctionSpace/HexEdgeBasis.cpp
@@ -8,7 +8,7 @@ HexEdgeBasis::HexEdgeBasis(int order){
   /*
   // Set Basis Type //
   this->order = order;
-  
+
   type = 1;
   dim  = 3;
 
@@ -60,43 +60,43 @@ HexEdgeBasis::HexEdgeBasis(int order){
   int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
 
 
-  // Lagrange // 
-  lagrange[0] = 
+  // Lagrange //
+  lagrange[0] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lagrange[1] = 
+  lagrange[1] =
     (Polynomial(1, 1, 0, 0))                          *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lagrange[2] = 
+  lagrange[2] =
     (Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 1, 0)) *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lagrange[3] = 
+  lagrange[3] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 1, 0))                          *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lagrange[4] = 
+  lagrange[4] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
      Polynomial(1, 0, 0, 1);
 
-  lagrange[5] = 
+  lagrange[5] =
     (Polynomial(1, 1, 0, 0))                          *
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
      Polynomial(1, 0, 0, 1);
 
-  lagrange[6] = 
+  lagrange[6] =
     (Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 1, 0)) *
      Polynomial(1, 0, 0, 1);
 
-  lagrange[7] = 
+  lagrange[7] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
     (Polynomial(1, 0, 1, 0))                          *
      Polynomial(1, 0, 0, 1);
@@ -105,44 +105,44 @@ HexEdgeBasis::HexEdgeBasis(int order){
   for(int i = 0; i < 12; i++)
     lagrangeSum[i] = lagrange[edge1[i]] + lagrange[edge2[i]];
 
-    
+
   // Lifting //
-  lifting[0] = 
+  lifting[0] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[1] = 
+  lifting[1] =
     (Polynomial(1, 1, 0, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[2] = 
+  lifting[2] =
     (Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[3] = 
+  lifting[3] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[4] = 
+  lifting[4] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[5] = 
+  lifting[5] =
     (Polynomial(1, 1, 0, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[6] = 
+  lifting[6] =
     (Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[7] = 
+  lifting[7] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0))                          +
      Polynomial(1, 0, 0, 1);
@@ -156,14 +156,14 @@ HexEdgeBasis::HexEdgeBasis(int order){
   std::vector<std::vector<Polynomial>*> basis(size);
 
 
-  // Edge Based (Nedelec) // 
+  // Edge Based (Nedelec) //
   int i = 0; // Function Counter
   Polynomial oneHalf(0.5, 0, 0, 0);
 
   for(int e = 0; e < 12; e++){
-    basis[i] = 
+    basis[i] =
       new std::vector<Polynomial>((liftingSub[e]).gradient());
-    
+
     basis[i]->at(0).mul(lagrangeSum[e]);
     basis[i]->at(1).mul(lagrangeSum[e]);
     basis[i]->at(2).mul(lagrangeSum[e]);
@@ -179,9 +179,9 @@ HexEdgeBasis::HexEdgeBasis(int order){
   // Edge Based (High Order) //
   for(int l = 1; l < orderPlus; l++){
     for(int e = 0; e < 12; e++){
-      basis[i] = 
+      basis[i] =
 	new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
-     
+
       i++;
     }
   }
@@ -204,12 +204,12 @@ HexEdgeBasis::HexEdgeBasis(int order){
 
   // 'Lambda' Functions
   for(int f = 0; f < 6; f++)
-    lambda[f] = 
+    lambda[f] =
       lagrange[face1[f]] +
       lagrange[face2[f]] +
       lagrange[face3[f]] +
       lagrange[face4[f]];
-  
+
   // Gradients
   for(int f = 0; f < 6; f++){
     grXi[f]  = xi[f].gradient();
@@ -219,26 +219,26 @@ HexEdgeBasis::HexEdgeBasis(int order){
   // Compositions
   for(int l = 0; l < orderPlus; l++){
     iLegendreXi[l]  = new Polynomial[6];
-    iLegendreEta[l] = new Polynomial[6]; 
- 
-    legendreXi[l]  = new Polynomial[6];   
+    iLegendreEta[l] = new Polynomial[6];
+
+    legendreXi[l]  = new Polynomial[6];
     legendreEta[l] = new Polynomial[6];
 
     for(int f = 0; f < 6; f++){
       iLegendreXi[l][f]  = intLegendre[l].compose(xi[f]);
       iLegendreEta[l][f] = intLegendre[l].compose(eta[f]);
-      
+
       legendreXi[l][f]   = legendre[l].compose(xi[f]);
       legendreEta[l][f]  = legendre[l].compose(eta[f]);
     }
   }
- 
+
 
   // Face Based (Type 1) //
   for(int l1 = 1; l1 < orderPlus; l1++){
     for(int l2 = 1; l2 < orderPlus; l2++){
       for(int f = 0; f < 6; f++){
-	basis[i] = new std::vector<Polynomial>((iLegendreXi[l1][f]   * 
+	basis[i] = new std::vector<Polynomial>((iLegendreXi[l1][f]   *
 						iLegendreEta[l2][f]  *
 						lambda[f]).gradient());
 
@@ -254,10 +254,10 @@ HexEdgeBasis::HexEdgeBasis(int order){
       for(int f = 0; f < 6; f++){
 	basis[i] = new std::vector<Polynomial>(3);
 
-	Polynomial tmp1 = 
-	  legendreXi[l1][f]   * 
+	Polynomial tmp1 =
+	  legendreXi[l1][f]   *
 	  iLegendreEta[l2][f];
-		       
+
 	Polynomial tmp2 =
 	  iLegendreXi[l1][f] *
 	  legendreEta[l2][f];
@@ -270,8 +270,8 @@ HexEdgeBasis::HexEdgeBasis(int order){
 	vector<Polynomial> gr2 = grEta[f];
 	gr2[0].mul(tmp2);
 	gr2[1].mul(tmp2);
-	gr2[2].mul(tmp2);	
-	
+	gr2[2].mul(tmp2);
+
 	basis[i]->at(0) = (gr1[0] - gr2[0]) * lambda[f];
 	basis[i]->at(1) = (gr1[1] - gr2[1]) * lambda[f];
 	basis[i]->at(2) = (gr1[2] - gr2[2]) * lambda[f];
@@ -286,14 +286,14 @@ HexEdgeBasis::HexEdgeBasis(int order){
   for(int l = 1; l < orderPlus; l++){
     for(int f = 0; f < 6; f++){
       Polynomial tmp = iLegendreEta[l][f] * lambda[f];
-      
+
       basis[i] = new std::vector<Polynomial>(grXi[f]);
-      
+
       basis[i]->at(0).mul(tmp);
       basis[i]->at(1).mul(tmp);
       basis[i]->at(2).mul(tmp);
-      
-      i++;   
+
+      i++;
     }
   }
 
@@ -302,24 +302,24 @@ HexEdgeBasis::HexEdgeBasis(int order){
   for(int l = 1; l < orderPlus; l++){
     for(int f = 0; f < 6; f++){
       Polynomial tmp = iLegendreXi[l][f] * lambda[f];
-      
+
       basis[i] = new std::vector<Polynomial>(grEta[f]);
-      
+
       basis[i]->at(0).mul(tmp);
       basis[i]->at(1).mul(tmp);
       basis[i]->at(2).mul(tmp);
-      
-      i++;   
+
+      i++;
     }
   }
 
-  
+
   // Cell Based (Preliminary) //
   Polynomial px   = Polynomial(2, 1, 0, 0);
   Polynomial py   = Polynomial(2, 0, 1, 0);
   Polynomial pz   = Polynomial(2, 0, 0, 1);
   Polynomial zero = Polynomial(0, 0, 0, 0);
- 
+
   px = px - Polynomial(1, 0, 0, 0);
   py = py - Polynomial(1, 0, 0, 0);
   pz = pz - Polynomial(1, 0, 0, 0);
@@ -338,7 +338,7 @@ HexEdgeBasis::HexEdgeBasis(int order){
   for(int l1 = 1; l1 < orderPlus; l1++){
     for(int l2 = 1; l2 < orderPlus; l2++){
       for(int l3 = 1; l3 < orderPlus; l3++){
-	basis[i] = new std::vector<Polynomial>((iLegendreX[l1] * 
+	basis[i] = new std::vector<Polynomial>((iLegendreX[l1] *
 						iLegendreY[l2] *
 						iLegendreZ[l3]).gradient());
 
@@ -347,16 +347,16 @@ HexEdgeBasis::HexEdgeBasis(int order){
       }
     }
   }
-  
+
 
   // Cell Based (Type 2 -- First Part) //
   for(int j = 0; j < cellNumber; j++){
     basis[i] = new std::vector<Polynomial>(3);
 
     int off = j + cellStart;
- 
+
     basis[i]->at(0) = basis[off]->at(0);
-    basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0); 
+    basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
     basis[i]->at(2) = basis[off]->at(2);
 
     i++;
@@ -370,7 +370,7 @@ HexEdgeBasis::HexEdgeBasis(int order){
     int off = j + cellStart;
 
     basis[i]->at(0) = basis[off]->at(0);
-    basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0); 
+    basis[i]->at(1) = basis[off]->at(1) * Polynomial(-1, 0, 0, 0);
     basis[i]->at(2) = basis[off]->at(2) * Polynomial(-1, 0, 0, 0);
 
     i++;
@@ -379,13 +379,13 @@ HexEdgeBasis::HexEdgeBasis(int order){
 
   // Cell Based (Type 3 -- First Part) //
   for(int l2 = 1; l2 < orderPlus; l2++){
-    for(int l3 = 1; l3 < orderPlus; l3++){ 
+    for(int l3 = 1; l3 < orderPlus; l3++){
       basis[i] = new std::vector<Polynomial>(3);
 
       basis[i]->at(0) = iLegendreY[l2] * iLegendreZ[l3];
       basis[i]->at(1) = zero;
       basis[i]->at(2) = zero;
-      
+
       i++;
     }
   }
@@ -393,30 +393,30 @@ HexEdgeBasis::HexEdgeBasis(int order){
 
   // Cell Based (Type 3 -- Second Part) //
   for(int l1 = 1; l1 < orderPlus; l1++){
-    for(int l3 = 1; l3 < orderPlus; l3++){      
+    for(int l3 = 1; l3 < orderPlus; l3++){
       basis[i] = new std::vector<Polynomial>(3);
-      
+
       basis[i]->at(0) = zero;
       basis[i]->at(1) = iLegendreX[l1] * iLegendreZ[l3];
       basis[i]->at(2) = zero;
-      
+
       i++;
     }
-  }  
+  }
 
 
   // Cell Based (Type 3 -- Thrid Part) //
   for(int l1 = 1; l1 < orderPlus; l1++){
     for(int l2 = 1; l2 < orderPlus; l2++){
       basis[i] = new std::vector<Polynomial>(3);
-      
+
       basis[i]->at(0) = zero;
       basis[i]->at(1) = zero;
       basis[i]->at(2) = iLegendreX[l1] * iLegendreY[l2];
-      
+
       i++;
     }
-  }  
+  }
 
 
   // Free Temporary Sapce //
@@ -438,21 +438,21 @@ HexEdgeBasis::HexEdgeBasis(int order){
 
   for(int l = 0; l < orderPlus; l++){
     delete[] iLegendreXi[l];
-    delete[] iLegendreEta[l]; 
-    delete[] legendreXi[l];   
+    delete[] iLegendreEta[l];
+    delete[] legendreXi[l];
     delete[] legendreEta[l];
   }
 
   delete[] iLegendreXi;
-  delete[] iLegendreEta; 
-  delete[] legendreXi;   
+  delete[] iLegendreEta;
+  delete[] legendreXi;
   delete[] legendreEta;
 
   delete[] iLegendreX;
   delete[] iLegendreY;
   delete[] iLegendreZ;
 
-  
+
   // Set Basis //
   this->basis = new std::vector<const std::vector<Polynomial>*>
     (basis.begin(), basis.end());
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 5cc68e7f0cc0401270948fc6a9edcb035dee1375..6bb1691e6366aaad091b167b2366df5d362f39c4 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -33,16 +33,16 @@ HexNodeBasis::HexNodeBasis(int order){
   Legendre::integrated(legendre, order);
 
   // Vertices definig Edges & Permutations //
-  const int edgeV[2][12][2] = 
+  const int edgeV[2][12][2] =
     {
-      { {0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 5}, {2, 3}, 
+      { {0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 5}, {2, 3},
 	{2, 6},	{3, 7}, {4, 5}, {4, 7}, {5, 6}, {6, 7} },
 
-      { {1, 0}, {3, 0}, {4, 0}, {2, 1}, {5, 1}, {3, 2}, 
+      { {1, 0}, {3, 0}, {4, 0}, {2, 1}, {5, 1}, {3, 2},
 	{6, 2},	{7, 3}, {5, 4}, {7, 4}, {6, 5}, {7, 6} },
     };
 
-  const int faceV[6][4] = 
+  const int faceV[6][4] =
     {
       {0, 3, 2, 1},
       {0, 1, 5, 4},
@@ -52,44 +52,44 @@ HexNodeBasis::HexNodeBasis(int order){
       {4, 5, 6, 7}
     };
 
-  
+
   // Lifting //
-  lifting[0] = 
+  lifting[0] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[1] = 
+  lifting[1] =
     (Polynomial(1, 1, 0, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[2] = 
+  lifting[2] =
     (Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[3] = 
+  lifting[3] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
 
-  lifting[4] = 
+  lifting[4] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[5] = 
+  lifting[5] =
     (Polynomial(1, 1, 0, 0))                          +
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[6] = 
+  lifting[6] =
     (Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0)) +
      Polynomial(1, 0, 0, 1);
 
-  lifting[7] = 
+  lifting[7] =
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
     (Polynomial(1, 0, 1, 0))                          +
      Polynomial(1, 0, 0, 1);
@@ -99,48 +99,48 @@ HexNodeBasis::HexNodeBasis(int order){
   basis = new std::vector<const Polynomial*>(size);
 
 
-  // Vertex Based (Lagrange) // 
-  (*basis)[0] = 
+  // Vertex Based (Lagrange) //
+  (*basis)[0] =
     new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
 
-  (*basis)[1] = 
+  (*basis)[1] =
     new Polynomial((Polynomial(1, 1, 0, 0))                          *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
 
-  (*basis)[2] = 
+  (*basis)[2] =
     new Polynomial((Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 1, 0)) *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
 
-  (*basis)[3] = 
+  (*basis)[3] =
     new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 1, 0))                          *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1)));
 
-  (*basis)[4] = 
+  (*basis)[4] =
     new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
 		   Polynomial(1, 0, 0, 1));
 
-  (*basis)[5] = 
+  (*basis)[5] =
     new Polynomial((Polynomial(1, 1, 0, 0))                          *
 		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
 		   Polynomial(1, 0, 0, 1));
 
-  (*basis)[6] = 
+  (*basis)[6] =
     new Polynomial((Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 1, 0)) *
 		   Polynomial(1, 0, 0, 1));
 
-  (*basis)[7] = 
+  (*basis)[7] =
     new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
 		   (Polynomial(1, 0, 1, 0))                          *
 		   Polynomial(1, 0, 0, 1));
 
-  
+
   // Edge Based //
   // Keep counting
   int i = 8;
@@ -152,14 +152,14 @@ HexNodeBasis::HexNodeBasis(int order){
   for(int l = 1; l < order; l++){
     for(int e = 0; e < 12; e++){
       (*basis)[i] = new Polynomial(
-	legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) * 
+	legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) *
 	(*(*basis)[edge1[e]] + *(*basis)[edge2[e]]));
-            
+
       i++;
     }
   }
 
-  
+
   // Face Based (Preliminary) //
   // Points definig Faces
   const int face1[6] = {0, 3, 2, 1, 5, 4};
@@ -177,7 +177,7 @@ HexNodeBasis::HexNodeBasis(int order){
 
   // 'Lambda' Functions
   for(int f = 0; f < 6; f++)
-    lambda[f] = 
+    lambda[f] =
       *(*basis)[face1[f]] +
       *(*basis)[face2[f]] +
       *(*basis)[face3[f]] +
@@ -187,18 +187,18 @@ HexNodeBasis::HexNodeBasis(int order){
   // Face Based //
   for(int l1 = 1; l1 < order; l1++){
     for(int l2 = 1; l2 < order; l2++){
-      for(int f = 0; f < 6; f++){	
+      for(int f = 0; f < 6; f++){
 	(*basis)[i] = new Polynomial(
 	  legendre[l1].compose(xi[f])  *
 	  legendre[l2].compose(eta[f]) *
 	  lambda[f]);
-	  
+
 	i++;
       }
     }
   }
- 
-  
+
+
   // Cell Based //
   Polynomial px = Polynomial(2, 1, 0, 0);
   Polynomial py = Polynomial(2, 0, 1, 0);
@@ -211,11 +211,11 @@ HexNodeBasis::HexNodeBasis(int order){
   for(int l1 = 1; l1 < order; l1++){
     for(int l2 = 1; l2 < order; l2++){
       for(int l3 = 1; l3 < order; l3++){
-	(*basis)[i] = 
-	  new Polynomial(legendre[l1].compose(px) * 
+	(*basis)[i] =
+	  new Polynomial(legendre[l1].compose(px) *
 			 legendre[l2].compose(py) *
 			 legendre[l3].compose(pz));
-	
+
 	i++;
       }
     }
@@ -237,7 +237,7 @@ HexNodeBasis::~HexNodeBasis(void){
   // Vertex Based //
   for(int i = 0; i < nVertex; i++)
     delete (*node)[i];
-  
+
   delete node;
 
 
@@ -245,10 +245,10 @@ HexNodeBasis::~HexNodeBasis(void){
   for(int c = 0; c < 2; c++){
     for(int i = 0; i < nEdge; i++)
       delete (*(*edge)[c])[i];
-    
+
     delete (*edge)[c];
   }
-  
+
   delete edge;
 
 
@@ -256,7 +256,7 @@ HexNodeBasis::~HexNodeBasis(void){
   for(int c = 0; c < 6; c++){
     for(int i = 0; i < nFace; i++)
       delete (*(*face)[c])[i];
-    
+
     delete (*face)[c];
   }
 
diff --git a/FunctionSpace/Legendre.cpp b/FunctionSpace/Legendre.cpp
index e7ad4f9bbb1a8bebf5ed6d70bc0d910491c832bb..a4e7f72620511fa1763bb3cd524548150e692abe 100644
--- a/FunctionSpace/Legendre.cpp
+++ b/FunctionSpace/Legendre.cpp
@@ -1,6 +1,6 @@
 #include "Legendre.h"
 
-Polynomial Legendre::legendre(const int n, const Polynomial& l, 
+Polynomial Legendre::legendre(const int n, const Polynomial& l,
 			      const Polynomial& lMinus){
   const double nPlus = n + 1;
 
@@ -8,16 +8,16 @@ Polynomial Legendre::legendre(const int n, const Polynomial& l,
     l * Polynomial(1, 1, 0, 0) * ((2 * n + 1) / nPlus) - lMinus * (n / nPlus);
 }
 
-Polynomial Legendre::scaled(const int n, const Polynomial& l, 
+Polynomial Legendre::scaled(const int n, const Polynomial& l,
 			    const Polynomial& lMinus){
   const double nPlus = n + 1;
 
   return
-    l      * Polynomial(1, 1, 0, 0) * ((2 * n + 1) / nPlus) - 
+    l      * Polynomial(1, 1, 0, 0) * ((2 * n + 1) / nPlus) -
     lMinus * Polynomial(1, 0, 2, 0) * ((    n    ) / nPlus);
 }
 
-Polynomial Legendre::integrated(const int n, const Polynomial& l, 
+Polynomial Legendre::integrated(const int n, const Polynomial& l,
 				const Polynomial& lMinus){
   const double nPlus = n + 1;
 
@@ -25,12 +25,12 @@ Polynomial Legendre::integrated(const int n, const Polynomial& l,
     l * Polynomial(1, 1, 0, 0) * ((2 * n - 1) / nPlus) - lMinus * ((n - 2) / nPlus);
 }
 
-Polynomial Legendre::intScaled(const int n, const Polynomial& l, 
+Polynomial Legendre::intScaled(const int n, const Polynomial& l,
 			       const Polynomial& lMinus){
   const double nPlus = n + 1;
 
   return
-    l      * Polynomial(1, 1, 0, 0) * ((2 * n - 1) / nPlus) - 
+    l      * Polynomial(1, 1, 0, 0) * ((2 * n - 1) / nPlus) -
     lMinus * Polynomial(1, 0, 2, 0) * ((    n - 2) / nPlus);
 }
 
@@ -41,7 +41,7 @@ void Legendre::legendre(Polynomial* polynomial, const int order){
 
   if(order >= 0)
     polynomial[0] = Polynomial(1, 0, 0, 0);
-  
+
   if(order >= 1)
     polynomial[1] = Polynomial(1, 1, 0, 0);
 
@@ -49,7 +49,7 @@ void Legendre::legendre(Polynomial* polynomial, const int order){
     for(k = 2; k <= order; k++){
       i = k - 1;
       j = k - 2;
-      
+
       polynomial[k] = legendre(i, polynomial[i], polynomial[j]);
     }
   }
@@ -68,7 +68,7 @@ void Legendre::integrated(Polynomial* polynomial, const int order){
     for(k = 2; k < order; k++){
       i = k - 1;
       j = k - 2;
-      
+
       polynomial[k] = integrated(k, polynomial[i], polynomial[j]);
     }
   }
@@ -79,7 +79,7 @@ void Legendre::scaled(Polynomial* polynomial, const int order){
 
   if(order >= 0)
     polynomial[0] = Polynomial(1, 0, 0, 0);
-  
+
   if(order >= 1)
     polynomial[1] = Polynomial(1, 1, 0, 0);
 
@@ -87,7 +87,7 @@ void Legendre::scaled(Polynomial* polynomial, const int order){
     for(k = 2; k <= order; k++){
       i = k - 1;
       j = k - 2;
-      
+
       polynomial[k] = scaled(i, polynomial[i], polynomial[j]);
     }
   }
@@ -106,7 +106,7 @@ void Legendre::intScaled(Polynomial* polynomial, const int order){
     for(k = 2; k < order; k++){
       i = k - 1;
       j = k - 2;
-      
+
       polynomial[k] = intScaled(k, polynomial[i], polynomial[j]);
     }
   }
diff --git a/FunctionSpace/Legendre.h b/FunctionSpace/Legendre.h
index 4120b65adfcdbe50af71b773a5bb354915f21550..25924558f094f6cc74e3e3dc9e218962f6286c86 100644
--- a/FunctionSpace/Legendre.h
+++ b/FunctionSpace/Legendre.h
@@ -6,8 +6,8 @@
 /**
    @class Legendre
    @brief Generators for legendre Polynomial%s
-   
-   This class handles the generation of legendre 
+
+   This class handles the generation of legendre
    Polynomial%s of many types:
    @li Classical legendre (Legendre::legendre)
    @li Integrated legendre (Legendre::integrated)
@@ -29,50 +29,50 @@ class Legendre{
 
  private:
   static Polynomial legendre(const int n,
-			     const Polynomial& l, 
+			     const Polynomial& l,
 			     const Polynomial& lMinus);
-  
+
   static Polynomial integrated(const int n,
-			       const Polynomial& l, 
+			       const Polynomial& l,
 			       const Polynomial& lMinus);
 
   static Polynomial scaled(const int n,
-			   const Polynomial& l, 
+			   const Polynomial& l,
 			   const Polynomial& lMinus);
 
   static Polynomial intScaled(const int n,
-			      const Polynomial& l, 
+			      const Polynomial& l,
 			      const Polynomial& lMinus);
 };
 
 /**
    @fn void Legendre::legendre(Polynomial*, const int)
-   @param polynomial An @em allocated array (of size @c 'order' @c + @c 1) 
+   @param polynomial An @em allocated array (of size @c 'order' @c + @c 1)
    for storing the requested legendre Polynomial%s
    @param order The @em maximal order of the requested Polynomial%s
    @return Stores in @c 'polynomial' all the @em classical legendre Polynomial%s
-   of order [@c 0, @c 'order']  
+   of order [@c 0, @c 'order']
 
    @fn void Legendre::integrated(Polynomial*, const int)
    @param polynomial An @em allocated array (of size @c 'order')
    for storing the requested legendre Polynomial%s
    @param order The @em maximal order of the requested Polynomial%s
    @return Stores in @c 'polynomial' all the @em integrated legendre Polynomial%s
-   of order [@c 1, @c 'order']  
+   of order [@c 1, @c 'order']
 
    @fn void Legendre::scaled(Polynomial*, const int)
-   @param polynomial An @em allocated array (of size @c 'order' @c + @c 1) 
+   @param polynomial An @em allocated array (of size @c 'order' @c + @c 1)
    for storing the requested legendre Polynomial%s
    @param order The @em maximal order of the requested Polynomial%s
    @return Stores in @c 'polynomial' all the @em scaled legendre Polynomial%s
-   of order [@c 0, @c 'order']  
+   of order [@c 0, @c 'order']
 
    @fn void Legendre::intScaled(Polynomial*, const int)
    @param polynomial An @em allocated array (of size @c 'order')
    for storing the requested legendre Polynomial%s
    @param order The @em maximal order of the requested Polynomial%s
-   @return Stores in @c 'polynomial' all the @em scaled @em integrated 
-   legendre Polynomial%s of order [@c 1, @c 'order']  
+   @return Stores in @c 'polynomial' all the @em scaled @em integrated
+   legendre Polynomial%s of order [@c 1, @c 'order']
  */
 
 #endif
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 6e9c8933f0978c9afbcd56de8240ca23bafc8a11..3839f9b0fc6c8da00a4d2457938e452a2c2fc11e 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -14,7 +14,7 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
 
   // Set Basis Type //
   this->order = order;
-  
+
   type = 0;
   dim  = 1;
 
@@ -43,23 +43,23 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
 
   // Vertex Based (Lagrange) //
   for(unsigned int s = 0; s < nRefSpace; s++){
-    basis[s][0] = 
-      new Polynomial(Polynomial(0.5, 0, 0, 0) - 
+    basis[s][0] =
+      new Polynomial(Polynomial(0.5, 0, 0, 0) -
 		     Polynomial(0.5, 1, 0, 0));
-    
-    basis[s][1] = 
-      new Polynomial(Polynomial(0.5, 0, 0, 0) + 
+
+    basis[s][1] =
+      new Polynomial(Polynomial(0.5, 0, 0, 0) +
 		     Polynomial(0.5, 1, 0, 0));
   }
 
   // Edge Based //
   for(unsigned int s = 0; s < nRefSpace; s++){
     unsigned int i = nVertex;
-    
+
     for(unsigned int l = 1; l < order; l++){
-      basis[s][i] = 
+      basis[s][i] =
 	new Polynomial(intLegendre[l].compose(x[(*(*edgeV[s])[0])[0]]));
-      
+
       i++;
     }
   }
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
index 4b2be2283b293487998ca9c4064305463ae0ccff..d4d41069dd02ad999328ea232a3f37a2400363b6 100644
--- a/FunctionSpace/LineReferenceSpace.cpp
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -19,7 +19,7 @@ LineReferenceSpace::LineReferenceSpace(void){
   // Face Definition //
   nFace   = 0;
   refFace = NULL;
-  
+
   // Init All //
   init();
 }
@@ -28,7 +28,7 @@ LineReferenceSpace::~LineReferenceSpace(void){
   // Delete Ref Edge //
   for(unsigned int i = 0; i < nEdge; i++)
     delete[] refEdge[i];
-  
+
   delete[] refEdge;
 }
 
@@ -36,7 +36,7 @@ string LineReferenceSpace::toLatex(void) const{
   stringstream stream;
 
   stream << "\\documentclass{article}" << endl << endl
-	 
+
 	 << "\\usepackage{longtable}"  << endl
 	 << "\\usepackage{tikz}"       << endl
 	 << "\\usetikzlibrary{arrows}" << endl << endl
@@ -52,10 +52,10 @@ string LineReferenceSpace::toLatex(void) const{
 
 	   << "\\node[vertex] (n0) at(0, 0) {$0$};" << endl
 	   << "\\node[vertex] (n1) at(3, 0) {$1$};" << endl << endl
-      
-	   << "\\path[line]" 
+
+	   << "\\path[line]"
 	   << " (n" << (*(*(*edge)[p])[0])[0] << ")"
-	   << " -- " 
+	   << " -- "
 	   << " (n" << (*(*(*edge)[p])[0])[1] << ");"
 	   << endl
 
@@ -64,6 +64,6 @@ string LineReferenceSpace::toLatex(void) const{
 
   stream << "\\end{longtable}" << endl
 	 << "\\end{document}"  << endl;
-  
+
   return stream.str();
 }
diff --git a/FunctionSpace/LineReferenceSpace.h b/FunctionSpace/LineReferenceSpace.h
index 9059d23ff07162dcb494e96168bf505a76bf3062..ab613949acc073fd1c0e6dc5a170f506e2c15df3 100644
--- a/FunctionSpace/LineReferenceSpace.h
+++ b/FunctionSpace/LineReferenceSpace.h
@@ -16,7 +16,7 @@ class LineReferenceSpace: public ReferenceSpace{
   LineReferenceSpace(void);
   virtual ~LineReferenceSpace(void);
 
-  virtual std::string toLatex(void) const;  
+  virtual std::string toLatex(void) const;
 };
 
 /**
diff --git a/FunctionSpace/Polynomial.cpp b/FunctionSpace/Polynomial.cpp
index fe952cc856924a8b60542b38924673aa2378b4b0..c5a3ca8ee24082007b80c26ecca10e752db41547 100644
--- a/FunctionSpace/Polynomial.cpp
+++ b/FunctionSpace/Polynomial.cpp
@@ -13,16 +13,16 @@ Polynomial::Polynomial(const double coef, const int powerX,
                                           const int powerZ){
   nMon = 1;
   mon  = new monomial_t[1];
-  
+
   mon[0].coef     = coef;
   mon[0].power[0] = powerX;
   mon[0].power[1] = powerY;
   mon[0].power[2] = powerZ;
 }
- 
+
 Polynomial::Polynomial(const Polynomial& other){
   nMon = other.nMon;
-  mon  = copyMonomial(other.mon, nMon); 
+  mon  = copyMonomial(other.mon, nMon);
 }
 
 Polynomial::Polynomial(void){
@@ -35,7 +35,7 @@ Polynomial::~Polynomial(void){
     delete[] mon;
 }
 
-void Polynomial::derivative(const int dim){ 
+void Polynomial::derivative(const int dim){
   // Take derivative //
   for(int i = 0; i < nMon; i++){
     mon[i].coef *= mon[i].power[dim];
@@ -56,19 +56,19 @@ void Polynomial::derivative(const int dim){
   // If no monomial any more ---> return zero polynomial
   if(!N){
     delete[] mon;
-    
+
     mon  = zeroPolynomial();
     nMon = 1;
     return;
   }
-  
+
   // If no zero found ---> return;
-  if(N == nMon) 
+  if(N == nMon)
     return;
 
   // Else, remove them //
   monomial_t* tmp = mon;
-  
+
   mon  = new monomial_t[N];
   nMon = N;
 
@@ -79,7 +79,7 @@ void Polynomial::derivative(const int dim){
 
   delete[] tmp;
 
-  // Sort resulting monomial and return // 
+  // Sort resulting monomial and return //
   sort(mon, nMon);
 
   return;
@@ -97,7 +97,7 @@ vector<Polynomial> Polynomial::gradient(void) const{
   grad[0].derivative(0);
   grad[1].derivative(1);
   grad[2].derivative(2);
-  
+
   return grad;
 }
 
@@ -144,11 +144,11 @@ Polynomial Polynomial::divergence(const vector<Polynomial>& p){
 
 double Polynomial::at
   (const double x, const double y, const double z) const{
-  
+
   double val = 0;
   for(int i = 0; i < nMon; i++){
-    val += mon[i].coef * pow(x, mon[i].power[0]) 
-                       * pow(y, mon[i].power[1]) 
+    val += mon[i].coef * pow(x, mon[i].power[0])
+                       * pow(y, mon[i].power[1])
                        * pow(z, mon[i].power[2]);
   }
 
@@ -160,7 +160,7 @@ fullVector<double> Polynomial::at(const vector<Polynomial>& P,
 				  const double y,
 				  const double z){
   fullVector<double> val(3);
-  
+
   val(0) = P[0].at(x, y, z);
   val(1) = P[1].at(x, y, z);
   val(2) = P[2].at(x, y, z);
@@ -171,7 +171,7 @@ fullVector<double> Polynomial::at(const vector<Polynomial>& P,
 
 Polynomial Polynomial::operator+(const Polynomial& other) const{
   Polynomial newP;
-  
+
   newP.nMon = mergeMon(mon, nMon, other.mon, other.nMon, &newP.mon);
 
   return newP;
@@ -179,11 +179,11 @@ Polynomial Polynomial::operator+(const Polynomial& other) const{
 
 Polynomial Polynomial::operator-(const Polynomial& other) const{
   const int   otherNMon  = other.nMon;
-  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);   
-  
+  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);
+
   Polynomial  newP;
 
-  mult(otherMinus, otherNMon, -1); 
+  mult(otherMinus, otherNMon, -1);
   newP.nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &newP.mon);
 
   delete[] otherMinus;
@@ -192,7 +192,7 @@ Polynomial Polynomial::operator-(const Polynomial& other) const{
 
 Polynomial Polynomial::operator*(const Polynomial& other) const{
   Polynomial newP;
-  
+
   newP.nMon = mult(mon, nMon, other.mon, other.nMon, &newP.mon);
 
   return newP;
@@ -200,7 +200,7 @@ Polynomial Polynomial::operator*(const Polynomial& other) const{
 
 Polynomial Polynomial::operator*(const double alpha) const{
   Polynomial newP;
-  
+
   newP.mon  = copyMonomial(mon, nMon);
   newP.nMon = nMon;
 
@@ -212,7 +212,7 @@ Polynomial Polynomial::operator*(const double alpha) const{
 
 void Polynomial::add(const Polynomial& other){
   monomial_t* tmp = mon;
-  
+
   nMon = mergeMon(mon, nMon, other.mon, other.nMon, &mon);
 
   delete[] tmp;
@@ -220,10 +220,10 @@ void Polynomial::add(const Polynomial& other){
 
 void Polynomial::sub(const Polynomial& other){
   const int   otherNMon  = other.nMon;
-  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);   
+  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);
   monomial_t* tmp        = mon;
-  
-  mult(otherMinus, otherNMon, -1); 
+
+  mult(otherMinus, otherNMon, -1);
   nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &mon);
 
   delete[] otherMinus;
@@ -232,7 +232,7 @@ void Polynomial::sub(const Polynomial& other){
 
 void Polynomial::mul(const Polynomial& other){
   monomial_t* tmp = mon;
-  
+
   nMon = mult(mon, nMon, other.mon, other.nMon, &mon);
 
   delete[] tmp;
@@ -250,21 +250,21 @@ void Polynomial::power(const int n){
   switch(n){
   case 0:
     delete[] mon;
-    
+
     mon  = unitPolynomial();
     nMon = 1;
-    
+
     break;
-    
+
   case 1:
     break;
 
   default:
     Polynomial old = *this;
-    
+
     for(int i = 1; i < n; i++)
       mul(old);
-    
+
     break;
   }
 }
@@ -291,7 +291,7 @@ Polynomial Polynomial::compose(const Polynomial& otherA,
 void Polynomial::operator=(const Polynomial& other){
   if(mon)
     delete[] mon;
-  
+
   nMon = other.nMon;
   mon  = copyMonomial(other.mon, nMon);
 }
@@ -308,7 +308,7 @@ string Polynomial::toString(const Polynomial::monomial_t* mon, const bool isAbs)
   }
 
   // If we're here, we do not have a constant term
-  
+
   // If we have a coefficient of '1', we don't display it
   if(notUnitCoef && isAbs)
     stream << abs(mon->coef);
@@ -328,7 +328,7 @@ string Polynomial::toString(const Polynomial::monomial_t* mon, const bool isAbs)
 	stream << " * ";
 
       stream << coefName[i];
-      
+
       if(mon->power[i] != 1)
 	stream << "^" << mon->power[i];
 
@@ -344,7 +344,7 @@ string Polynomial::toString(void) const{
   bool isAbs = false;
 
   stream << toString(&mon[0], isAbs);
-    
+
   for(int i = 1; i < nMon; i++){
     if(mon[i].coef < 0.0){
       stream << " - ";
@@ -362,17 +362,17 @@ string Polynomial::toString(void) const{
   return stream.str();
 }
 
-bool Polynomial::isSmaller(const Polynomial::monomial_t* a, 
+bool Polynomial::isSmaller(const Polynomial::monomial_t* a,
 			   const Polynomial::monomial_t* b){
-  // GRevLex order: 
+  // GRevLex order:
   // http://www.math.uiuc.edu/Macaulay2/doc/Macaulay2-1.4/share/doc/Macaulay2/Macaulay2Doc/html/___G__Rev__Lex.html
-  
+
   int dif[3];
   int last = 0;
 
   if(isSmallerPower(a, b))
     return true;
-  
+
   if(isEqualPower(a, b)){
     for(int i = 0, j = 2; i < 3; i++, j--)
       dif[i] = b->power[j] - a->power[j];
@@ -383,11 +383,11 @@ bool Polynomial::isSmaller(const Polynomial::monomial_t* a,
 
     if(last < 0)
       return true;
-    
+
     else
       return false;
   }
-    
+
   return false;
 }
 
@@ -405,7 +405,7 @@ void Polynomial::swap(monomial_t* mon, const int i, const int j){
 }
 
 
-int Polynomial::mergeMon(monomial_t* sourceA, const int sizeA, 
+int Polynomial::mergeMon(monomial_t* sourceA, const int sizeA,
 			 monomial_t* sourceB, const int sizeB,
 			 monomial_t** dest){
   stack<monomial_t> s;
@@ -421,20 +421,20 @@ int Polynomial::mergeMon(monomial_t* sourceA, const int sizeA,
 
     else if(sourceB[j].coef == 0.0)
       j++;
-    
+
     else if(isEqual(&sourceA[i], &sourceB[j])){
       tmp       = sourceA[i];
-      tmp.coef += sourceB[j].coef; 
-      
+      tmp.coef += sourceB[j].coef;
+
       if(tmp.coef != 0.0){
 	s.push(tmp);
 	N++;
       }
-      
+
       i++;
       j++;
     }
-    
+
     else if(isSmaller(&sourceA[i], &sourceB[j])){
       s.push(sourceA[i]);
       i++;
@@ -488,11 +488,11 @@ int Polynomial::mult(const monomial_t* sourceA, const int sizeA,
 
   if(sizeA < sizeB){
     a     = sourceA;
-    b     = sourceB; 
+    b     = sourceB;
     nDist = sizeA;
     size  = sizeB;
   }
-  
+
   else{
     a     = sourceB;
     b     = sourceA;
@@ -511,7 +511,7 @@ int Polynomial::mult(const monomial_t* sourceA, const int sizeA,
 
   for(int i = 0; i < nDist; i++){
     dist[i] = copyMonomial(b, size);
-  
+
     distribute(dist[i], size, &a[i]);
   }
 
@@ -522,8 +522,8 @@ int Polynomial::mult(const monomial_t* sourceA, const int sizeA,
 
   for(int i = 1, j = 0; i < nDist; i++, j++){
     tmp[j] = dist[0];
-    
-    finalSize = mergeMon(dist[0], finalSize, 
+
+    finalSize = mergeMon(dist[0], finalSize,
 			 dist[i], size,
 			 &dist[0]);
   }
@@ -536,7 +536,7 @@ int Polynomial::mult(const monomial_t* sourceA, const int sizeA,
     delete[] dist[i];
     delete[] tmp[j];
   }
-  
+
   delete[] dist;
   delete[] tmp;
 
@@ -549,31 +549,31 @@ void Polynomial::mult(monomial_t* source, const int size, const double alpha){
 }
 
 
-void Polynomial::distribute(monomial_t* src, const int size, const monomial_t* m){ 
+void Polynomial::distribute(monomial_t* src, const int size, const monomial_t* m){
   for(int i = 0; i < size; i++){
     src[i].coef *= m->coef;
-    
+
     src[i].power[0] += m->power[0];
     src[i].power[1] += m->power[1];
-    src[i].power[2] += m->power[2];    
+    src[i].power[2] += m->power[2];
   }
 }
 
 void Polynomial::compose(const Polynomial::monomial_t* src,
 			 Polynomial comp,
 			 stack<Polynomial::monomial_t>* stk){
-  
+
   comp.power(src->power[0]);
   comp.mul(src->coef);
-  
+
   const int size = comp.nMon;
 
   for(int i = 0; i < size; i++){
     if(comp.mon[i].coef != 0){
-      
+
       comp.mon[i].power[1] += src->power[1];
       comp.mon[i].power[2] += src->power[2];
-      
+
       stk->push(comp.mon[i]);
     }
   }
@@ -585,17 +585,17 @@ void Polynomial::compose(const Polynomial::monomial_t* src,
 
   compA.power(src->power[0]);
   compB.power(src->power[1]);
-  
+
   compA.mul(compB);
   compA.mul(src->coef);
-  
+
   const int size = compA.nMon;
 
   for(int i = 0; i < size; i++){
     if(compA.mon[i].coef != 0){
-      
+
       compA.mon[i].power[2] += src->power[2];
-      
+
       stk->push(compA.mon[i]);
     }
   }
@@ -617,12 +617,12 @@ Polynomial Polynomial::polynomialFromStack(std::stack<Polynomial::monomial_t>& s
     newMon[0] = stk.top();
     newNMon   = 1;
     stk.pop();
-    
+
     while(!stk.empty()){
       tmp     = newMon;
       newNMon = mergeMon(newMon, newNMon, &stk.top(), 1, &newMon);
       stk.pop();
-      
+
       delete[] tmp;
     }
   }
@@ -648,23 +648,23 @@ Polynomial::monomial_t* Polynomial::copyMonomial(const monomial_t* src, const in
 
 Polynomial::monomial_t* Polynomial::zeroPolynomial(void){
   monomial_t* zero = new monomial_t[1];
-  
+
   zero->coef     = 0;
   zero->power[0] = 0;
   zero->power[1] = 0;
   zero->power[2] = 0;
- 
+
   return zero;
 }
 
 Polynomial::monomial_t* Polynomial::unitPolynomial(void){
   monomial_t* unit = new monomial_t[1];
-  
+
   unit->coef     = 1;
   unit->power[0] = 0;
   unit->power[1] = 0;
   unit->power[2] = 0;
- 
+
   return unit;
 }
 
@@ -705,7 +705,7 @@ int main(void){
   Polynomial p4 = p3 * p3;
   cout << "p3 = " << p3.toString() << endl;
   cout << "p4 = " << p4.toString() << endl;
-  
+
   p3.mul(p3);
   cout << "p3 = " << p3.toString() << endl;
 
@@ -715,10 +715,10 @@ int main(void){
   cout << "p5 = " << p5.toString() << endl;
   cout << "p6 = " << p6.toString() << endl;
   cout << "p7 = " << p7.toString() << endl;
-  
+
   p6.sub(p6);
   cout << "p6 = " << p6.toString() << endl;
-  
+
 
   Polynomial p8(p4);
   Polynomial p9(p4);
diff --git a/FunctionSpace/Polynomial.h b/FunctionSpace/Polynomial.h
index 322f665fcd4d6b65b27b61cbd31bcf615db43037..ba52b8e75eef2cafb9ba259d7665884d151d69f2 100644
--- a/FunctionSpace/Polynomial.h
+++ b/FunctionSpace/Polynomial.h
@@ -16,14 +16,14 @@
 
 // We suppose 3D Polynomial
 class Polynomial{
- private: 
+ private:
   static const char coefName[3];
 
   struct monomial_t{
     double coef;
     int power[3];
   };
-  
+
   int         nMon;
   monomial_t*  mon;
 
@@ -42,10 +42,10 @@ class Polynomial{
   static Polynomial divergence(const std::vector<Polynomial>& p);
 
   double operator()
-    (const double x, const double y, const double z) const;  
-  
+    (const double x, const double y, const double z) const;
+
   double at
-    (const double x, const double y, const double z) const;  
+    (const double x, const double y, const double z) const;
 
   static fullVector<double> at(const std::vector<Polynomial>& P,
 			       const double x,
@@ -65,7 +65,7 @@ class Polynomial{
   void power(const int n);
 
   Polynomial compose(const Polynomial& other) const;
-  Polynomial compose(const Polynomial& otherA, 
+  Polynomial compose(const Polynomial& otherA,
 		     const Polynomial& otherB) const;
 
   void operator=(const Polynomial& other);
@@ -75,31 +75,31 @@ class Polynomial{
  private:
   static std::string toString(const monomial_t* mon, const bool isAbs);
 
-  static bool isSmaller(const monomial_t* a, const monomial_t*b); 
-  static bool isEqual(const monomial_t* a, const monomial_t*b); 
-  static bool isSmallerPower(const monomial_t* a, const monomial_t* b);  
+  static bool isSmaller(const monomial_t* a, const monomial_t*b);
+  static bool isEqual(const monomial_t* a, const monomial_t*b);
+  static bool isSmallerPower(const monomial_t* a, const monomial_t* b);
   static bool isEqualPower(const monomial_t* a, const monomial_t* b);
-  
+
   static void sort(monomial_t* mon, const int size);
   static void swap(monomial_t* mon, const int i, const int j);
 
-  static int mergeMon(monomial_t* sourceA, const int sizeA, 
+  static int mergeMon(monomial_t* sourceA, const int sizeA,
 		      monomial_t* sourceB, const int sizeB,
 		      monomial_t** dest);
-  
+
   static int mult(const monomial_t* sourceA, const int sizeA,
 		  const monomial_t* sourceB, const int sizeB,
 		  monomial_t** dest);
-  
+
   static void mult(monomial_t* source, const int size, const double alpha);
 
   static void distribute(monomial_t* src, const int size, const monomial_t* m);
 
-  static void compose(const monomial_t* src, 
+  static void compose(const monomial_t* src,
 		      Polynomial comp,
 		      std::stack<monomial_t>* stk);
 
-  static void compose(const monomial_t* src, 
+  static void compose(const monomial_t* src,
 		      Polynomial compA, Polynomial compB,
 		      std::stack<monomial_t>* stk);
 
@@ -114,11 +114,11 @@ class Polynomial{
 /**
    @fn Polynomial::Polynomial(const double, const int, const int, const int)
    @param coef The coeficient of the futur monomial
-   @param powerX The power of the '@c x' coordinate 
+   @param powerX The power of the '@c x' coordinate
    of the futur monomial
-   @param powerY The power of the '@c y' coordinate 
+   @param powerY The power of the '@c y' coordinate
    of the futur monomial
-   @param powerZ The power of the '@c z' coordinate 
+   @param powerZ The power of the '@c z' coordinate
    of the futur monomial
    @return Returns a new Monomial with the given
    parameters
@@ -140,7 +140,7 @@ class Polynomial{
    @em with @em no @em monomials.@n
    In particular, the empty Polynomial is @em not
    the @em zero @em Polynomial.@n
-   Indeed, the @em zero @em Polynomial has one monomial, 
+   Indeed, the @em zero @em Polynomial has one monomial,
    @em @c 0.
    **
 
@@ -151,7 +151,7 @@ class Polynomial{
    @fn Polynomial::derivative
    @param dim The dimention to use for the
    derivation
-   @returns Derivates this Polynomial with 
+   @returns Derivates this Polynomial with
    respect to the given dimention
    @note
    Dimention:
@@ -161,16 +161,16 @@ class Polynomial{
    **
 
    @fn Polynomial::gradient
-   @return Returns a Vector with the gradient 
+   @return Returns a Vector with the gradient
    of this Polynomial
    **
-   
+
    @fn Polynomial::curl
    @param p A vector of Polynomial%s
    @return Returns the @em curl of the given
    vector of Polynomial%s
    **
-   
+
    @fn Polynomial::divergence
    @param p A vector of Polynomial%s
    @return Returns the @em divergence of the given
@@ -222,35 +222,35 @@ class Polynomial{
 
    @fn Polynomial Polynomial::operator*(const double) const
    @param alpha A value
-   @return Returns a @em new Polynomial, 
+   @return Returns a @em new Polynomial,
    which is this Polynomial @em multiplied by @c alpha
    **
 
    @fn Polynomial::add
    @param other An other Polynomial
-   @return The given Polynomial is 
+   @return The given Polynomial is
    @em added to this Polynomial
    @note
-   The result of this operation is stored in 
-   this Polynomial  
+   The result of this operation is stored in
+   this Polynomial
    **
 
    @fn Polynomial::sub
    @param other An other Polynomial
-   @return The given Polynomial is 
+   @return The given Polynomial is
    @em substracted to this Polynomial
    @note
-   The result of this operation is stored in 
-   this Polynomial  
+   The result of this operation is stored in
+   this Polynomial
    **
 
    @fn void Polynomial::mul(const Polynomial&)
    @param other An other Polynomial
-   @return The given Polynomial is 
+   @return The given Polynomial is
    @em multiplied with this Polynomial
    @note
-   The result of this operation is stored in 
-   this Polynomial  
+   The result of this operation is stored in
+   this Polynomial
    **
 
    @fn void Polynomial::mul(const double)
@@ -258,8 +258,8 @@ class Polynomial{
    @return This Polynomial is @em multiplied
    by the given value
    @note
-   The result of this operation is stored in 
-   this Polynomial  
+   The result of this operation is stored in
+   this Polynomial
    **
 
    @fn Polynomial::power
@@ -268,7 +268,7 @@ class Polynomial{
    **
 
    @fn Polynomial Polynomial::compose(const Polynomial&) const
-   @param other An other Polynomial, 
+   @param other An other Polynomial,
    called @c Q(x, @c y, @c z)
    @return
    Let this Polynomial be @c P(x, @c y, @c z).@n
@@ -277,13 +277,13 @@ class Polynomial{
    **
 
    @fn Polynomial Polynomial::compose(const Polynomial&, const Polynomial&) const
-   @param otherA An other Polynomial, 
+   @param otherA An other Polynomial,
    called @c Q(x, @c y, @c z)
-   @param otherB An other Polynomial, 
+   @param otherB An other Polynomial,
    called @c R(x, @c y, @c z)
    @return
    Let this Polynomial be @c P(x, @c y, @c z).@n
-   This method returns a @em new Polynomial, representing 
+   This method returns a @em new Polynomial, representing
    @c P(Q(x, @c y, @c z), @c R(x, @c y, @c z), @c z)
    **
 
@@ -301,8 +301,8 @@ class Polynomial{
 // Inline Functions //
 //////////////////////
 
-inline double Polynomial::operator() (const double x, 
-				      const double y, 
+inline double Polynomial::operator() (const double x,
+				      const double y,
 				      const double z) const{
   return at(x, y, z);
 }
@@ -316,20 +316,20 @@ inline bool Polynomial::isEqual(const Polynomial::monomial_t* a,
 
 inline bool Polynomial::isSmallerPower(const Polynomial::monomial_t* a,
 				       const Polynomial::monomial_t* b){
-  
-  return 
-    a->power[0] + a->power[1] + a->power[2] 
-    < 
-    b->power[0] + b->power[1] + b->power[2] ;    
+
+  return
+    a->power[0] + a->power[1] + a->power[2]
+    <
+    b->power[0] + b->power[1] + b->power[2] ;
 }
 
 inline bool Polynomial::isEqualPower(const Polynomial::monomial_t* a,
 				     const Polynomial::monomial_t* b){
-  
-  return 
-    a->power[0] + a->power[1] + a->power[2] 
+
+  return
+    a->power[0] + a->power[1] + a->power[2]
     ==
-    b->power[0] + b->power[1] + b->power[2] ;    
+    b->power[0] + b->power[1] + b->power[2] ;
 }
 
 #endif
diff --git a/FunctionSpace/ReferenceSpace.cpp b/FunctionSpace/ReferenceSpace.cpp
index a56ab022bd65b9dc8b33ad9c8f448e6006763a6a..fca8c8cfac571f54361c7b4e38e144aba22a0af8 100644
--- a/FunctionSpace/ReferenceSpace.cpp
+++ b/FunctionSpace/ReferenceSpace.cpp
@@ -16,25 +16,25 @@ ReferenceSpace::~ReferenceSpace(void){
   // Destroy Tree //
   destroy(&pTreeRoot);
   delete[] perm;
- 
+
   // Delete Permutated Edge //
   for(unsigned int p = 0; p < nPerm; p++){
     for(unsigned int i = 0; i < nEdge; i++)
-      delete (*(*edge)[p])[i];      
-    
+      delete (*(*edge)[p])[i];
+
     delete (*edge)[p];
   }
-  
+
   delete edge;
 
   // Delete Permutated Face //
   for(unsigned int p = 0; p < nPerm; p++){
     for(unsigned int i = 0; i < nFace; i++)
-      delete (*(*face)[p])[i];      
-    
+      delete (*(*face)[p])[i];
+
     delete (*face)[p];
   }
-  
+
   delete face;
 }
 
@@ -44,11 +44,11 @@ void ReferenceSpace::init(void){
   nextLeafId = 0;
 
   pTreeRoot.depth    = 0;
-  pTreeRoot.last     = NULL; 
+  pTreeRoot.last     = NULL;
   pTreeRoot.number   = nVertex;
   pTreeRoot.possible = new unsigned int[pTreeRoot.number];
   pTreeRoot.next     = NULL;
-  
+
   for(unsigned int i = 0; i < pTreeRoot.number; i++)
     pTreeRoot.possible[i] = i;
 
@@ -64,7 +64,7 @@ void ReferenceSpace::init(void){
     perm[i] = lPerm->front();
     lPerm->pop_front();
   }
-  
+
   delete lPerm;
 
   // Get Edges & Faces //
@@ -88,12 +88,12 @@ void ReferenceSpace::populate(node* pTreeRoot){
     // Init Permutation
     pTreeRoot->next   = NULL;
     pTreeRoot->leafId = nextLeafId;
-    
+
     // Value for Next Permutation
     nextLeafId++;
     nPerm++;
-    
-    // Put this Permutation in queue 
+
+    // Put this Permutation in queue
     //                (AND IN ORDER)
     lPerm->push_back(pTreeRoot->last);
   }
@@ -103,27 +103,27 @@ void ReferenceSpace::populate(node* pTreeRoot){
     // We got 'number' child nodes
     pTreeRoot->next = new node[number];
 
-    // Init each child node 
+    // Init each child node
     for(unsigned int i = 0; i < number; i++){
       nextLast = pTreeRoot->possible[i];
       offset   = 0;
-      
+
       // Depth and Last Choices of child nodes
       pTreeRoot->next[i].depth       = nextDepth;
       pTreeRoot->next[i].last        = new unsigned int[nextDepth];
       pTreeRoot->next[i].last[depth] = nextLast;
-      
+
       for(unsigned int j = 0; j < depth; j++)
 	pTreeRoot->next[i].last[j] = pTreeRoot->last[j];
 
       // Possibilities of child node
       pTreeRoot->next[i].number   = nextNumber;
       pTreeRoot->next[i].possible = new unsigned int[nextNumber];
-      
+
       for(unsigned int j = 0; j < nextNumber; j++){
 	  if(pTreeRoot->possible[j] == nextLast)
 	    offset = 1;
-	  
+
 	  pTreeRoot->next[i].possible[j] = pTreeRoot->possible[j + offset];
       }
 
@@ -140,7 +140,7 @@ void ReferenceSpace::destroy(node* node){
     destroy(&node->next[i]);
     node->number--;
   }
-    
+
   delete[] node->possible;
   delete[] node->last;
   delete[] node->next;
@@ -154,9 +154,9 @@ void ReferenceSpace::getEdge(void){
     tmp = new vector<const vector<unsigned int>*>(nEdge);
 
     for(unsigned int i = 0; i < nEdge; i++)
-      (*tmp)[i] = inOrder(p, 
-			  refEdge[i][0], 
-			  refEdge[i][1]);   
+      (*tmp)[i] = inOrder(p,
+			  refEdge[i][0],
+			  refEdge[i][1]);
     (*edge)[p] = tmp;
   }
 }
@@ -169,57 +169,57 @@ void ReferenceSpace::getFace(void){
     tmp = new vector<const vector<unsigned int>*>(nFace);
 
     for(unsigned int i = 0; i < nFace; i++)
-      (*tmp)[i] = inOrder(p, 
-			  refFace[i][0], 
+      (*tmp)[i] = inOrder(p,
+			  refFace[i][0],
 			  refFace[i][1],
-			  refFace[i][2]);   
+			  refFace[i][2]);
     (*face)[p] = tmp;
   }
 }
 
 const vector<unsigned int>* ReferenceSpace::
-inOrder(unsigned int permutation, 
-	unsigned int a, 
+inOrder(unsigned int permutation,
+	unsigned int a,
 	unsigned int b){
-  
+
   unsigned int v;
   unsigned int k = 0;
-  vector<unsigned int>* inorder = 
+  vector<unsigned int>* inorder =
     new vector<unsigned int>(2);
 
   for(unsigned int i = 0; i < nVertex; i++){
     v = perm[permutation][i];
-    
-    if(v == a || v == b){   
+
+    if(v == a || v == b){
       (*inorder)[k] = v;
       k++;
     }
   }
-  
+
   return inorder;
 }
 
 const std::vector<unsigned int>* ReferenceSpace::
-inOrder(unsigned int permutation, 
-	unsigned int a, 
+inOrder(unsigned int permutation,
+	unsigned int a,
 	unsigned int b,
 	unsigned int c){
-  
+
   unsigned int v;
   unsigned int k = 0;
-  vector<unsigned int>* inorder = 
+  vector<unsigned int>* inorder =
     new vector<unsigned int>(3);
 
   for(unsigned int i = 0; i < nVertex; i++){
     v = perm[permutation][i];
-    
-    if(v == a || v == b || v == c){   
+
+    if(v == a || v == b || v == c){
       (*inorder)[k] = v;
       k++;
     }
   }
-  
-  return inorder;  
+
+  return inorder;
 }
 
 unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
@@ -232,12 +232,12 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
 
   for(int i = 0; i < nVertex; i++){
     vertex[i].first  = i;
-    vertex[i].second = element.getVertex(i);  
+    vertex[i].second = element.getVertex(i);
   }
-  
-  // Sort Them with repsect to Vertex Global ID // 
+
+  // Sort Them with repsect to Vertex Global ID //
   //                 (vertex[i].second->getNum) //
-  std::sort(vertex.begin(), vertex.end(), sortPredicate);  
+  std::sort(vertex.begin(), vertex.end(), sortPredicate);
 
   // Tree Lookup //
   try{
@@ -251,7 +251,7 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
 }
 
 unsigned int ReferenceSpace::treeLookup(const node* root,
-					vector<pair<unsigned int, MVertex*> >& 
+					vector<pair<unsigned int, MVertex*> >&
 					sortedArray){
   // Temp Data //
   unsigned int choice;
@@ -277,22 +277,22 @@ unsigned int ReferenceSpace::treeLookup(const node* root,
 
   // Else: Return Leaf ID //
   else
-    return root->leafId; 
+    return root->leafId;
 }
 
 string ReferenceSpace::toString(void) const{
   stringstream  stream;
-  
+
   // Tree //
   stream << "Tree:"              << endl;
   stream << toString(&pTreeRoot) << endl;
 
   // ReferenceSpaces //
   stream << "Reference Spaces:" << endl;
-  
+
   for(unsigned int i = 0; i < nPerm; i++){
     stream << "  * ";
-    
+
     for(unsigned int j = 0; j < nVertex; j++)
       stream << perm[i][j] << " ";
 
@@ -300,23 +300,23 @@ string ReferenceSpace::toString(void) const{
   }
 
   stream << "Edges Permutations:" << endl;
-  
+
   for(unsigned int i = 0; i < nPerm; i++){
     stream << "  * RefSpace #" << i + 1 << ":" << endl;
-    
+
     for(unsigned int j = 0; j < nEdge; j++)
-      stream << "      -- [" 
+      stream << "      -- ["
 	     << edge->at(i)->at(j)->at(0) << ", "
 	     << edge->at(i)->at(j)->at(1) << "]" << endl;
   }
 
   stream << "Faces Permutations:" << endl;
-  
+
   for(unsigned int i = 0; i < nPerm; i++){
     stream << "  * RefSpace #" << i + 1 << ":" << endl;
-    
+
     for(unsigned int j = 0; j < nFace; j++)
-      stream << "      -- [" 
+      stream << "      -- ["
 	     << face->at(i)->at(j)->at(0) << ", "
 	     << face->at(i)->at(j)->at(1) << ", "
 	     << face->at(i)->at(j)->at(2) << "]" << endl;
@@ -338,6 +338,6 @@ string ReferenceSpace::toString(const node* node) const{
     stream << toString(&node->next[i]);
 
   stream << ")";
-  
+
   return stream.str();
 }
diff --git a/FunctionSpace/ReferenceSpace.h b/FunctionSpace/ReferenceSpace.h
index 580d56239f7c2c474f04134dac109e692eb4eabf..88e86a097839f6ce35fc3ec08bb0aac00bc19cfa 100644
--- a/FunctionSpace/ReferenceSpace.h
+++ b/FunctionSpace/ReferenceSpace.h
@@ -25,7 +25,7 @@ class ReferenceSpace{
     unsigned int  number;   // Number of Next Choise
     unsigned int* possible; // Possible  Next Choise
     node_s*       next;     // Next           Choise
-  
+
     unsigned int  leafId;   // If leaf: this leaf number
                             //     (from 0 to nLeaf - 1)
                             // Else: no meaning
@@ -57,12 +57,12 @@ class ReferenceSpace{
   virtual ~ReferenceSpace(void);
 
   unsigned int getNPermutation(void) const;
-  
+
   unsigned int getPermutation(const MElement& element) const;
 
   std::vector<const std::vector<const std::vector<unsigned int>*>*>&
     getAllEdge(void) const;
-  
+
   std::vector<const std::vector<const std::vector<unsigned int>*>*>&
     getAllFace(void) const ;
 
@@ -76,24 +76,24 @@ class ReferenceSpace{
   void init(void);
   void populate(node* pTreeRoot);
   void destroy(node* node);
-  
+
   void getEdge(void);
   void getFace(void);
 
-  const std::vector<unsigned int>* inOrder(unsigned int permutation, 
-					   unsigned int a, 
+  const std::vector<unsigned int>* inOrder(unsigned int permutation,
+					   unsigned int a,
 					   unsigned int b);
 
-  const std::vector<unsigned int>* inOrder(unsigned int permutation, 
-					   unsigned int a, 
+  const std::vector<unsigned int>* inOrder(unsigned int permutation,
+					   unsigned int a,
 					   unsigned int b,
 					   unsigned int c);
 
-  static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a, 
+  static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a,
 			    const std::pair<unsigned int, MVertex*>& b);
 
   static unsigned int treeLookup(const node* root,
-				 std::vector<std::pair<unsigned int, MVertex*> >& 
+				 std::vector<std::pair<unsigned int, MVertex*> >&
 				 sortedArray);
 
   std::string toString(const node* node) const;
@@ -110,7 +110,7 @@ class ReferenceSpace{
    @fn ReferenceSpace::~ReferenceSpace
    Deletes this ReferenceSpace
    **
-   
+
    @fn ReferenceSpace::getNPermutation
    @returns Returns the number of permutation of this
    ReferenceSpace
@@ -118,35 +118,35 @@ class ReferenceSpace{
 
    @fn ReferenceSpace::getPermutation
    @param element A MElement
-   @returns Returns the a natural number defining 
+   @returns Returns the a natural number defining
    the permutation of the given element
-   
+
    @note If no permutation is found (e.g. the given
-   element does not belong the @em same @em geometrical 
+   element does not belong the @em same @em geometrical
    entity as this ReferenceSpace) an Exception is thrown
    **
 
    @fn ReferenceSpace::getAllEdge
    @return Returns every Edge permutation of this ReferenceSpace
-   
+
    @note
-   @li The fisrt vector represents a particular permutation 
+   @li The fisrt vector represents a particular permutation
    (see ReferenceSpace::getPermutation())
-   @li The second vector represents a particular edge 
-   (for a given permutation) 
-   @li The last vector represents the Vertex @c IDs of 
+   @li The second vector represents a particular edge
+   (for a given permutation)
+   @li The last vector represents the Vertex @c IDs of
    the given edge (in the @em geometrical reference space)
    **
 
    @fn ReferenceSpace::getAllFace
    @return Returns every Face permutation of this ReferenceSpace
-   
+
    @note
-   @li The fisrt vector represents a particular permutation 
+   @li The fisrt vector represents a particular permutation
    (see ReferenceSpace::getPermutation())
-   @li The second vector represents a particular face 
-   (for a given permutation) 
-   @li The last vector represents the Vertex @c IDs of 
+   @li The second vector represents a particular face
+   (for a given permutation)
+   @li The last vector represents the Vertex @c IDs of
    the given face (in the @em geometrical reference space)
    **
 
@@ -155,21 +155,21 @@ class ReferenceSpace{
    **
 
    @fn ReferenceSpace::toLatex
-   @return Returns a string (of a Latex file) 
-   describing this ReferenceSpace   
+   @return Returns a string (of a Latex file)
+   describing this ReferenceSpace
  */
 
 //////////////////////
 // Inline Functions //
 //////////////////////
 
-inline 
-unsigned int 
+inline
+unsigned int
 ReferenceSpace::getNPermutation(void) const{
   return nPerm;
 }
-  
-inline 
+
+inline
 std::vector<const std::vector<const std::vector<unsigned int>*>*>&
 ReferenceSpace::getAllEdge(void) const{
   return *edge;
@@ -181,9 +181,9 @@ ReferenceSpace::getAllFace(void) const{
   return *face;
 }
 
-inline 
-bool 
-ReferenceSpace::sortPredicate(const std::pair<unsigned int, MVertex*>& a, 
+inline
+bool
+ReferenceSpace::sortPredicate(const std::pair<unsigned int, MVertex*>& a,
 			      const std::pair<unsigned int, MVertex*>& b){
   return a.second->getNum() < b.second->getNum();
 }
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 4f36f7a8d21a31d6da279e9c48bbe3eaa86a7a61..37971f323cde740a0898ed833b8bd1674a104789 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -17,7 +17,7 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
 
   // Set Basis Type //
   this->order = order;
-  
+
   type = 0;
   dim  = 3;
 
@@ -42,17 +42,17 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
   Legendre::intScaled(intLegendre, order);
 
   // Lagrange Polynomial //
-  const Polynomial lagrange[4] = 
+  const Polynomial lagrange[4] =
     {
-      Polynomial(Polynomial(1, 0, 0, 0) - 
-		 Polynomial(1, 1, 0, 0) - 
-		 Polynomial(1, 0, 1, 0) - 
+      Polynomial(Polynomial(1, 0, 0, 0) -
+		 Polynomial(1, 1, 0, 0) -
+		 Polynomial(1, 0, 1, 0) -
 		 Polynomial(1, 0, 0, 1)),
- 
+
       Polynomial(Polynomial(1, 1, 0, 0)),
 
       Polynomial(Polynomial(1, 0, 1, 0)),
-      
+
       Polynomial(Polynomial(1, 0, 0, 1))
     };
 
@@ -77,12 +77,12 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
 
     for(int e = 0; e < 6; e++){
       for(unsigned int l = 1; l < order; l++){
-	basis[s][i] = 
+	basis[s][i] =
 	  new Polynomial(intLegendre[l].compose
-			 (lagrange[(*(*edgeV[s])[e])[0]] - 
+			 (lagrange[(*(*edgeV[s])[e])[0]] -
 			  lagrange[(*(*edgeV[s])[e])[1]]
-			  , 
-			  lagrange[(*(*edgeV[s])[e])[0]] + 
+			  ,
+			  lagrange[(*(*edgeV[s])[e])[0]] +
 			  lagrange[(*(*edgeV[s])[e])[1]]));
 	i++;
       }
@@ -96,23 +96,23 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
     for(int f = 0; f < 4; f++){
       for(int l1 = 1; l1 < orderMinus; l1++){
 	for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
-	  Polynomial sum = 
-	    lagrange[(*(*faceV[s])[f])[0]] + 
-	    lagrange[(*(*faceV[s])[f])[1]] + 
-	    lagrange[(*(*faceV[s])[f])[2]];	  
-	  
-	  basis[s][i] = 
+	  Polynomial sum =
+	    lagrange[(*(*faceV[s])[f])[0]] +
+	    lagrange[(*(*faceV[s])[f])[1]] +
+	    lagrange[(*(*faceV[s])[f])[2]];
+
+	  basis[s][i] =
 	    new Polynomial(intLegendre[l1].compose
-			   (lagrange[(*(*faceV[s])[f])[0]] - 
+			   (lagrange[(*(*faceV[s])[f])[0]] -
 			    lagrange[(*(*faceV[s])[f])[1]]
-			    , 
-			    lagrange[(*(*faceV[s])[f])[0]] + 
-			    lagrange[(*(*faceV[s])[f])[1]]) 
-			   
-			   * 
-			   
-			   lagrange[(*(*faceV[s])[f])[2]] 
-			   * 
+			    ,
+			    lagrange[(*(*faceV[s])[f])[0]] +
+			    lagrange[(*(*faceV[s])[f])[1]])
+
+			   *
+
+			   lagrange[(*(*faceV[s])[f])[2]]
+			   *
 			   sclLegendre[l2].compose
 			   (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum));
 	  i++;
@@ -120,7 +120,7 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
       }
     }
   }
-  
+
   // Cell Based //
   const Polynomial oneMinusFour         = Polynomial(1, 0, 0, 0) - lagrange[3];
   const Polynomial twoThreeOneMinusFour = lagrange[2] * 2 - oneMinusFour;
@@ -131,14 +131,14 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
 
   for(unsigned int s = 0; s < nRefSpace; s++){
     unsigned int i = nVertex + nEdge + nFace;
-  
+
     for(int l1 = 1; l1 < orderMinusTwo; l1++){
       for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
 	for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
-	  basis[s][i] = 
-	    new Polynomial(intLegendre[l1].compose(sub, add)             * 
-			   lagrange[2]                                   * 		 
-			   sclLegendre[l2].compose(twoThreeOneMinusFour, 
+	  basis[s][i] =
+	    new Polynomial(intLegendre[l1].compose(sub, add)             *
+			   lagrange[2]                                   *
+			   sclLegendre[l2].compose(twoThreeOneMinusFour,
 						   oneMinusFour)         *
 			   lagrange[3]                                   *
 			   legendre[l3].compose(twoFourMinusOne));
diff --git a/FunctionSpace/TetReferenceSpace.cpp b/FunctionSpace/TetReferenceSpace.cpp
index dc8d039c3a1a87b8f9e59f1dc9e2d27260348a5d..247172bb05bf39fdf8b7d2f61e7aec9cc9f128d7 100644
--- a/FunctionSpace/TetReferenceSpace.cpp
+++ b/FunctionSpace/TetReferenceSpace.cpp
@@ -15,17 +15,17 @@ TetReferenceSpace::TetReferenceSpace(void){
   for(unsigned int i = 0; i < nEdge; i++){
     refEdge[i]    = new unsigned int[2];
     refEdge[i][0] = MTetrahedron::edges_tetra(i, 0);
-    refEdge[i][1] = MTetrahedron::edges_tetra(i, 1); 
+    refEdge[i][1] = MTetrahedron::edges_tetra(i, 1);
   }
-  
+
   // Face Definition //
   nFace   = 4;
   refFace = new unsigned int*[nFace];
 
   for(unsigned int i = 0; i < nFace; i++){
     refFace[i]    = new unsigned int[3];
-    refFace[i][0] = MTetrahedron::faces_tetra(i, 0); 
-    refFace[i][1] = MTetrahedron::faces_tetra(i, 1); 
+    refFace[i][0] = MTetrahedron::faces_tetra(i, 0);
+    refFace[i][1] = MTetrahedron::faces_tetra(i, 1);
     refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
   }
 
@@ -37,13 +37,13 @@ TetReferenceSpace::~TetReferenceSpace(void){
   // Delete Ref Edge //
   for(unsigned int i = 0; i < nEdge; i++)
     delete[] refEdge[i];
-  
+
   delete[] refEdge;
 
   // Delete Ref Face //
   for(unsigned int i = 0; i < nFace; i++)
     delete[] refFace[i];
-  
+
   delete[] refFace;
 }
 
@@ -51,7 +51,7 @@ string TetReferenceSpace::toLatex(void) const{
   stringstream stream;
 
   stream << "\\documentclass{article}" << endl << endl
-	 
+
 	 << "\\usepackage{longtable}"  << endl
 	 << "\\usepackage{tikz}"       << endl
 	 << "\\usetikzlibrary{arrows}" << endl << endl
@@ -67,26 +67,26 @@ string TetReferenceSpace::toLatex(void) const{
 
 	   << "\\node[vertex] (n0) at(0, 0) {$0$};" << endl
 	   << "\\node[vertex] (n1) at(3, 0) {$1$};" << endl
-	   << "\\node[vertex] (n2) at(0, 3) {$2$};" << endl 
+	   << "\\node[vertex] (n2) at(0, 3) {$2$};" << endl
 	   << "\\node[vertex] (n3) at(1, 1) {$3$};" << endl << endl;
 
     for(unsigned int i = 0; i < 6; i++)
-      stream << "\\path[line]" 
+      stream << "\\path[line]"
 	     << " (n" << (*(*(*edge)[p])[i])[0] << ")"
-	     << " -- " 
+	     << " -- "
 	     << " (n" << (*(*(*edge)[p])[i])[1] << ");"
 	     << endl;
 
     if((p + 1) % 3)
       stream << "\\end{tikzpicture} & "        << endl << endl;
-    
+
     else
       stream << "\\end{tikzpicture} \\\\ \\\\" << endl << endl;
   }
 
   stream << "\\end{longtable}" << endl
 	 << "\\end{document}"  << endl;
-  
+
   return stream.str();
 }
 
@@ -96,7 +96,7 @@ int main(void){
   TetReferenceSpace p;
 
   cout << p.toLatex() << endl;
-  
+
   return 0;
 }
 */
diff --git a/FunctionSpace/TetReferenceSpace.h b/FunctionSpace/TetReferenceSpace.h
index 4895fa849091d7c619d7f33cca208c9533894dd0..cec36e26ad534d90a6eec4859492b293ad891da3 100644
--- a/FunctionSpace/TetReferenceSpace.h
+++ b/FunctionSpace/TetReferenceSpace.h
@@ -16,7 +16,7 @@ class TetReferenceSpace: public ReferenceSpace{
   TetReferenceSpace(void);
   virtual ~TetReferenceSpace(void);
 
-  virtual std::string toLatex(void) const;  
+  virtual std::string toLatex(void) const;
 };
 
 /**
diff --git a/FunctionSpace/TriLagrangeBasis.h b/FunctionSpace/TriLagrangeBasis.h
index a008a69fb29911515050d6bb1fbea7671cafd163..272a21f32d84e27dea05dbe852389e507260aa40 100644
--- a/FunctionSpace/TriLagrangeBasis.h
+++ b/FunctionSpace/TriLagrangeBasis.h
@@ -6,11 +6,11 @@
 /**
    @class TriLagrangeBasis
    @brief Lagrange Basis for Triangles
- 
-   This class can instantiate a @em Lagrange @em Basis 
+
+   This class can instantiate a @em Lagrange @em Basis
    for a Triangle and for a given Order.@n
-   
-   It uses 
+
+   It uses
    <a href="http://geuz.org/gmsh/">gmsh</a> Basis.
  */
 
@@ -18,21 +18,21 @@ class TriLagrangeBasis: public BasisLagrange{
  public:
   //! @param order A natural number
   //!
-  //! Returns a new TriLagrangeBasis 
+  //! Returns a new TriLagrangeBasis
   //! of the given Order
   TriLagrangeBasis(unsigned int order);
-  
+
   //! Deletes this Basis
   //!
   virtual ~TriLagrangeBasis(void);
 
- private:  
+ private:
   //! @param order A natural number
-  //! @return Returns the @em tag of a @em Triangle of 
-  //! the given order 
+  //! @return Returns the @em tag of a @em Triangle of
+  //! the given order
   static unsigned int getTag(unsigned int order);
 
-  //! @param order A natural number 
+  //! @param order A natural number
   //! @return Returns Lagrangian Points on a Triangle
   //! for the given Order
   static fullMatrix<double> triPoint(unsigned int order);
diff --git a/FunctionSpace/TriReferenceSpace.cpp b/FunctionSpace/TriReferenceSpace.cpp
index f3b40cd9b598ae136da7d4653f189d75e2f7ab26..e4a1d15850f16d43224d50b872964812253f13b5 100644
--- a/FunctionSpace/TriReferenceSpace.cpp
+++ b/FunctionSpace/TriReferenceSpace.cpp
@@ -15,13 +15,13 @@ TriReferenceSpace::TriReferenceSpace(void){
   for(unsigned int i = 0; i < nEdge; i++){
     refEdge[i]    = new unsigned int[2];
     refEdge[i][0] = MTriangle::edges_tri(i, 0);
-    refEdge[i][1] = MTriangle::edges_tri(i, 1); 
+    refEdge[i][1] = MTriangle::edges_tri(i, 1);
   }
 
   // Face Definition //
   nFace   = 0;
   refFace = NULL;
-  
+
   // Init All //
   init();
 }
@@ -30,7 +30,7 @@ TriReferenceSpace::~TriReferenceSpace(void){
   // Delete Ref Edge //
   for(unsigned int i = 0; i < nEdge; i++)
     delete[] refEdge[i];
-  
+
   delete[] refEdge;
 }
 
@@ -38,7 +38,7 @@ string TriReferenceSpace::toLatex(void) const{
   stringstream stream;
 
   stream << "\\documentclass{article}" << endl << endl
-	 
+
 	 << "\\usepackage{longtable}"  << endl
 	 << "\\usepackage{tikz}"       << endl
 	 << "\\usetikzlibrary{arrows}" << endl << endl
@@ -57,22 +57,22 @@ string TriReferenceSpace::toLatex(void) const{
 	   << "\\node[vertex] (n2) at(0, 3) {$2$};" << endl << endl;
 
     for(unsigned int i = 0; i < 3; i++)
-      stream << "\\path[line]" 
+      stream << "\\path[line]"
 	     << " (n" << (*(*(*edge)[p])[i])[0] << ")"
-	     << " -- " 
+	     << " -- "
 	     << " (n" << (*(*(*edge)[p])[i])[1] << ");"
 	     << endl;
 
     if((p + 1) % 3)
       stream << "\\end{tikzpicture} & "        << endl << endl;
-    
+
     else
       stream << "\\end{tikzpicture} \\\\ \\\\" << endl << endl;
   }
 
   stream << "\\end{longtable}" << endl
 	 << "\\end{document}"  << endl;
-  
+
   return stream.str();
 }
 
@@ -82,7 +82,7 @@ int main(void){
   TriReferenceSpace p;
 
   cout << p.toLatex() << endl;
-  
+
   return 0;
 }
 */
diff --git a/FunctionSpace/TriReferenceSpace.h b/FunctionSpace/TriReferenceSpace.h
index 0ac863a4162c7704c8d10033f13e54b2cdaf2bae..9f52b5d64eb9fe56dd5fff7e840c95d15bc77b25 100644
--- a/FunctionSpace/TriReferenceSpace.h
+++ b/FunctionSpace/TriReferenceSpace.h
@@ -16,7 +16,7 @@ class TriReferenceSpace: public ReferenceSpace{
   TriReferenceSpace(void);
   virtual ~TriReferenceSpace(void);
 
-  virtual std::string toLatex(void) const;  
+  virtual std::string toLatex(void) const;
 };
 
 /**