diff --git a/FunctionSpace/Basis.cpp b/FunctionSpace/Basis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..153bd0243b884937777374746fac8adec574a5b2
--- /dev/null
+++ b/FunctionSpace/Basis.cpp
@@ -0,0 +1,7 @@
+#include "Basis.h"
+
+Basis::Basis(void){
+}
+
+Basis::~Basis(void){
+}
diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
new file mode 100644
index 0000000000000000000000000000000000000000..c9094edb561e1e935c9f195fe96efc83b3671510
--- /dev/null
+++ b/FunctionSpace/Basis.h
@@ -0,0 +1,109 @@
+#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
diff --git a/FunctionSpace/BasisScalar.cpp b/FunctionSpace/BasisScalar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e902527ff83168f8c2fb0593af5f0422edb8bcbf
--- /dev/null
+++ b/FunctionSpace/BasisScalar.cpp
@@ -0,0 +1,8 @@
+#include "BasisScalar.h"
+
+BasisScalar::BasisScalar(void){
+ scalar = true;
+}
+
+BasisScalar::~BasisScalar(void){
+}
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
new file mode 100644
index 0000000000000000000000000000000000000000..753e743d2997f09f704f101bd34253b303725fe2
--- /dev/null
+++ b/FunctionSpace/BasisScalar.h
@@ -0,0 +1,47 @@
+#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
diff --git a/FunctionSpace/BasisVector.cpp b/FunctionSpace/BasisVector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..97c61e72098bc8775b3788d651b92bcd7a68582d
--- /dev/null
+++ b/FunctionSpace/BasisVector.cpp
@@ -0,0 +1,9 @@
+#include "BasisVector.h"
+
+BasisVector::BasisVector(void){
+ scalar = false;
+}
+
+BasisVector::~BasisVector(void){
+}
+
diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..a932f8e938626412655450fda845550cc870bd55
--- /dev/null
+++ b/FunctionSpace/BasisVector.h
@@ -0,0 +1,48 @@
+#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
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9b73fd520ec26848183e15cee588da099096ea5f
--- /dev/null
+++ b/FunctionSpace/CMakeLists.txt
@@ -0,0 +1,27 @@
+# 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
diff --git a/FunctionSpace/Legendre.cpp b/FunctionSpace/Legendre.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7ad4f9bbb1a8bebf5ed6d70bc0d910491c832bb
--- /dev/null
+++ b/FunctionSpace/Legendre.cpp
@@ -0,0 +1,163 @@
+#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;
+}
+*/
diff --git a/FunctionSpace/Legendre.h b/FunctionSpace/Legendre.h
new file mode 100644
index 0000000000000000000000000000000000000000..4120b65adfcdbe50af71b773a5bb354915f21550
--- /dev/null
+++ b/FunctionSpace/Legendre.h
@@ -0,0 +1,78 @@
+#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
diff --git a/FunctionSpace/Matrix.cpp b/FunctionSpace/Matrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..752345e14d506af4b8c32d5458caf692fe1304e0
--- /dev/null
+++ b/FunctionSpace/Matrix.cpp
@@ -0,0 +1,67 @@
+#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();
+}
diff --git a/FunctionSpace/Matrix.h b/FunctionSpace/Matrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..fe8425ff90b16887556a7d61dea326d8fd307e8a
--- /dev/null
+++ b/FunctionSpace/Matrix.h
@@ -0,0 +1,131 @@
+#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
diff --git a/FunctionSpace/Polynomial.cpp b/FunctionSpace/Polynomial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e277d748a7ad0f3b3c3e57c437c9c5c78c1ef63c
--- /dev/null
+++ b/FunctionSpace/Polynomial.cpp
@@ -0,0 +1,714 @@
+#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;
+}
+
+*/
diff --git a/FunctionSpace/Polynomial.h b/FunctionSpace/Polynomial.h
new file mode 100644
index 0000000000000000000000000000000000000000..02d4968b163f5175a9f902989b8255c089ab0c67
--- /dev/null
+++ b/FunctionSpace/Polynomial.h
@@ -0,0 +1,287 @@
+#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
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd072beb30d601685b3e078e0fab7425cc62acd9
--- /dev/null
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -0,0 +1,243 @@
+#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;
+}
+*/
diff --git a/FunctionSpace/QuadEdgeBasis.h b/FunctionSpace/QuadEdgeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f65d404c668ed4270e3fb9a94ee51adc678ffe4
--- /dev/null
+++ b/FunctionSpace/QuadEdgeBasis.h
@@ -0,0 +1,33 @@
+#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
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..73cecded561e54d8c0183e645558f2e6efb210da
--- /dev/null
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -0,0 +1,157 @@
+#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;
+}
+*/
diff --git a/FunctionSpace/QuadNodeBasis.h b/FunctionSpace/QuadNodeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..e0cbf3c89c67650d7861e5cc1f7f253da277be6d
--- /dev/null
+++ b/FunctionSpace/QuadNodeBasis.h
@@ -0,0 +1,33 @@
+#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
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2004cc53a61aaf5a5cd5ade78ba207b98f96dae
--- /dev/null
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -0,0 +1,208 @@
+#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;
+}
+*/
diff --git a/FunctionSpace/TriEdgeBasis.h b/FunctionSpace/TriEdgeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..12f1271792d4780d12e61b6548ffb78f092f5944
--- /dev/null
+++ b/FunctionSpace/TriEdgeBasis.h
@@ -0,0 +1,34 @@
+#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
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fd2805125c03e71694ab78659e4846d07fcf07eb
--- /dev/null
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -0,0 +1,113 @@
+#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;
+}
+*/
diff --git a/FunctionSpace/TriNedelecBasis.h b/FunctionSpace/TriNedelecBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..65179c78ad0f2fdf2bdaa699bf6ed198058fee89
--- /dev/null
+++ b/FunctionSpace/TriNedelecBasis.h
@@ -0,0 +1,29 @@
+#ifndef _TRINEDELECBASIS_H_
+#define _TRINEDELECBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+ @class TriNedelecBasis
+ @brief Nedelec Basis for Triangles
+
+ This class can instantiate a Nedelec Basis
+ for Triangles.@n
+*/
+
+class TriNedelecBasis: public BasisVector{
+ public:
+ TriNedelecBasis(void);
+ virtual ~TriNedelecBasis(void);
+};
+
+
+/**
+ @fn TriNedelecBasis::TriNedelecBasis
+ @return Returns a new Nedelec Basis for Triangles
+
+ @fn TriNedelecBasis::~TriNedelecBasis(void)
+ @return Deletes this Basis
+*/
+
+#endif
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef000f573acc80806be24c0c061d8599aa3c6ca7
--- /dev/null
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -0,0 +1,148 @@
+#include "TriNodeBasis.h"
+#include "Legendre.h"
+
+TriNodeBasis::TriNodeBasis(const int order){
+ // Set Basis Type //
+ this->order = order;
+
+ type = 0;
+ size = (order + 1) * (order + 2) / 2;
+ nodeNbr = 3;
+ dim = 2;
+
+ // Alloc Temporary Space //
+ Polynomial* legendre = new Polynomial[order];
+ Polynomial* intLegendre = new Polynomial[order];
+ Polynomial* lagrangeSub = new Polynomial[3];
+ Polynomial* lagrangeSum = new Polynomial[3];
+
+ // Classical and Intrated-Scaled Legendre Polynomial //
+ const int orderMinus = order - 1;
+
+ Legendre::legendre(legendre, orderMinus);
+ Legendre::intScaled(intLegendre, order);
+
+
+ // Basis //
+ basis = new Polynomial[size];
+
+ // Vertex Based (Lagrange) //
+ basis[0] =
+ Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
+
+ basis[1] =
+ Polynomial(1, 1, 0, 0);
+
+ basis[2] =
+ Polynomial(1, 0, 1, 0);
+
+
+ // Lagrange Sum //
+ for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
+ lagrangeSum[i] = basis[i] + basis[j];
+
+ // Lagrange Sub //
+ for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
+ lagrangeSub[i] = basis[j] - basis[i];
+
+
+ // Edge Based //
+ int i = 3;
+
+ for(int l = 1; l < order; l++){
+ for(int e = 0; e < 3; e++){
+ basis[i] =
+ intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]);
+
+ i++;
+ }
+ }
+
+ // Cell Based //
+ Polynomial p = basis[2] * 2 - Polynomial(1, 0, 0, 0);
+ const int orderMinusTwo = order - 2;
+
+ for(int l1 = 1; l1 < orderMinus; l1++){
+ for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
+ basis[i] =
+ intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) *
+ legendre[l2].compose(p) * basis[2];
+
+ i++;
+ }
+ }
+
+ // Free Temporary Sapce //
+ delete[] legendre;
+ delete[] intLegendre;
+ delete[] lagrangeSub;
+ delete[] lagrangeSum;
+}
+
+TriNodeBasis::~TriNodeBasis(void){
+ delete[] basis;
+}
+
+/*
+#include <cstdio>
+int main(void){
+ const int P = 5;
+ const double d = 0.01;
+
+ TriNodeBasis 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 - i + 1\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 = 0; i < b.getSize(); i++)
+ printf("figure;\ncontourf(x, y, p%d);\ncolorbar;\n", i + 1, i + 1);
+
+ printf("\n");
+
+ return 0;
+}
+*/
diff --git a/FunctionSpace/TriNodeBasis.h b/FunctionSpace/TriNodeBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f39e7969dbd80dcb948e8ccf338e627b65ac69d
--- /dev/null
+++ b/FunctionSpace/TriNodeBasis.h
@@ -0,0 +1,33 @@
+#ifndef _TRINODEBASIS_H_
+#define _TRINODEBASIS_H_
+
+#include "BasisScalar.h"
+
+/**
+ @class TriNodeBasis
+ @brief A Node-Basis for Triangles
+
+ This class can instantiate a Node-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 TriNodeBasis: public BasisScalar{
+ public:
+ TriNodeBasis(const int order);
+ virtual ~TriNodeBasis(void);
+};
+
+/**
+ @fn TriNodeBasis::TriNodeBasis(const int order)
+ @param order The order of the Basis
+ @return Returns a new Node-Basis for Triangles of the given order
+
+ @fn TriNodeBasis::~TriNodeBasis(void)
+ @return Deletes this Basis
+*/
+
+#endif
diff --git a/FunctionSpace/Vector.h b/FunctionSpace/Vector.h
new file mode 100644
index 0000000000000000000000000000000000000000..b808b73a74967b883ae78bf174f4a715d6ed8b43
--- /dev/null
+++ b/FunctionSpace/Vector.h
@@ -0,0 +1,235 @@
+#ifndef _VECTOR_H_
+#define _VECTOR_H_
+
+#include <string>
+
+extern "C"{
+#include <cblas.h>
+}
+
+#include "Exception.h"
+
+/**
+ @class Vector
+ @brief Handles vectors
+
+ This class represents a vector of type < T >
+ and of dimention @c n.
+
+ @todo
+ Use Inheritance instead of template@n
+ |@n
+ |@n
+ ---> We have methods that changes between types
+ (e.g: at(...), 3D automatic allocation, scalar mult)
+ */
+
+class Solver;
+
+template<class T>
+class Vector{
+ private:
+ int N;
+ T* v;
+
+ friend class Solver;
+ friend class Matrix;
+
+ public:
+ Vector(const int a);
+ Vector(void);
+ ~Vector(void);
+
+ int dim(void) const;
+
+ T& operator()(const int i);
+ T operator()(const int i) const;
+ void operator=(const Vector<T>& other);
+
+ T get(const int i) const;
+ void set(const int i, const T a);
+
+ Vector<double> at(const double x, const double y, const double z) const;
+
+ Vector<T> operator+(const Vector<T>& other);
+ Vector<T> operator-(const Vector<T>& other);
+ Vector<T> operator*(const T& other);
+ //Vector<T> operator*(const double alpha);
+
+ void add(const Vector<T>& b);
+ void sub(const Vector<T>& b);
+ void mul(const T& other);
+ //void mul(const double alpha);
+ double dot(const Vector<T>& v) const;
+
+ void allToZero(void);
+
+ std::string toString(void) const;
+};
+
+
+/**
+ @fn Vector::Vector(int)
+ @param a Size of the futur Vector
+ @return Returns a new Vector of size @c a
+
+ @fn Vector::Vector(void)
+ @return Returns a new Vector of size @em @c 3
+
+ @fn Vector::~Vector
+ @return Deletes this Vector
+
+ @fn Vector::dim
+ @return Returns the @em dimention of this vector
+
+ @fn T& Vector::operator()(const int)
+ @param i An index of this Vector
+ @return Returns a @em reference to
+ the element at position @c i
+
+ @fn T Vector::operator()(const int) const
+ @param i An index of this Vector
+ @return Returns the @em value of
+ the element at position @c i
+
+ @fn Vector::get
+ @param i An index in this Vector
+ @return Returns the @em value of the element
+ at position @c i
+
+ @fn Vector::set
+ @param i An index in this Vector
+ @param a A value
+ @return Sets the given value at position @c i
+ in this Vector
+
+ @fn Vector::at
+ @param x The @em first dimesion of a @c 3D Polynomial
+ @param y The @em second dimesion of a @c 3D Polynomial
+ @param z The @em third dimesion of a @c 3D Polynomial
+ @return Retuns the @c 3D Vector resulting of the @em evaluation
+ of this @em @c 3D Vector @em of @em Polynomial%s
+
+ @warning
+ This works only for @em @c 3D Vector @em of @em Polynomial%s
+
+ @fn Vector<T> Vector::operator+(const Vector<T>&)
+ @param other An other Vector
+ @return Returns the Vector containing
+ the @em sum of this Vector with the other
+
+ @fn Vector<T> Vector::operator-(const Vector<T>&)
+ @param other An other Vector
+ @return Returns the Vector containing
+ the @em substraction of this Vector with the other
+
+ @fn Vector<T> Vector::operator*(const T&)
+ @param other A value
+ @return Returns the Vector containing
+ the @em product (@em element @em by @em element)
+ of this Vector with the given value
+
+ @fn Vector::add
+ @param b An other Vector
+ @return Adds the given Vector with this one
+ @note
+ The result is stored in this Vector
+
+ @fn Vector::sub
+ @param b An other Vector
+ @return Substracts the given Vector with this one
+ @note
+ The result is stored in this Vector
+
+ @fn Vector::mul
+ @param other A value
+ @return Multiplies all the elements of this Vector
+ with the given value
+ @note
+ The result is stored in this Vector
+
+ @fn Vector::dot
+ @param v An Vector
+ @return Returns the @em dot @em product of this
+ Vector with the given one
+
+ @fn Vector::allToZero
+ @return Sets all the elements of this Vector
+ to the @em zero @em element
+
+ @fn Vector::toString
+ @return Returns a string representing this Vector
+*/
+
+//////////////////////////////////////////////////////////////////////
+// Templated Implementations //
+//////////////////////////////////////////////////////////////////////
+
+template<class T>
+Vector<T>::Vector(const int a){
+ if(!a)
+ throw Exception("Vector must by of dimension bigger than 0");
+
+ N = a;
+ v = new T[N];
+}
+
+template<class T>
+Vector<T>::Vector(void){
+ N = 3;
+ v = new T[N];
+}
+
+template<class T>
+Vector<T>::~Vector(void){
+ delete[] v;
+}
+
+//////////////////////////////////////////////////////////////////////
+// Inline Templated Implementations //
+//////////////////////////////////////////////////////////////////////
+
+template<class T>
+inline int Vector<T>::dim(void) const{
+ return N;
+}
+
+template<class T>
+inline T& Vector<T>::operator()(const int i){
+ return v[i];
+}
+
+template<class T>
+inline T Vector<T>::operator()(const int i) const{
+ return v[i];
+}
+
+template<class T>
+void Vector<T>::operator=(const Vector<T>& other){
+ if(N != other.N)
+ throw Exception("Vectors must be of the same dimension");
+
+ for(int i = 0; i < N; i++)
+ v[i] = other.v[i];
+}
+
+template<class T>
+inline T Vector<T>::get(const int i) const{
+ return v[i];
+}
+
+template<class T>
+inline void Vector<T>::set(const int i, const T a){
+ v[i] = a;
+}
+
+//////////////////////////////////////////////////////////////////////
+// Inline Vector<double> Implementations //
+//////////////////////////////////////////////////////////////////////
+
+template<>
+inline double Vector<double>::dot(const Vector<double>& v) const{
+ return cblas_ddot(N, (*this).v, 1, v.v, 1);
+}
+
+#endif
diff --git a/FunctionSpace/VectorDouble.cpp b/FunctionSpace/VectorDouble.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0c1dcbe98d975cb5a352d609b4ad3af0e47dfcd
--- /dev/null
+++ b/FunctionSpace/VectorDouble.cpp
@@ -0,0 +1,44 @@
+#include <sstream>
+#include "Vector.h"
+
+using namespace std;
+
+template<>
+void Vector<double>::allToZero(void){
+ for(int i = 0; i < N; i++)
+ v[i] = 0.0;
+}
+
+template<>
+void Vector<double>::add(const Vector<double>& b){
+ if(this->N != b.N)
+ throw Exception("Vectors must be with of the same size");
+
+ for(int i = 0; i < N; i++)
+ v[i] += b.v[i];
+}
+
+template<>
+void Vector<double>::sub(const Vector<double>& b){
+ if(this->N != b.N)
+ throw Exception("Vectors must be with of the same size");
+
+ for(int i = 0; i < N; i++)
+ v[i] -= b.v[i];
+}
+
+template<>
+void Vector<double>::mul(const double& alpha){
+ for(int i = 0; i < N; i++)
+ v[i] *= alpha;
+}
+
+template<>
+string Vector<double>::toString(void) const{
+ stringstream s;
+
+ for(int i = 0; i < N; i++)
+ s << scientific << showpos << v[i] << endl;
+
+ return s.str();
+}
diff --git a/FunctionSpace/VectorPolynomial.cpp b/FunctionSpace/VectorPolynomial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7cd83e08dc5e1b678737e20406763ba1715e3683
--- /dev/null
+++ b/FunctionSpace/VectorPolynomial.cpp
@@ -0,0 +1,53 @@
+#include <sstream>
+#include "Vector.h"
+#include "Polynomial.h"
+
+using namespace std;
+
+template<>
+Vector<double> Vector<Polynomial>::at(const double x,
+ const double y,
+ const double z) const{
+ Vector<double> val(N);
+
+ for(int i = 0; i < N; i++)
+ val(i) = v[i].at(x, y, z);
+
+ return val;
+}
+
+template<>
+void Vector<Polynomial>::add(const Vector<Polynomial>& other){
+ for(int i = 0; i < N; i++)
+ v[i].add(other.v[i]);
+}
+
+template<>
+void Vector<Polynomial>::sub(const Vector<Polynomial>& other){
+ for(int i = 0; i < N; i++)
+ v[i].sub(other.v[i]);
+}
+
+template<>
+void Vector<Polynomial>::mul(const Polynomial& other){
+ for(int i = 0; i < N; i++)
+ v[i].mul(other);
+}
+/*
+template<>
+void Vector<Polynomial>::mul(const double alpha){
+ for(int i = 0; i < N; i++)
+ v[i].mul(alpha);
+}
+*/
+template<>
+string Vector<Polynomial>::toString(void) const{
+ stringstream s;
+
+ s << endl;
+
+ for(int i = 0; i < N; i++)
+ s << "[" << v[i].toString() << "]" << endl;
+
+ return s.str();
+}