diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp index 58376a0404fc47b421613e1e141482ea2f84a7fd..299201fe79acc3f79265d542670987249de19334 100644 --- a/FunctionSpace/BasisGenerator.cpp +++ b/FunctionSpace/BasisGenerator.cpp @@ -3,6 +3,8 @@ #include "Exception.h" #include "LineNodeBasis.h" +#include "LineEdgeBasis.h" +#include "LineNedelecBasis.h" #include "QuadNodeBasis.h" #include "QuadEdgeBasis.h" @@ -43,7 +45,10 @@ Basis* BasisGenerator::linGen(int basisType, int order){ switch(basisType){ case 0: return new LineNodeBasis(order); - case 1: throw Exception("1-form not implemented on Lines"); + case 1: + if (order == 0) return new LineNedelecBasis(); + else return new LineEdgeBasis(order); + case 2: throw Exception("2-form not implemented on Lines"); case 3: throw Exception("3-form not implemented on Lines"); diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt index 04d39e2f1e67a3653c37d2d4721d84f36ee629ae..d02168cc8d459e84a7d773d3d2448e91be037a6f 100644 --- a/FunctionSpace/CMakeLists.txt +++ b/FunctionSpace/CMakeLists.txt @@ -21,6 +21,8 @@ set(SRC TriLagrangeBasis.cpp LineNodeBasis.cpp + LineEdgeBasis.cpp + LineNedelecBasis.cpp QuadNodeBasis.cpp QuadEdgeBasis.cpp diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e88105404248ca2da641a1bc1e1546f8b9a39ba6 --- /dev/null +++ b/FunctionSpace/LineEdgeBasis.cpp @@ -0,0 +1,106 @@ +#include "LineEdgeBasis.h" +#include "Legendre.h" + +using namespace std; + +LineEdgeBasis::LineEdgeBasis(int order){ + // Set Basis Type // + this->order = order; + + type = 1; + dim = 1; + + nVertex = 0; + nEdge = (order + 1); + nFace = 0; + nCell = 0; + + nEdgeClosure = 2; + nFaceClosure = 0; + + size = nVertex + nEdge + nFace + nCell; + + // Alloc Temporary Space // + const int orderPlus = order + 1; + Polynomial* intLegendre = new Polynomial[orderPlus]; + + const Polynomial plusOneHalf(+0.5, 0, 0, 0); + const Polynomial minusOneHalf(-0.5, 0, 0, 0); + const Polynomial zero(0, 0, 0, 0); + + const Polynomial x[2] = { + Polynomial(+1, 1, 0, 0), + Polynomial(-1, 1, 0, 0) + }; + + // Legendre Polynomial // + Legendre::integrated(intLegendre, orderPlus); + + // Permutations // + const int permutation[2] = {0, 1}; + + // Basis // + node = new vector<vector<Polynomial>*>(nVertex); + edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure); + face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure); + cell = new vector<vector<Polynomial>*>(nCell); + + (*edge)[0] = new vector<vector<Polynomial>*>(nEdge); + (*edge)[1] = new vector<vector<Polynomial>*>(nEdge); + + + // Edge Based (Nedelec) // + (*(*edge)[0])[0] = new vector<Polynomial>(3); + (*(*edge)[0])[0]->at(0) = plusOneHalf; + (*(*edge)[0])[0]->at(1) = zero; + (*(*edge)[0])[0]->at(2) = zero; + + (*(*edge)[1])[0] = new vector<Polynomial>(3); + (*(*edge)[1])[0]->at(0) = minusOneHalf; + (*(*edge)[1])[0]->at(1) = zero; + (*(*edge)[1])[0]->at(2) = zero; + + + // Edge Based (High Order) // + for(int c = 0; c < 2; c++){ + unsigned int i = 0; + + for(int l = 1; l < orderPlus; l++){ + (*(*edge)[c])[i + 1] = + new vector<Polynomial>((intLegendre[l].compose(x[permutation[c]])).gradient()); + + i++; + } + } + + + // Free Temporary Space // + delete[] intLegendre; +} + +LineEdgeBasis::~LineEdgeBasis(void){ + // Vertex Based // + for(int i = 0; i < nVertex; i++) + delete (*node)[i]; + + delete node; + + + // Edge Based // + for(int c = 0; c < 2; c++){ + for(int i = 0; i < nEdge; i++) + delete (*(*edge)[c])[i]; + + delete (*edge)[c]; + } + + delete edge; + + + // Face Based // + delete face; + + + // Cell Based // + delete cell; +} diff --git a/FunctionSpace/LineEdgeBasis.h b/FunctionSpace/LineEdgeBasis.h new file mode 100644 index 0000000000000000000000000000000000000000..c5e7e794f46e1983fc334e01860103d2af55e9d8 --- /dev/null +++ b/FunctionSpace/LineEdgeBasis.h @@ -0,0 +1,34 @@ +#ifndef _LINEEDGEBASIS_H_ +#define _LINEEDGEBASIS_H_ + +#include "BasisVector.h" + +/** + @class LineEdgeBasis + @brief An Edge Basis for Lines + + This class can instantiate an Edge-Based Basis + (high or low order) for Lines.@n + + It uses an @em adaptation of + <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a> + Basis for @em high @em order Polynomial%s generation.@n + + This Basis is a restriction of a Quad Basis to @f$y = 0@f$.@n + + It also uses the following mapping: @f$x = \frac{u + 1}{2}@f$. +*/ + +class LineEdgeBasis: public BasisVector{ + public: + //! @param order The order of the Basis + //! + //! Returns a new Edge-Basis for Lines of the given order + LineEdgeBasis(int order); + + //! Deletes this Basis + //! + virtual ~LineEdgeBasis(void); +}; + +#endif diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp new file mode 100644 index 0000000000000000000000000000000000000000..04b9a79533136ad8cd90ab0edd19be080718dc3f --- /dev/null +++ b/FunctionSpace/LineNedelecBasis.cpp @@ -0,0 +1,75 @@ +#include "LineNedelecBasis.h" +#include "Legendre.h" + +using namespace std; + +LineNedelecBasis::LineNedelecBasis(void){ + // Set Basis Type // + order = 1; + + type = 1; + dim = 1; + + nVertex = 0; + nEdge = 1; + nFace = 0; + nCell = 0; + + nEdgeClosure = 2; + nFaceClosure = 0; + + size = 1; + + // Alloc Temporary Space // + const Polynomial plusOneHalf(+0.5, 0, 0, 0); + const Polynomial minusOneHalf(-0.5, 0, 0, 0); + const Polynomial zero(0, 0, 0, 0); + + // Basis // + node = new vector<vector<Polynomial>*>(nVertex); + edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure); + face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure); + cell = new vector<vector<Polynomial>*>(nCell); + + (*edge)[0] = new vector<vector<Polynomial>*>(nEdge); + (*edge)[1] = new vector<vector<Polynomial>*>(nEdge); + + + // Nedelec // + (*(*edge)[0])[0] = new vector<Polynomial>(3); + (*(*edge)[0])[0]->at(0) = plusOneHalf; + (*(*edge)[0])[0]->at(1) = zero; + (*(*edge)[0])[0]->at(2) = zero; + + (*(*edge)[1])[0] = new vector<Polynomial>(3); + (*(*edge)[1])[0]->at(0) = minusOneHalf; + (*(*edge)[1])[0]->at(1) = zero; + (*(*edge)[1])[0]->at(2) = zero; +} + +LineNedelecBasis::~LineNedelecBasis(void){ + // Vertex Based // + for(int i = 0; i < nVertex; i++) + delete (*node)[i]; + + delete node; + + + // Edge Based // + for(int c = 0; c < 2; c++){ + for(int i = 0; i < nEdge; i++) + delete (*(*edge)[c])[i]; + + delete (*edge)[c]; + } + + delete edge; + + + // Face Based // + delete face; + + + // Cell Based // + delete cell; +} diff --git a/FunctionSpace/LineNedelecBasis.h b/FunctionSpace/LineNedelecBasis.h new file mode 100644 index 0000000000000000000000000000000000000000..8579b1a47098ec87ffd39a00947478f935dc22c3 --- /dev/null +++ b/FunctionSpace/LineNedelecBasis.h @@ -0,0 +1,25 @@ +#ifndef _LINENEDELECBASIS_H_ +#define _LINENEDELECBASIS_H_ + +#include "BasisVector.h" + +/** + @class LineNedelecBasis + @brief Nedelec Basis for Lines + + This class can instantiate a Nedelec Basis + for Lines.@n +*/ + +class LineNedelecBasis: public BasisVector{ + public: + //! Returns a new Nedelec Basis for Lines + //! + LineNedelecBasis(void); + + //! Deletes this Basis + //! + virtual ~LineNedelecBasis(void); +}; + +#endif diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp index 762d68b31fb9ce4054d69b2352ceac5af706c9e1..4b94c8c5c2456e3810891619c3703ed05482eb7f 100644 --- a/FunctionSpace/LineNodeBasis.cpp +++ b/FunctionSpace/LineNodeBasis.cpp @@ -20,10 +20,20 @@ LineNodeBasis::LineNodeBasis(int order){ size = nVertex + nEdge + nFace + nCell; - // Legendre Polynomial // + // Alloc Temporary Space // Polynomial* intLegendre = new Polynomial[order]; + + const Polynomial x[2] = { + Polynomial(+1, 1, 0, 0), + Polynomial(-1, 1, 0, 0) + }; + + // Legendre Polynomial // Legendre::integrated(intLegendre, order); + // Permutations // + const int permutation[2] = {0, 1}; + // Basis // node = new vector<Polynomial*>(nVertex); edge = new vector<vector<Polynomial*>*>(nEdgeClosure); @@ -45,12 +55,6 @@ LineNodeBasis::LineNodeBasis(int order){ // Edge Based // - const int permutation[2] = {0, 1}; - const Polynomial x[2] = { - Polynomial(+1, 1, 0, 0), - Polynomial(-1, 1, 0, 0) - }; - for(int c = 0; c < 2; c++){ unsigned int i = 0;