diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index bc9f4a22cd5d3cf16a3759c5484393410e9f9dfa..6582c54ef22abc2fb47b76490798134fe589500e 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -15,6 +15,7 @@ set(SRC
TriNodeBasis.cpp
TriEdgeBasis.cpp
TriNedelecBasis.cpp
+ HexNodeBasis.cpp
LocalFunctionSpace.cpp
LocalFunctionSpaceScalar.cpp
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b78cf12e4ec0aa45bdb868336e1fd3592b8b241f
--- /dev/null
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -0,0 +1,247 @@
+#include "HexNodeBasis.h"
+#include "Legendre.h"
+
+HexNodeBasis::HexNodeBasis(const int order){
+ // Set Basis Type //
+ this->order = order;
+
+ type = 0;
+ size = (order + 1) * (order + 1) * (order + 1);
+ nodeNbr = 8;
+ dim = 3;
+
+ // Alloc Temporary Space //
+ Polynomial* legendre = new Polynomial[order];
+ Polynomial* lifting = new Polynomial[8];
+
+ // Legendre Polynomial //
+ Legendre::integrated(legendre, order);
+
+ // Lifting //
+ 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] =
+ (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] =
+ (Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+
+ 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] =
+ (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] =
+ (Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) +
+ Polynomial(1, 0, 0, 1);
+
+ lifting[6] =
+ (Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ Polynomial(1, 0, 0, 1);
+
+ lifting[7] =
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+ (Polynomial(1, 0, 1, 0)) +
+ Polynomial(1, 0, 0, 1);
+
+
+ // Basis //
+ basis = new std::vector<Polynomial>(size);
+
+
+ // Vertex Based (Lagrange) //
+ (*basis)[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));
+
+ (*basis)[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));
+
+ (*basis)[2] =
+ (Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 0, 1));
+
+ (*basis)[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));
+
+ (*basis)[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);
+
+ (*basis)[5] =
+ (Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1);
+
+ (*basis)[6] =
+ (Polynomial(1, 1, 0, 0)) *
+ (Polynomial(1, 0, 1, 0)) *
+ Polynomial(1, 0, 0, 1);
+
+ (*basis)[7] =
+ (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;
+
+ // Points definig Edges
+ int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
+ int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
+
+ for(int l = 1; l < order; l++){
+ for(int e = 0; e < 12; e++){
+ (*basis)[i] =
+ legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) *
+ ((*basis)[edge1[e]] + (*basis)[edge2[e]]);
+
+ i++;
+ }
+ }
+
+
+ // Face Based //
+ // Points definig Faces
+ int face1[6] = {0, 3, 2, 1, 5, 4};
+ int face2[6] = {1, 7, 6, 0, 6, 7};
+ int face3[6] = {2, 6, 5, 4, 7, 3};
+ int face4[6] = {3, 2, 1, 5, 4, 0};
+
+ for(int l1 = 1; l1 < order; l1++){
+ for(int l2 = 1; l2 < order; l2++){
+ for(int f = 0; f < 6; f++){
+ (*basis)[i] = Polynomial(42, 0, 0, 0);
+
+ i++;
+ }
+ }
+ }
+
+
+ // Cell Based //
+ Polynomial px = Polynomial(2, 1, 0, 0);
+ Polynomial py = Polynomial(2, 0, 1, 0);
+ Polynomial pz = Polynomial(2, 0, 0, 1);
+
+ px = px - Polynomial(1, 0, 0, 0);
+ py = py - Polynomial(1, 0, 0, 0);
+ pz = pz - Polynomial(1, 0, 0, 0);
+
+ for(int l1 = 1; l1 < order; l1++){
+ for(int l2 = 1; l2 < order; l2++){
+ for(int l3 = 1; l3 < order; l3++){
+ (*basis)[i] = legendre[l1].compose(px) *
+ legendre[l2].compose(py) *
+ legendre[l3].compose(pz);
+
+ i++;
+ }
+ }
+ }
+
+
+ // Free Temporary Sapce //
+ delete[] legendre;
+ delete[] lifting;
+}
+
+HexNodeBasis::~HexNodeBasis(void){
+ delete basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+
+ const int P = 2;
+ const double d = 0.05;
+
+ HexNodeBasis b(P);
+
+ const std::vector<Polynomial>& basis = b.getBasis();
+
+ printf("\n");
+ printf("clear all;\n");
+ printf("\n");
+
+ printf("\n");
+ printf("Order = %d\n", b.getOrder());
+ printf("Type = %d\n", b.getType());
+ printf("Size = %d\n", b.getSize());
+ printf("NodeNumber = %d\n", b.getNodeNbr());
+ printf("Dimension = %d\n", b.getDim());
+ printf("\n");
+
+ printf("function r = p(i, x, y, z)\n");
+ printf("p = zeros(%d, 1);\n", b.getSize());
+ printf("\n");
+
+ for(int i = 0; i < b.getSize(); i++)
+ printf("p(%d) = %s;\n", i + 1, basis[i].toString().c_str());
+
+ printf("\n");
+ printf("r = p(i, 1);\n");
+ printf("end\n");
+ printf("\n");
+ /*
+ printf("d = %f;\nx = [0:d:1];\ny = x;\nz = x;\n\nlx = length(x);\nly = length(y);\nlz = length(z);\n\n", d);
+
+ for(int i = 0; i < b.getSize(); i++)
+ printf("p%d = zeros(lx, ly, lz);\n", i + 1);
+
+ printf("\n");
+ printf("for i = 1:lx\n");
+ printf("for j = 1:ly\n");
+ printf("for k = 1:lz\n");
+ printf("\n");
+
+ for(int i = 0; i < b.getSize(); i++)
+ printf("p%d(j, i, k) = p(%d, x(i), y(j), z(k));\n", i + 1, i + 1);
+
+ printf("end\n");
+ printf("end\n");
+ printf("end\n");
+
+ printf("\n");
+ printf("SizeOfBasis = %lu\n", sizeof(b) + sizeof(basis) * b.getSize());
+ printf("\n");
+
+ printf("\n");
+ for(int i = b.getSize(); i > 0; i--){
+ printf("figure;\n");
+ printf("hold on;\n");
+ printf("for k = 1:lz\n");
+ printf("surf(x, y, p%d(:, :, k) + z(k));\n", i);
+ printf("end\n");
+ printf("hold off;\n");
+ }
+
+ printf("\n");
+ */
+ return 0;
+}
+*/
diff --git a/FunctionSpace/HexNodeBasis.h b/FunctionSpace/HexNodeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..d25a05617cb96a86c196556dc38094d39b69ee84
--- /dev/null
+++ b/FunctionSpace/HexNodeBasis.h
@@ -0,0 +1,29 @@
+#ifndef _HEXNODEBASIS_H_
+#define _HEXNODEBASIS_H_
+
+#include "BasisScalar.h"
+
+/**
+ @class HexNodeBasis
+ @brief A Node-Basis for Hexahedrons
+
+ This class can instantiate a Node-Based Basis
+ (high or low order) for Hexahedrons.@n
+
+ It uses
+ <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
+ Basis for @em high @em order Polynomial%s generation.@n
+ */
+
+class HexNodeBasis: public BasisScalar{
+ public:
+ //! Returns a new Node-Basis for Hexahedrons of the given order
+ //! @param order The order of the Basis
+ HexNodeBasis(const int order);
+
+ //! @return Deletes this Basis
+ //!
+ virtual ~HexNodeBasis(void);
+};
+
+#endif