Edge and Node Basis (Tri and Quad)

FunctionSpace/Basis.cpp

+#include "Basis.h"
+
+Basis::Basis(void){
+}
+
+Basis::~Basis(void){
+}

FunctionSpace/Basis.h

+#ifndef _BASIS_H_
+#define _BASIS_H_
+
+/**
+   @class Basis
+   @brief Mother class of all Basis
+
+   This class is the @em mother (by @em inheritence) of all Basis.@n
+
+   A Basis is @em set of @em linearly @em independent Polynomial%s 
+   (or Vector%s of Polynomial%s).@n
+ */
+
+class Basis{
+ protected:
+  bool scalar;
+
+  int order;
+  int type;
+  int size;
+  int nodeNbr;
+  int dim;
+
+ public:
+  virtual ~Basis(void);
+  
+  int getOrder(void) const;
+  int getType(void) const;
+  int getSize(void) const;
+  int getNodeNbr(void) const;
+  int getDim(void) const;
+  
+  bool isScalar(void) const;
+
+ protected:
+  Basis(void);
+};
+
+/**
+   @fn Basis::~Basis(void)
+   @return Deletes this Basis
+
+   @fn int Basis::getOrder(void) const
+   @return Returns the @em polynomial @em order of the Basis
+
+   @fn int Basis::getType(void) const
+   @return Returns the @em type of the Basis:
+   @li 0 for 0-form
+   @li 1 for 1-form
+   @li 2 for 2-form
+   @li 3 for 3-form
+
+   @todo Check if the 'form numbering' is good
+   
+   @fn int Basis::getSize(void) const
+   @return Returns the number of Polynomial%s 
+   (or Vector%s of Polynomial%s) in the Basis
+
+   @fn int Basis::getNodeNbr(void) const
+   @return Returns the node number of 
+   the @em geometric @em reference @em element
+ 
+   @fn int Basis::getDim(void) const
+   @return Returns the @em dimension 
+   (1D, 2D or 3D) of the Basis
+
+   @fn bool Basis::isScalar(void) const
+   @return Returns:
+   @li @c true, if this is a 
+   @em scalar Basis (BasisScalar)
+   @li @c false, if this is a
+   @em vectorial Basis (BasisVector)
+
+   @note
+   Scalar basis are sets of 
+   Polynomial%s@n
+   Vectorial basis are sets of 
+   Vector%s of Polynomial%s
+*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline int Basis::getOrder(void) const{
+  return order;
+}
+
+inline int Basis::getType(void) const{
+  return type;
+}
+
+inline int Basis::getSize(void) const{
+  return size;
+}
+
+inline int Basis::getNodeNbr(void) const{
+  return nodeNbr;
+}
+
+inline int Basis::getDim(void) const{
+  return dim;
+}
+
+inline bool Basis::isScalar(void) const{
+  return scalar;
+}
+
+#endif

FunctionSpace/BasisScalar.cpp

+#include "BasisScalar.h"
+
+BasisScalar::BasisScalar(void){
+  scalar = true;
+}
+
+BasisScalar::~BasisScalar(void){
+}

FunctionSpace/BasisScalar.h

+#ifndef _BASISSCALAR_H_
+#define _BASISSCALAR_H_
+
+#include "Basis.h"
+#include "Polynomial.h"
+
+/**
+   @class BasisScalar
+   @brief Mother class of all 
+   @em scalar Basis
+
+   This class is the @em mother (by @em inheritence) 
+   of all @em scalar Basis.@n
+*/
+
+class BasisScalar: public Basis{
+ protected:
+  Polynomial* basis;
+
+ public:
+  virtual ~BasisScalar(void);
+
+  Polynomial* getBasis(void) const;
+
+ protected:
+  BasisScalar(void);
+};
+
+/**
+   @fn BasisScalar::~BasisScalar
+   @return Deletes this Basis
+
+   @fn BasisScalar::getBasis
+   @return Returns the set of 
+   @em Polynomial%s
+   that defines this Basis
+*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline Polynomial* BasisScalar::getBasis(void) const{
+  return basis;
+}
+
+#endif

FunctionSpace/BasisVector.cpp

+#include "BasisVector.h"
+
+BasisVector::BasisVector(void){
+  scalar = false;
+}
+
+BasisVector::~BasisVector(void){
+}
+

FunctionSpace/BasisVector.h

+#ifndef _BASISVECTOR_H_
+#define _BASISVECTOR_H_
+
+#include "Basis.h"
+#include "Polynomial.h"
+#include "Vector.h"
+
+/**
+   @class BasisVector
+   @brief Mother class of all 
+   @em vectorial Basis
+
+   This class is the @em mother (by @em inheritence) 
+   of all @em vectorial Basis.@n
+*/
+
+class BasisVector: public Basis{
+ protected:
+  Vector<Polynomial>* basis;
+
+ public:
+  virtual ~BasisVector(void);
+
+  Vector<Polynomial>* getBasis(void) const;
+
+ protected:
+  BasisVector(void);
+};
+
+/**
+   @fn BasisVector::~BasisVector
+   @return Deletes this Basis
+
+   @fn BasisVector::getBasis
+   @return Returns the set of 
+   @em Vector%s @em of @em Polynomial%s
+   that defines this Basis
+*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline Vector<Polynomial>* BasisVector::getBasis(void) const{
+  return basis;
+}
+
+#endif

FunctionSpace/CMakeLists.txt

+# Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+#
+# See the LICENSE.txt file for license information. Please report all
+# bugs and problems to <gmsh@geuz.org>.
+
+set(SRC
+  Matrix.cpp
+  VectorDouble.cpp
+  VectorPolynomial.cpp
+  Polynomial.cpp
+  Legendre.cpp
+
+  Basis.cpp
+  BasisScalar.cpp
+  BasisVector.cpp
+  QuadNodeBasis.cpp
+  QuadEdgeBasis.cpp
+  TriNodeBasis.cpp
+  TriEdgeBasis.cpp
+  TriNedelecBasis.cpp
+)
+
+file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) 
+append_gmsh_src(FunctionSpace "${SRC};${HDR}")
+
+## Compatibility with SmallFEM (TO BE REMOVED !!!)
+add_sources_in_gmsh(FunctionSpace "${SRC}")
\ No newline at end of file

FunctionSpace/Legendre.cpp

+#include "Legendre.h"
+
+Polynomial Legendre::legendre(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) - lMinus * (n / nPlus);
+}
+
+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) - 
+    lMinus * Polynomial(1, 0, 2, 0) * ((    n    ) / nPlus);
+}
+
+Polynomial Legendre::integrated(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) - lMinus * ((n - 2) / nPlus);
+}
+
+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) - 
+    lMinus * Polynomial(1, 0, 2, 0) * ((    n - 2) / nPlus);
+}
+
+
+
+void Legendre::legendre(Polynomial* polynomial, const int order){
+  int i, j, k;
+
+  if(order >= 0)
+    polynomial[0] = Polynomial(1, 0, 0, 0);
+  
+  if(order >= 1)
+    polynomial[1] = Polynomial(1, 1, 0, 0);
+
+  if(order >= 2){
+    for(k = 2; k <= order; k++){
+      i = k - 1;
+      j = k - 2;
+      
+      polynomial[k] = legendre(i, polynomial[i], polynomial[j]);
+    }
+  }
+}
+
+void Legendre::integrated(Polynomial* polynomial, const int order){
+  int i, j, k;
+
+  if(order >= 1)
+    polynomial[0] = Polynomial(1, 1, 0, 0);
+
+  if(order >= 2)
+    polynomial[1] = (Polynomial(1, 2, 0, 0) + Polynomial(-1, 0, 0, 0)) * 0.5;
+
+  if(order >= 3){
+    for(k = 2; k < order; k++){
+      i = k - 1;
+      j = k - 2;
+      
+      polynomial[k] = integrated(k, polynomial[i], polynomial[j]);
+    }
+  }
+}
+
+void Legendre::scaled(Polynomial* polynomial, const int order){
+  int i, j, k;
+
+  if(order >= 0)
+    polynomial[0] = Polynomial(1, 0, 0, 0);
+  
+  if(order >= 1)
+    polynomial[1] = Polynomial(1, 1, 0, 0);
+
+  if(order >= 2){
+    for(k = 2; k <= order; k++){
+      i = k - 1;
+      j = k - 2;
+      
+      polynomial[k] = scaled(i, polynomial[i], polynomial[j]);
+    }
+  }
+}
+
+void Legendre::intScaled(Polynomial* polynomial, const int order){
+  int i, j, k;
+
+  if(order >= 1)
+    polynomial[0] = Polynomial(1, 1, 0, 0);
+
+  if(order >= 2)
+    polynomial[1] = (Polynomial(1, 2, 0, 0) + Polynomial(-1, 0, 2, 0)) * 0.5;
+
+  if(order >= 3){
+    for(k = 2; k < order; k++){
+      i = k - 1;
+      j = k - 2;
+      
+      polynomial[k] = intScaled(k, polynomial[i], polynomial[j]);
+    }
+  }
+}
+
+/*
+#include <iostream>
+int main(void){
+  Polynomial p[13];
+  Legendre::legendre(p, 12);
+
+  for(int i = 0; i < 13; i++)
+    std::cout << "l(" << i << ")= " << (p[i]).toString() << std::endl;
+
+
+  std::cout << std::endl;
+
+  Polynomial pi[12];
+  Legendre::integrated(pi, 12);
+
+  for(int i = 0; i < 12; i++)
+    std::cout << "L(" << i + 1 << ")= " << pi[i].toString() << std::endl;
+
+
+  std::cout << std::endl;
+
+  Polynomial ps[13];
+  Legendre::scaled(ps, 12);
+
+  for(int i = 0; i < 13; i++)
+    std::cout << "ls(" << i << ")= " << (ps[i]).toString() << std::endl;
+
+
+  std::cout << std::endl;
+
+  Polynomial pis[12];
+  Legendre::intScaled(pis, 12);
+
+  for(int i = 0; i < 12; i++)
+    std::cout << "Ls(" << i + 1<< ")= " << (pis[i]).toString() << std::endl;
+
+
+  Polynomial p1 = Polynomial(1, 1, 1, 0);// - Polynomial(1, 0, 0, 0);
+  Polynomial p2 = Polynomial(1, 0, 1, 1);// - Polynomial(1, 0, 0, 0);
+
+  Polynomial p3 = pis[11].compose(p1, p2);
+
+  std::cout << std::endl << "p1 = " << p1.toString() << std::endl;
+  std::cout << "p2 = " << p2.toString() << std::endl;
+  std::cout << "Ls(12, p1, p2) = " << p3.toString() << std::endl;
+
+  return 0;
+}
+*/

FunctionSpace/Legendre.h

+#ifndef _LEGENDRE_H_
+#define _LEGENDRE_H_
+
+#include "Polynomial.h"
+
+/**
+   @class Legendre
+   @brief Generators for legendre Polynomial%s
+   
+   This class handles the generation of legendre 
+   Polynomial%s of many types:
+   @li Classical legendre (Legendre::legendre)
+   @li Integrated legendre (Legendre::integrated)
+   @li Scaled legendre (Legendre::scaled)
+   @li Integrated Scaled legendre (Legendre::intScaled)
+
+   @note
+   It is @em not @em requiered to instantiate a Legendre class.@n
+   Indeed, all its methods are @em static.@n
+   Each method is actualy a @em specific Polynomial @em generator.
+ */
+
+class Legendre{
+ public:
+  static void legendre(Polynomial* polynomial, const int order);
+  static void integrated(Polynomial* polynomial, const int order);
+  static void scaled(Polynomial* polynomial, const int order);
+  static void intScaled(Polynomial* polynomial, const int order);
+
+ private:
+  static Polynomial legendre(const int n,
+			     const Polynomial& l, 
+			     const Polynomial& lMinus);
+  
+  static Polynomial integrated(const int n,
+			       const Polynomial& l, 
+			       const Polynomial& lMinus);
+
+  static Polynomial scaled(const int n,
+			   const Polynomial& l, 
+			   const Polynomial& lMinus);
+
+  static Polynomial intScaled(const int n,
+			      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) 
+   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']  
+
+   @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']  
+
+   @fn void Legendre::scaled(Polynomial*, const int)
+   @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']  
+
+   @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']  
+ */
+
+#endif

FunctionSpace/Matrix.cpp

+#include <sstream>
+#include "Matrix.h"
+#include "Exception.h"
+
+extern "C"{
+#include <cblas.h>
+}
+
+Matrix::Matrix(const int n, const int m){
+  if(!(n || m))
+    throw Exception("Matrix can't be of dimension (0, 0)");
+    
+  nRow = n;
+  nCol = m;
+  nElem = nRow * nCol;
+  matrix = new double[nElem];
+}
+
+Matrix::~Matrix(void){
+  if(matrix)
+    delete[] matrix;
+}
+
+void Matrix::allToZero(void){
+  for(int i = 0; i < nElem; i++)
+    matrix[i] = 0.0;
+}
+
+Vector<double> Matrix::mult(const Vector<double>& v) const{
+  Vector<double> s(v.N);
+  
+  cblas_dgemv(CblasRowMajor, CblasNoTrans, 
+	      nRow, nCol, 
+	      1.0, matrix, nCol,
+	      v.v, 1,
+	      0.0, s.v, 1);
+
+  return s;
+}
+
+std::string Matrix::toString(void) const{
+  std::stringstream s;
+  
+  for(int i = 0; i < nRow; i++){
+    for(int j = 0; j < nCol; j++){
+      s << std::scientific << std::showpos << get(i, j) << "\t";
+    }
+    s << std::endl;
+  }
+
+  return s.str();
+}
+
+std::string Matrix::toStringMatlab(void) const{
+  std::stringstream s;
+  
+  s << "[";
+  for(int i = 0; i < nRow; i++){
+    for(int j = 0; j < nCol; j++){
+      s << std::scientific << std::showpos << get(j, i) << " ";
+    }
+    s << ";";
+  }
+  s << "]";
+  
+  return s.str();
+}

FunctionSpace/Matrix.h

+#ifndef _MATRIX_H_
+#define _MATRIX_H_
+
+#include <string>
+#include "Vector.h"
+
+/**
+   @class Matrix
+   @brief Handles matrices
+   
+   This class represents a @c n by @c m matrix.
+
+   @todo
+   Other than Matrix of double
+ */
+
+class Solver;
+
+class Matrix{
+ private:
+  int nRow;
+  int nCol;
+  int nElem;
+  double *matrix;
+  friend class Solver;
+
+ public:
+   Matrix(const int n, const int m);
+  ~Matrix(void);
+
+  int row(void) const;
+  int col(void) const;
+
+  void   set(const int i, const int j, const double v);
+  void   allToZero(void);
+  double get(const int i, const int j) const;
+  
+  double& operator()(const int i, const int j);
+  double  operator()(const int i, const int j) const;
+  
+  Vector<double> mult(const Vector<double>& v) const;
+
+  std::string toString(void) const;
+  std::string toStringMatlab(void) const;
+};
+
+/**
+   @fn Matrix::Matrix
+   @param n The number of @em rows of the future Matrix
+   @param m The number of @em columns of the future Matrix
+   @return Returns a new @c n by @c m Matrix
+
+   @fn Matrix::~Matrix
+   @return Deletes this Matrix
+
+   @fn Matrix::row
+   @return Returns the number of @c rows of this Matrix
+   
+   @fn Matrix::col
+   @return Returns the number of @c columns of this Matrix
+
+   @fn Matrix::set
+   @param i A @em row index
+   @param j A @em column index
+   @param v A value
+   @returns Sets the given value at the position (@c i, @c j)
+   in this Matrix
+
+   @fn Matrix::allToZero
+   @returns Sets all the Matrix elements to @em @c 0
+ 
+   @fn Matrix::get
+   @param i A @em row index
+   @param j A @em column index
+   @returns Retuns the value at the position (@c i, @c j)
+   in this Matrix
+   
+   @fn double& Matrix::operator()(const int, const int)
+   @param i A @em row index
+   @param j A @em column index
+   @return Returns a @em reference to the element
+   at position (@c i, @c j) in this Matrix
+
+   @fn double Matrix::operator()(const int, const int) const
+   @param i A @em row index
+   @param j A @em column index
+   @return Returns the @em value of the element
+   at position (@c i, @c j) in this Matrix
+
+   @fn Matrix::mult
+   @param v A Vector of size @c m (for a @c n by @c m Matrix)
+   @return Returns the Vector (of size @c n) resulting 
+   of the product of @em this Matrix with the @em given Vector
+
+   @fn Matrix::toString
+   @return Retuns a string correspondig to this Matrix
+
+   @fn Matrix::toStringMatlab
+   @return Same as Matrix::toString, but the string is in 
+   Matlab / Octave format
+*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline int Matrix::row(void) const{
+  return nRow;
+}
+
+inline int Matrix::col(void) const{
+  return nCol;
+}
+
+inline void Matrix::set(const int i, const int j, const double v){
+  matrix[j + i * nCol] = v;
+}
+
+inline double Matrix::get(const int i, const int j) const{
+  return matrix[j + i * nCol];
+}
+
+inline double& Matrix::operator()(const int i, const int j){
+  return matrix[j + i * nCol];
+}
+
+inline double Matrix::operator()(const int i, const int j) const{
+  return matrix[j + i * nCol];
+}
+
+#endif

FunctionSpace/Polynomial.cpp

+#include <cmath>
+#include <sstream>
+#include <stack>
+#include "Polynomial.h"
+
+using namespace std;
+
+const char Polynomial::coefName[3] = {'x', 'y', 'z'};
+
+Polynomial::Polynomial(const double coef, const int powerX,
+ 	                                  const int powerY,
+                                          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); 
+}
+
+Polynomial::Polynomial(void){
+  nMon = 0;
+  mon  = NULL;
+}
+
+Polynomial::~Polynomial(void){
+  if(mon)
+    delete[] mon;
+}
+
+void Polynomial::derivative(const int dim){ 
+  // Take derivative //
+  for(int i = 0; i < nMon; i++){
+    mon[i].coef *= mon[i].power[dim];
+    mon[i].power[dim] -= 1;
+  }
+
+  // Remove zero monomials //
+  int N = 0;
+  stack<monomial_t*> s;
+
+  for(int i = 0; i < nMon; i++){
+    if(mon[i].coef != 0.0){
+      s.push(&mon[i]);
+      N++;
+    }
+  }
+
+  // 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) 
+    return;
+
+  // Else, remove them //
+  monomial_t* tmp = mon;
+  
+  mon  = new monomial_t[N];
+  nMon = N;
+
+  for(int i = N - 1; i >= 0; i--){
+    mon[i] = *(s.top());
+    s.pop();
+  }
+
+  delete[] tmp;
+
+  // Sort resulting monomial and return // 
+  sort(mon, nMon);
+
+  return;
+}
+
+Vector<Polynomial> Polynomial::gradient(void) const{
+  Vector<Polynomial> grad(3);
+
+  // Copy Polynomial //
+  grad(0) = *this;
+  grad(1) = *this;
+  grad(2) = *this;
+
+  // Derivative with respect to each direction //
+  grad(0).derivative(0);
+  grad(1).derivative(1);
+  grad(2).derivative(2);
+  
+  return grad;
+}
+
+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]) 
+                       * pow(z, mon[i].power[2]);
+  }
+
+  return val;
+}
+
+Polynomial Polynomial::operator+(const Polynomial& other) const{
+  Polynomial newP;
+  
+  newP.nMon = mergeMon(mon, nMon, other.mon, other.nMon, &newP.mon);
+
+  return newP;
+}
+
+Polynomial Polynomial::operator-(const Polynomial& other) const{
+  const int   otherNMon  = other.nMon;
+  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);   
+  
+  Polynomial  newP;
+
+  mult(otherMinus, otherNMon, -1); 
+  newP.nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &newP.mon);
+
+  delete[] otherMinus;
+  return newP;
+}
+
+Polynomial Polynomial::operator*(const Polynomial& other) const{
+  Polynomial newP;
+  
+  newP.nMon = mult(mon, nMon, other.mon, other.nMon, &newP.mon);
+
+  return newP;
+}
+
+Polynomial Polynomial::operator*(const double alpha) const{
+  Polynomial newP;
+  
+  newP.mon  = copyMonomial(mon, nMon);
+  newP.nMon = nMon;
+
+  newP.mul(alpha);
+
+  return newP;
+}
+
+
+void Polynomial::add(const Polynomial& other){
+  monomial_t* tmp = mon;
+  
+  nMon = mergeMon(mon, nMon, other.mon, other.nMon, &mon);
+
+  delete[] tmp;
+}
+
+void Polynomial::sub(const Polynomial& other){
+  const int   otherNMon  = other.nMon;
+  monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);   
+  monomial_t* tmp        = mon;
+  
+  mult(otherMinus, otherNMon, -1); 
+  nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &mon);
+
+  delete[] otherMinus;
+  delete[] tmp;
+}
+
+void Polynomial::mul(const Polynomial& other){
+  monomial_t* tmp = mon;
+  
+  nMon = mult(mon, nMon, other.mon, other.nMon, &mon);
+
+  delete[] tmp;
+}
+
+void Polynomial::mul(const double alpha){
+  for(int i = 0; i < nMon; i++)
+    mon[i].coef *= alpha;
+}
+
+void Polynomial::power(const int n){
+  if (n < 0)
+    return;
+
+  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;
+  }
+}
+
+Polynomial Polynomial::compose(const Polynomial& other) const{
+  stack<monomial_t> stk;
+
+  for(int i = 0; i < nMon; i++)
+    compose(&mon[i], other, &stk);
+
+  return polynomialFromStack(stk);
+}
+
+Polynomial Polynomial::compose(const Polynomial& otherA,
+			       const Polynomial& otherB) const{
+  stack<monomial_t> stk;
+
+  for(int i = 0; i < nMon; i++)
+    compose(&mon[i], otherA, otherB, &stk);
+
+  return polynomialFromStack(stk);
+}
+
+void Polynomial::operator=(const Polynomial& other){
+  if(mon)
+    delete[] mon;
+  
+  nMon = other.nMon;
+  mon  = copyMonomial(other.mon, nMon);
+}
+
+string Polynomial::toString(const Polynomial::monomial_t* mon, const bool isAbs){
+  stringstream stream;
+  const bool minusOne    = mon->coef == -1.0;
+  const bool notUnitCoef = mon->coef !=  1.0 && !minusOne;
+
+  // If we have a constant term
+  if(!mon->power[0] && !mon->power[1] && !mon->power[2]){
+    stream << mon->coef;
+    return stream.str();
+  }
+
+  // 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);
+
+  if(notUnitCoef && !isAbs)
+    stream << mon->coef;
+
+  if(minusOne && !isAbs)
+    stream << "-";
+
+  // We look for each power
+  bool notOnce = false;
+  for(int i = 0; i < 3; i++){
+    // If we have a non zero power, we display it
+    if(mon->power[i]){
+      if(notUnitCoef || notOnce)
+	stream << " * ";
+
+      stream << coefName[i];
+      
+      if(mon->power[i] != 1)
+	stream << "^" << mon->power[i];
+
+      notOnce = true;
+    }
+  }
+
+  return stream.str();
+}
+
+string Polynomial::toString(void) const{
+  stringstream stream;
+  bool isAbs = false;
+
+  stream << toString(&mon[0], isAbs);
+    
+  for(int i = 1; i < nMon; i++){
+    if(mon[i].coef < 0.0){
+      stream << " - ";
+      isAbs = true;
+    }
+
+    else{
+      stream << " + ";
+      isAbs = false;
+    }
+
+    stream << toString(&mon[i], isAbs);
+  }
+
+  return stream.str();
+}
+
+bool Polynomial::isSmaller(const Polynomial::monomial_t* a, 
+			   const Polynomial::monomial_t* b){
+  // 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];
+
+    for(int i = 0; i < 3; i++)
+      if(dif[i])
+	last = dif[i];
+
+    if(last < 0)
+      return true;
+    
+    else
+      return false;
+  }
+    
+  return false;
+}
+
+void Polynomial::sort(monomial_t* mon, const int size){
+  for(int i = 0; i < size; i++)
+    for(int j = i; j < size; j++)
+      if(isSmaller(&mon[j], &mon[i]))
+        swap(mon, i, j);
+}
+
+void Polynomial::swap(monomial_t* mon, const int i, const int j){
+  monomial_t tmp = mon[i];
+  mon[i] = mon[j];
+  mon[j] = tmp;
+}
+
+
+int Polynomial::mergeMon(monomial_t* sourceA, const int sizeA, 
+			 monomial_t* sourceB, const int sizeB,
+			 monomial_t** dest){
+  stack<monomial_t> s;
+  monomial_t tmp;
+
+  int i = 0;
+  int j = 0;
+  int N = 0;
+
+  while(i < sizeA && j < sizeB){
+    if(sourceA[i].coef == 0.0)
+      i++;
+
+    else if(sourceB[j].coef == 0.0)
+      j++;
+    
+    else if(isEqual(&sourceA[i], &sourceB[j])){
+      tmp       = sourceA[i];
+      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++;
+      N++;
+    }
+
+    else{
+      s.push(sourceB[j]);
+      j++;
+      N++;
+    }
+  }
+
+  while(i == sizeA && j < sizeB){
+    s.push(sourceB[j]);
+    j++;
+    N++;
+  }
+
+  while(i < sizeA && j == sizeB){
+    s.push(sourceA[i]);
+    i++;
+    N++;
+  }
+
+  if(!N){
+    *dest = zeroPolynomial();
+    N++;
+  }
+
+  else{
+    *dest = new monomial_t[N];
+
+    for(int k = N - 1; k >= 0; k--){
+      (*dest)[k] = s.top();
+      s.pop();
+    }
+  }
+
+  return N;
+}
+
+int Polynomial::mult(const monomial_t* sourceA, const int sizeA,
+		     const monomial_t* sourceB, const int sizeB,
+		     monomial_t** dest){
+
+  const monomial_t* a; // smaller polynomial
+  const monomial_t* b; // bigger polynomial
+  int nDist;
+  int size;
+
+  if(sizeA < sizeB){
+    a     = sourceA;
+    b     = sourceB; 
+    nDist = sizeA;
+    size  = sizeB;
+  }
+  
+  else{
+    a     = sourceB;
+    b     = sourceA;
+    nDist = sizeB;
+    size  = sizeA;
+  }
+
+  // Check if zero //
+  if(a[0].coef == 0 || b[0].coef == 0){
+    *dest = zeroPolynomial();
+    return 1;
+  }
+
+  // Distrubute all monomials //
+  monomial_t** dist = new monomial_t*[nDist];
+
+  for(int i = 0; i < nDist; i++){
+    dist[i] = copyMonomial(b, size);
+  
+    distribute(dist[i], size, &a[i]);
+  }
+
+  // Merge //
+  int finalSize    = size;
+  int nDistMinus   = nDist - 1;
+  monomial_t** tmp = new monomial_t*[nDistMinus]; // Temp array for all dist[0];
+
+  for(int i = 1, j = 0; i < nDist; i++, j++){
+    tmp[j] = dist[0];
+    
+    finalSize = mergeMon(dist[0], finalSize, 
+			 dist[i], size,
+			 &dist[0]);
+  }
+
+  // Keep distributed polynomial //
+  *dest = dist[0];
+
+  // Free Temporary Resources and Return //
+  for(int i = 1, j = 0; i < nDist; i++, j++){
+    delete[] dist[i];
+    delete[] tmp[j];
+  }
+  
+  delete[] dist;
+  delete[] tmp;
+
+  return finalSize;
+}
+
+void Polynomial::mult(monomial_t* source, const int size, const double alpha){
+  for(int i = 0; i < size; i++)
+    source[i].coef *= alpha;
+}
+
+
+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];    
+  }
+}
+
+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]);
+    }
+  }
+}
+
+void Polynomial::compose(const Polynomial::monomial_t* src,
+			 Polynomial compA, Polynomial compB,
+			 stack<Polynomial::monomial_t>* stk){
+
+  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]);
+    }
+  }
+}
+
+Polynomial Polynomial::polynomialFromStack(std::stack<Polynomial::monomial_t>& stk){
+  Polynomial  newP;
+  monomial_t* tmp;
+  monomial_t* newMon;
+  int         newNMon;
+
+  if(!stk.size()){
+    newMon  = zeroPolynomial();
+    newNMon = 1;
+  }
+
+  else{
+    newMon    = new monomial_t[1];
+    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;
+    }
+  }
+
+  newP.nMon = newNMon;
+  newP.mon  = newMon;
+
+  return newP;
+}
+
+Polynomial::monomial_t* Polynomial::copyMonomial(const monomial_t* src, const int size){
+  monomial_t* dest = new monomial_t[size];
+
+  for(int i = 0; i < size; i++){
+    dest[i].coef     = src[i].coef;
+    dest[i].power[0] = src[i].power[0];
+    dest[i].power[1] = src[i].power[1];
+    dest[i].power[2] = src[i].power[2];
+  }
+
+  return dest;
+}
+
+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;
+}
+
+
+
+
+
+
+
+/*
+#include <iostream>
+int main(void){
+  Polynomial m0(-1  , 0, 0, 0);
+  Polynomial m1(4.2, 1, 0, 0);
+  Polynomial m2(4.2, 1, 1, 0);
+  Polynomial m3(4.2, 1, 0, 1);
+
+  cout << "m0 = " << m0.toString() << endl;
+  cout << "m1 = " << m1.toString() << endl;
+  cout << "m1(4, 5, 6) = " << m1.at(4, 5, 6) << endl;
+  cout << "m2 = " << m2.toString() << endl;
+  cout << "m3 = " << m3.toString() << endl;
+
+  Polynomial p0 = m1 + m0;
+  Polynomial p1 = m3 + m3;
+  Polynomial p2 = p1 + p0;
+
+  cout << "p0 = " << p0.toString() << endl;
+  cout << "p1 = " << p1.toString() << endl;
+  cout << "p2 = " << p2.toString() << endl;
+  cout << "p2(1.1, 2.2, 3.3) = " << p2.at(1.1, 2.2, 3.3) << endl;
+
+  p2.add(p0);
+  cout << "p2 = " << p2.toString() << endl;
+
+
+  Polynomial p3 = p2 * p0;
+  Polynomial p4 = p3 * p3;
+  cout << "p3 = " << p3.toString() << endl;
+  cout << "p4 = " << p4.toString() << endl;
+  
+  p3.mul(p3);
+  cout << "p3 = " << p3.toString() << endl;
+
+  Polynomial p5 = p2 - p3;
+  Polynomial p6 = p3 - p2;
+  Polynomial p7 = p6 - p6;
+  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);
+  Polynomial p10(p4);
+
+  p8.derivative(0);
+  p9.derivative(1);
+  p10.derivative(2);
+
+  cout << "p8 = "  << p8.toString()  << endl;
+  cout << "p9 = "  << p9.toString()  << endl;
+  cout << "p10 = " << p10.toString() << endl;
+
+  Polynomial p11 = p8 * 4;
+  cout << "p11 = " << p11.toString() << endl;
+  cout << "p11(1, 2, 3) = " << p11.at(1, 2, 3.1) << endl;
+
+  Polynomial m4(1, 1, 0, 0);
+  Polynomial m5(1, 0, 1, 0);
+  Polynomial m6(1, 0, 0, 1);
+
+  Polynomial p12 = m4 + m5 + m6;
+
+  cout << "p12 = " << p12.toString() << endl;
+
+  p12.power(4);
+  cout << "p12^4 = " << p12.toString() << endl;
+
+ Polynomial m7(1, 1, 0, 0);
+ Polynomial m8(2, 2, 0, 0);
+ Polynomial m9(3, 3, 0, 0);
+
+ Polynomial p13 = m9 + m7;
+ cout << "p13 = " << p13.toString() << endl;
+
+ Polynomial m10(1, 1, 1, 0);
+ Polynomial m11(2, 2, 1, 2);
+ Polynomial m12(1, 1, 3, 1);
+
+ Polynomial p14 = m10 + m11 + m12;
+ cout << "p14 = " << p14.toString() << endl;
+
+ Polynomial p15 = p14.compose(p13);
+ cout << "p15 = p14(p13) = " << p15.toString() << endl;
+
+ return 0;
+}
+
+*/

FunctionSpace/Polynomial.h

+#ifndef _POLYNOMIAL_H_
+#define _POLYNOMIAL_H_
+
+#include <string>
+#include <stack>
+#include "Vector.h"
+
+/**
+   @class Polynomial
+   @brief Represents @c 3D polynomials
+
+   This class represents @c 3D polynomials.@n
+   A @c 3D polynomials uses monomials of '@c xyz' type.
+ */
+
+// We suppose 3D Polynomial
+class Polynomial{
+ private: 
+  static const char coefName[3];
+
+  struct monomial_t{
+    double coef;
+    int power[3];
+  };
+  
+  int         nMon;
+  monomial_t*  mon;
+
+ public:
+   Polynomial(const double coef, const int powerX,
+ 	                         const int powerY,
+                                 const int powerZ);
+
+   Polynomial(const Polynomial& other);
+   Polynomial(void);
+  ~Polynomial(void);
+
+  void               derivative(const int dim);
+  Vector<Polynomial> gradient(void) const;
+
+  double operator()
+    (const double x, const double y, const double z) const;  
+  
+  double at
+    (const double x, const double y, const double z) const;  
+
+
+  Polynomial operator+(const Polynomial& other) const;
+  Polynomial operator-(const Polynomial& other) const;
+  Polynomial operator*(const Polynomial& other) const;
+  Polynomial operator*(const double alpha) const;
+
+  void add(const Polynomial& other);
+  void sub(const Polynomial& other);
+  void mul(const Polynomial& other);
+  void mul(const double alpha);
+
+  void power(const int n);
+
+  Polynomial compose(const Polynomial& other) const;
+  Polynomial compose(const Polynomial& otherA, 
+		     const Polynomial& otherB) const;
+
+  void operator=(const Polynomial& other);
+
+  std::string toString(void) const;
+
+ 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 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, 
+		      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, 
+		      Polynomial comp,
+		      std::stack<monomial_t>* stk);
+
+  static void compose(const monomial_t* src, 
+		      Polynomial compA, Polynomial compB,
+		      std::stack<monomial_t>* stk);
+
+  static Polynomial polynomialFromStack(std::stack<monomial_t>& stk);
+
+  static monomial_t* copyMonomial(const monomial_t* src, const int size);
+
+  static monomial_t* zeroPolynomial(void);
+  static monomial_t* unitPolynomial(void);
+};
+
+/**
+   @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 
+   of the futur monomial
+   @param powerY The power of the '@c y' coordinate 
+   of the futur monomial
+   @param powerZ The power of the '@c z' coordinate 
+   of the futur monomial
+   @return Returns a new Monomial with the given
+   parameters
+   @note
+   Note that Monomials are special case of Polynomial%s
+
+   @fn Polynomial::Polynomial(const Polynomial&)
+   @param other A Polynomial
+   @return Returns a new Polynomial, which is
+   the @em copy of the given one
+
+   @fn Polynomial::Polynomial(void)
+   @return Returns a new Polynomial, which is
+   @em empty
+   @warning
+   An @em empty Polynomial means: @em A @em 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, 
+   @em @c 0.
+
+   @fn Polynomial::~Polynomial
+   @return Deletes this Polynomial
+ 
+   @fn Polynomial::derivative
+   @param dim The dimention to use for the
+   derivation
+   @returns Derivates this Polynomial with 
+   respect to the given dimention
+   @note
+   Dimention:
+   @li @c 0 is for the @c x coordinate
+   @li @c 1 is for the @c y coordinate
+   @li @c 2 is for the @c z coordinate
+   
+   @fn Polynomial::gradient
+   @return Returns a Vector with the gradient 
+   of this Polynomial
+
+   @fn double Polynomial::operator()(const double, const double, const double)
+   @param x A value
+   @param y A value
+   @param z A value
+   @return Returns the @em evaluation of this
+   Polynomial at (@c x, @c y, @c z)
+
+   @fn Polynomial::at
+   @param x A value
+   @param y A value
+   @param z A value
+   @return Returns the @em evaluation of this
+   Polynomial at (@c x, @c y, @c z)
+
+   @fn Polynomial Polynomial::operator+(const Polynomial&) const
+   @param other An other Polynomial
+   @return Returns a @em new Polynomial, which is the
+   @em sum of this Polynomial and the given one
+
+   @fn Polynomial Polynomial::operator-(const Polynomial&) const
+   @param other An other Polynomial
+   @return Returns a @em new Polynomial, which is the
+   @em difference of this Polynomial and the given one
+
+   @fn Polynomial Polynomial::operator*(const Polynomial&) const
+   @param other An other Polynomial
+   @return Returns a @em new Polynomial, which is the
+   @em product of this Polynomial and the given one
+
+   @fn Polynomial Polynomial::operator*(const double) const
+   @param alpha A value
+   @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 
+   @em added to this Polynomial
+   @note
+   The result of this operation is stored in 
+   this Polynomial  
+
+   @fn Polynomial::sub
+   @param other An other Polynomial
+   @return The given Polynomial is 
+   @em substracted to this Polynomial
+   @note
+   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 
+   @em multiplied with this Polynomial
+   @note
+   The result of this operation is stored in 
+   this Polynomial  
+
+   @fn void Polynomial::mul(const double)
+   @param alpha A value
+   @return This Polynomial is @em multiplied
+   by the given value
+   @note
+   The result of this operation is stored in 
+   this Polynomial  
+
+   @fn Polynomial::power
+   @param n A @em natural number
+   @return Takes this Polynomial to the power @c n
+
+   @fn Polynomial Polynomial::compose(const Polynomial&) const
+   @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
+   This method returns a @em new Polynomial,
+   representing @c P(Q(x, @c y, @c z), @c y, @c z)
+
+   @fn Polynomial Polynomial::compose(const Polynomial&, const Polynomial&) const
+   @param otherA An other Polynomial, 
+   called @c Q(x, @c y, @c z)
+   @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 
+   @c P(Q(x, @c y, @c z), @c R(x, @c y, @c z), @c z)
+
+   @fn void Polynomial::operator=(const Polynomial&)
+   @param other A Polynomial
+   @return Sets this Polynomial to a @em copy
+   of the given one
+
+   @fn Polynomial::toString
+   @return Returns a string representing this Polynomial
+*/
+
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline double Polynomial::operator() (const double x, 
+				      const double y, 
+				      const double z) const{
+  return at(x, y, z);
+}
+
+inline bool Polynomial::isEqual(const Polynomial::monomial_t* a,
+				const Polynomial::monomial_t* b){
+  return a->power[0] == b->power[0] &&
+         a->power[1] == b->power[1] &&
+         a->power[2] == b->power[2];
+}
+
+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] ;    
+}
+
+inline bool Polynomial::isEqualPower(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] ;    
+}
+
+#endif

FunctionSpace/QuadEdgeBasis.cpp

+#include "QuadEdgeBasis.h"
+#include "Legendre.h"
+
+QuadEdgeBasis::QuadEdgeBasis(const int order){
+  // Set Basis Type //
+  this->order = order;
+  
+  type    = 1;
+  size    = 2 * (order + 2) * (order + 1);
+  nodeNbr = 4;
+  dim     = 2;
+
+  // Alloc Temporary Space //
+  const int  orderPlus = order + 1;
+  Polynomial* legendre    = new Polynomial[orderPlus];
+  Polynomial* intLegendre = new Polynomial[orderPlus];
+
+  Polynomial* iLegendreX  = new Polynomial[orderPlus];
+  Polynomial* iLegendreY  = new Polynomial[orderPlus];
+  Polynomial* legendreX   = new Polynomial[orderPlus];
+  Polynomial* legendreY   = new Polynomial[orderPlus];
+
+  Polynomial* lagrange    = new Polynomial[4];
+  Polynomial* lagrangeSum = new Polynomial[4];
+
+  Polynomial* lifting     = new Polynomial[4];
+  Polynomial* liftingSub  = new Polynomial[4];
+
+  // Integrated and classical Legendre Polynomial //
+  Legendre::integrated(intLegendre, orderPlus);
+  Legendre::legendre(legendre, order);
+
+  // Lagrange // 
+  lagrange[0] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  lagrange[1] = 
+    (Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  lagrange[2] = 
+    (Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 1, 0));
+
+  lagrange[3] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 1, 0));
+
+  // Lagrange Sum //
+  for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
+    lagrangeSum[i] = lagrange[i] + lagrange[j];
+    
+  // Lifting //
+  lifting[0] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  lifting[1] = 
+    (Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  lifting[2] = 
+    (Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 1, 0));
+
+  lifting[3] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 1, 0));
+
+  // Lifting Sub //
+  for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
+    liftingSub[i] = lifting[j] - lifting[i];
+
+
+  // Basis //
+  basis = new Vector<Polynomial>[size];
+  
+  // Edge Based (Nedelec) // 
+  int i = 0;
+  Polynomial oneHalf(0.5, 0, 0, 0);
+
+  for(int e = 0; e < 4; e++){
+    basis[i] = 
+      (liftingSub[e]).gradient();
+    
+    basis[i].mul(lagrangeSum[e]);
+    basis[i].mul(oneHalf);
+
+    i++;
+  }
+
+  // Edge Based (High Order) //
+  for(int l = 1; l < orderPlus; l++){
+    for(int e = 0; e < 4; e++){
+      basis[i] = 
+	(intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient();
+     
+      i++;
+    }
+  }
+
+  
+  // Cell Based (Preliminary) //
+  Polynomial px   = Polynomial(2, 1, 0, 0);
+  Polynomial py   = Polynomial(2, 0, 1, 0);
+  Polynomial zero = Polynomial(0, 0, 0, 0);
+
+  px = px - Polynomial(1, 0, 0, 0);
+  py = py - Polynomial(1, 0, 0, 0);
+
+  for(int l = 0; l < orderPlus; l++){
+    iLegendreX[l] = intLegendre[l].compose(px);
+    iLegendreY[l] = intLegendre[l].compose(py);
+     legendreX[l] =    legendre[l].compose(px);
+     legendreY[l] =    legendre[l].compose(py);
+  }
+
+  // Cell Based (Type 1) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l2 = 1; l2 < orderPlus; l2++){
+      basis[i] = (iLegendreX[l1] * iLegendreY[l2]).gradient();
+
+      i++;
+    }
+  }
+  
+  // Cell Based (Type 2) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l2 = 1; l2 < orderPlus; l2++){
+      basis[i](0) =  legendreX[l1] * iLegendreY[l2];
+      basis[i](1) = iLegendreX[l1] *  legendreY[l2] * -1;
+      basis[i](2) = zero;
+
+      i++;
+    }
+  }
+
+  // Cell Based (Type 3) //
+  for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
+    basis[i](0) = iLegendreY[l];
+    basis[i](1) = zero;
+    basis[i](2) = zero;
+
+    basis[iPlus](0) = zero;
+    basis[iPlus](1) = iLegendreX[l];
+    basis[iPlus](2) = zero;
+
+    i++;
+  }
+
+  // Free Temporary Sapce //
+  delete[] legendre;
+  delete[] intLegendre;
+
+  delete[] iLegendreX;
+  delete[] iLegendreY;
+  delete[] legendreX;
+  delete[] legendreY;
+
+  delete[] lagrange;
+  delete[] lagrangeSum;
+
+  delete[] lifting;
+  delete[] liftingSub;
+}
+
+QuadEdgeBasis::~QuadEdgeBasis(void){
+  delete[] basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+  const int P = 3;
+  const double d = 0.05;
+  const char x[2] = {'X', 'Y'};
+
+  QuadEdgeBasis b(P);
+  
+  const Vector<Polynomial>* basis = b.getBasis();
+  
+  printf("\n");
+  printf("clear all;\n");
+  printf("close 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 [rx ry] = p(i, x, y)\n");
+  printf("p = zeros(%d, 2);\n", b.getSize());
+  printf("\n");
+
+  for(int i = 0; i < b.getSize(); i++){
+    for(int j = 0; j < 2; j++)
+      printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i](j).toString().c_str());
+    //printf("p(%d) = %s", i, basis[i].toString().c_str());
+    printf("\n");
+  }
+
+  printf("\n");
+  printf("rx = p(i, 1);\n");
+  printf("ry = p(i, 2);\n");
+  printf("end\n");
+  printf("\n");
+  
+  printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
+  
+  for(int i = 0; i < b.getSize(); i++)
+    for(int j = 0; j < 2; j++)
+      printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
+
+  printf("\n");
+  printf("for i = 1:lx\n");
+  printf("for j = 1:ly\n");
+  printf("\n");
+
+  for(int i = 0; i < b.getSize(); i++)
+    printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
+
+  printf("\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 = 0; i < b.getSize(); i++)
+    printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
+  
+  printf("\n");
+
+  return 0;
+}
+*/

FunctionSpace/QuadEdgeBasis.h

+#ifndef _QUADEDGEBASIS_H_
+#define _QUADEDGEBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+   @class QuadEdgeBasis
+   @brief An Edge-Basis for Quads
+ 
+   This class can instantiate an Edge-Based Basis 
+   (high or low order) for Quads.@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 QuadEdgeBasis: public BasisVector{
+ public:
+  QuadEdgeBasis(const int order);
+  virtual ~QuadEdgeBasis(void);
+};
+
+/**
+   @fn QuadEdgeBasis::QuadEdgeBasis(const int order)
+   @param order The order of the Basis
+   @return Returns a new Edge-Basis for Quads of the given order
+
+   @fn QuadEdgeBasis::~QuadEdgeBasis(void)
+   @return Deletes this Basis
+*/
+
+#endif

FunctionSpace/QuadNodeBasis.cpp

+#include "QuadNodeBasis.h"
+#include "Legendre.h"
+
+QuadNodeBasis::QuadNodeBasis(const int order){
+  // Set Basis Type //
+  this->order = order;
+
+  type    = 0;
+  size    = (order + 1) * (order + 1);
+  nodeNbr = 4;
+  dim     = 2;
+
+  // Alloc Temporary Space //
+  Polynomial* legendre = new Polynomial[order];
+  Polynomial* lifting  = new Polynomial[4];
+
+  // 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));
+
+  lifting[1] = 
+    (Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  lifting[2] = 
+    (Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 1, 0));
+
+  lifting[3] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+    (Polynomial(1, 0, 1, 0));
+
+
+
+  // Basis //
+  basis = new 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));
+
+  basis[1] = 
+    (Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
+
+  basis[2] = 
+    (Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 1, 0));
+
+  basis[3] = 
+    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+    (Polynomial(1, 0, 1, 0));
+  
+  // Edge Based //
+  int i = 4;
+
+  for(int l = 1; l < order; l++){
+    for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){
+      basis[i] = 
+	legendre[l].compose(lifting[e2] - lifting[e1]) * (basis[e1] + basis[e2]);
+            
+      i++;
+    }
+  }
+
+  // Cell Based //
+  Polynomial px = Polynomial(2, 1, 0, 0);
+  Polynomial py = Polynomial(2, 0, 1, 0);
+
+  px = px - Polynomial(1, 0, 0, 0);
+  py = py - Polynomial(1, 0, 0, 0);
+
+  for(int l1 = 1; l1 < order; l1++){
+    for(int l2 = 1; l2 < order; l2++){
+      basis[i] = legendre[l1].compose(px) * legendre[l2].compose(py);
+
+      i++;
+    }
+  }
+
+  // Free Temporary Sapce //
+  delete[] legendre;
+  delete[] lifting;
+}
+
+QuadNodeBasis::~QuadNodeBasis(void){
+  delete[] basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+  const int P = 4;
+  const double d = 0.05;
+
+  QuadNodeBasis b(P);
+  
+  const 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)\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 = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
+  
+  for(int i = 0; i < b.getSize(); i++)
+    printf("p%d = zeros(lx, ly);\n", i + 1);
+
+  printf("\n");
+  printf("for i = 1:lx\n");
+  printf("for j = 1:ly\n");
+  printf("\n");
+
+  for(int i = 0; i < b.getSize(); i++)
+    printf("p%d(j, i) = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
+  
+  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;\ncontourf(x, y, p%d);\ncolorbar;\n", i, i);
+  
+  printf("\n");
+
+  return 0;
+}
+*/

FunctionSpace/QuadNodeBasis.h

+#ifndef _QUADNODEBASIS_H_
+#define _QUADNODEBASIS_H_
+
+#include "BasisScalar.h"
+
+/**
+   @class QuadNodeBasis
+   @brief A Node-Basis for Quads
+ 
+   This class can instantiate a Node-Based Basis 
+   (high or low order) for Quads.@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 QuadNodeBasis: public BasisScalar{
+ public:
+  QuadNodeBasis(const int order);
+  virtual ~QuadNodeBasis(void);
+};
+
+/**
+   @fn QuadNodeBasis::QuadNodeBasis(const int order)
+   @param order The order of the Basis
+   @return Returns a new Node-Basis for Quads of the given order
+
+   @fn QuadNodeBasis::~QuadNodeBasis(void)
+   @return Deletes this Basis
+*/
+
+#endif

FunctionSpace/TriEdgeBasis.cpp

+#include "TriEdgeBasis.h"
+#include "Legendre.h"
+
+TriEdgeBasis::TriEdgeBasis(const int order){
+  // Set Basis Type //
+  this->order = order;
+ 
+  type    = 1;
+  size    = (order + 1) * (order + 2);
+  nodeNbr = 3;
+  dim     = 2;
+
+  // Alloc Temporary Space //
+  const int orderPlus     = order + 1;
+  const int orderMinus    = order - 1;
+
+  Polynomial* legendre    = new Polynomial[orderPlus];
+  Polynomial* intLegendre = new Polynomial[orderPlus];
+  
+  Polynomial* lagrange    = new Polynomial[3];
+  Polynomial* lagrangeSub = new Polynomial[3];
+  Polynomial* lagrangeSum = new Polynomial[3];
+
+  Polynomial* u           = new Polynomial[orderPlus];
+  Polynomial* v           = new Polynomial[orderPlus];
+
+  // Classical and Intrated-Scaled Legendre Polynomial //
+  Legendre::legendre(legendre, order);
+  Legendre::intScaled(intLegendre, orderPlus);
+
+  // Lagrange //
+  lagrange[0] = 
+    Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
+
+  lagrange[1] = 
+    Polynomial(1, 1, 0, 0);
+
+  lagrange[2] = 
+    Polynomial(1, 0, 1, 0);
+
+  // Lagrange Sum //
+  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
+    lagrangeSum[i] = lagrange[j] + lagrange[i];
+
+  // Lagrange Sub //
+  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
+    lagrangeSub[i] = lagrange[i] - lagrange[j];
+
+
+  // Basis //
+  basis = new Vector<Polynomial>[size];
+
+  // Edge Based (Nedelec) //
+  int i = 0;
+
+  for(int j = 1; i < 3; j = (j + 1) % 3){
+    Vector<Polynomial> tmp = lagrange[j].gradient();
+    tmp.mul(lagrange[i]);
+
+    basis[i] = lagrange[i].gradient();
+    basis[i].mul(lagrange[j]);
+   
+    basis[i].sub(tmp);
+
+    i++;
+  }
+
+  // Edge Based (High Order) //
+  for(int l = 1; l < orderPlus; l++){
+    for(int e = 0; e < 3; e++){
+      basis[i] = 
+	(intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient();
+
+      i++;
+    }
+  }
+
+  // Cell Based (Preliminary) //
+  Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
+  
+  for(int l = 0; l < orderPlus; l++){
+    u[l] = intLegendre[l].compose(lagrangeSub[0] * -1, lagrangeSum[0]);
+    v[l] = legendre[l].compose(p);
+    v[l].mul(lagrange[2]);
+  }
+
+  // Cell Based (Type 1) //
+  for(int l1 = 1; l1 < order; l1++){
+    for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+      Vector<Polynomial> tmp = v[l2].gradient();
+      tmp.mul(u[l1]);
+
+      basis[i] = u[l1].gradient();
+      basis[i].mul(v[l2]);
+
+      basis[i].add(tmp);
+      
+      i++;
+    }
+  }
+
+  // Cell Based (Type 2) //
+  for(int l1 = 1; l1 < order; l1++){
+    for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+      Vector<Polynomial> tmp = v[l2].gradient();
+      tmp.mul(u[l1]);
+
+      basis[i] = u[l1].gradient();
+      basis[i].mul(v[l2]);
+
+      basis[i].sub(tmp);
+      
+      i++;
+    }
+  }
+
+  // Cell Based (Type 3) //
+  for(int l = 0; l < orderMinus; l++){
+    basis[i] = basis[0];
+    basis[i].mul(v[l]);
+    
+    i++;
+  }
+
+  // Free Temporary Sapce //
+  delete[] legendre;
+  delete[] intLegendre;
+  delete[] lagrange;
+  delete[] lagrangeSub;
+  delete[] lagrangeSum;
+  delete[] u;
+  delete[] v;
+}
+
+TriEdgeBasis::~TriEdgeBasis(void){
+  delete[] basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+  const int P = 1;
+  const double d = 0.05;
+  const char x[2] = {'X', 'Y'};
+
+  TriEdgeBasis b(P);
+  
+  const 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 [rx ry] = p(i, x, y)\n");
+  printf("p = zeros(%d, 2);\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("p(%d, 1) = %s;\n", i + 1, basis[i](0).toString().c_str());
+    printf("p(%d, 2) = %s;\n", i + 1, basis[i](1).toString().c_str());
+    printf("\n");
+  }
+  
+  printf("\n");
+  printf("rx = p(i, 1);\n");
+  printf("ry = p(i, 2);\n");
+  printf("end\n");
+  printf("\n");
+  
+  printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
+  
+  for(int i = 0; i < b.getSize(); i++)
+    for(int j = 0; j < 2; j++)
+      printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
+
+  printf("\n");
+  printf("for i = 1:lx\n");
+  printf("for j = 1:ly - i + 1\n");
+  printf("\n");
+
+  for(int i = 0; i < b.getSize(); i++)
+    printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
+  
+  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() - 1; i >= 0; i--)
+    printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
+  
+  printf("\n");
+  
+  return 0;
+}
+*/

FunctionSpace/TriEdgeBasis.h

+#ifndef _TRIEDGEBASIS_H_
+#define _TRIEDGEBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+   @class TriEdgeBasis
+   @brief An Edge-Basis for Triangles
+ 
+   This class can instantiate an Edge-Based Basis 
+   (high or low order) for Triangles.@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 TriEdgeBasis: public BasisVector{
+ public:
+  TriEdgeBasis(const int order);
+  virtual ~TriEdgeBasis(void);
+};
+
+
+/**
+   @fn TriEdgeBasis::TriEdgeBasis(const int order)
+   @param order The order of the Basis
+   @return Returns a new Edge-Basis for Triangles of the given order
+
+   @fn TriEdgeBasis::~TriEdgeBasis(void)
+   @return Deletes this Basis
+*/
+
+#endif

FunctionSpace/TriNedelecBasis.cpp

+#include "TriNedelecBasis.h"
+
+TriNedelecBasis::TriNedelecBasis(void){
+  // Set Basis Type //
+  order   = 1;
+  type    = 1;
+  size    = 3;
+  nodeNbr = 3;
+  dim     = 2;
+
+  // Lagrange //
+  Polynomial* lagrange    = new Polynomial[3];
+
+  lagrange[0] = 
+    Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
+
+  lagrange[1] = 
+    Polynomial(1, 1, 0, 0);
+
+  lagrange[2] = 
+    Polynomial(1, 0, 1, 0);
+
+  // Basis //
+  basis = new Vector<Polynomial>[size];
+
+  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
+    Vector<Polynomial> tmp = lagrange[j].gradient();
+    tmp.mul(lagrange[i]);
+
+    basis[i] = lagrange[i].gradient();
+    basis[i].mul(lagrange[j]);
+   
+    basis[i].sub(tmp);
+  }
+
+  // Free Temporary Sapce //
+  delete[] lagrange;
+}
+
+TriNedelecBasis::~TriNedelecBasis(void){
+  delete[] basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+  const int P = 1;
+  const double d = 0.05;
+  const char x[2] = {'X', 'Y'};
+
+  TriNedelec b;
+  
+  const 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 [rx ry] = p(i, x, y)\n");
+  printf("p = zeros(%d, 2);\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("p(%d, 1) = %s;\n", i + 1, basis[i](0).toString().c_str());
+    printf("p(%d, 2) = %s;\n", i + 1, basis[i](1).toString().c_str());
+    printf("\n");
+  }
+  
+  printf("\n");
+  printf("rx = p(i, 1);\n");
+  printf("ry = p(i, 2);\n");
+  printf("end\n");
+  printf("\n");
+  
+  printf("d = %lf;\nx = [0:d:1];\ny = x;\n\nlx = length(x);\nly = length(y);\n\n", d);
+  
+  for(int i = 0; i < b.getSize(); i++)
+    for(int j = 0; j < 2; j++)
+      printf("p%d%c = zeros(lx, ly);\n", i + 1, x[j]);
+
+  printf("\n");
+  printf("for i = 1:lx\n");
+  printf("for j = 1:ly - i + 1\n");
+  printf("\n");
+
+  for(int i = 0; i < b.getSize(); i++)
+    printf("[p%dX(j, i), p%dY(j, i)] = p(%d, x(i), y(j));\n", i + 1, i + 1, i + 1);
+  
+  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() - 1; i >= 0; i--)
+    printf("figure;\nquiver(x, y, p%dX, p%dY);\n", i + 1, i + 1);
+  
+  printf("\n");
+  
+  return 0;
+}
+*/