From 50e82020afbb26399220bcbe989ae55fe40dcd57 Mon Sep 17 00:00:00 2001 From: Nicolas Marsic <nicolas.marsic@gmail.com> Date: Fri, 19 Oct 2012 15:17:22 +0000 Subject: [PATCH] Preliminaries for New Interfaces of Vector Basis and TetEdgeBasis --- FunctionSpace/BasisScalar.h | 15 ++- FunctionSpace/BasisVector.h | 40 +++++++ FunctionSpace/CMakeLists.txt | 1 + FunctionSpace/HexEdgeBasis.h | 6 +- FunctionSpace/HexNodeBasis.h | 6 +- FunctionSpace/TetEdgeBasis.cpp | 211 +++++++++++++++++++++++++++++++++ FunctionSpace/TetEdgeBasis.h | 30 +++++ FunctionSpace/TetNodeBasis.h | 6 +- 8 files changed, 302 insertions(+), 13 deletions(-) create mode 100644 FunctionSpace/TetEdgeBasis.cpp create mode 100644 FunctionSpace/TetEdgeBasis.h diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h index 16b95fdd14..aafaef9e33 100644 --- a/FunctionSpace/BasisScalar.h +++ b/FunctionSpace/BasisScalar.h @@ -38,10 +38,17 @@ class BasisScalar: public Basis{ //! defining this (scalar) Basis, for the given closure const std::vector<const Polynomial*>& getFunctions(unsigned int closure) const; - const std::vector <const Polynomial*>& getNodeFunctions(void) const; - const std::vector<std::vector<const Polynomial*>*>& getEdgeFunctions(void) const; - const std::vector<std::vector<const Polynomial*>*>& getFaceFunctions(void) const; - const std::vector <const Polynomial*>& getCellFunctions(void) const; + const std::vector<const Polynomial*>& + getNodeFunctions(void) const; + + const std::vector<std::vector<const Polynomial*>*>& + getEdgeFunctions(void) const; + + const std::vector<std::vector<const Polynomial*>*>& + getFaceFunctions(void) const; + + const std::vector<const Polynomial*>& + getCellFunctions(void) const; virtual std::string toString(void) const; diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h index 66831266d4..6d74ec0dd6 100644 --- a/FunctionSpace/BasisVector.h +++ b/FunctionSpace/BasisVector.h @@ -20,6 +20,11 @@ class BasisVector: public Basis{ protected: + std::vector <const std::vector<Polynomial>*>* node; + std::vector<std::vector<const std::vector<Polynomial>*>*>* edge; + std::vector<std::vector<const std::vector<Polynomial>*>*>* face; + std::vector <const std::vector<Polynomial>*>* cell; + std::vector<const std::vector<Polynomial>*>* basis; std::vector<const std::vector<Polynomial>*>* revBasis; @@ -34,6 +39,18 @@ class BasisVector: public Basis{ const std::vector<const std::vector<Polynomial>*>& getFunctions(unsigned int closure) const; + const std::vector<const std::vector<Polynomial>*>& + getNodeFunctions(void) const; + + const std::vector<std::vector<const std::vector<Polynomial>*>*>& + getEdgeFunctions(void) const; + + const std::vector<std::vector<const std::vector<Polynomial>*>*>& + getFaceFunctions(void) const; + + const std::vector<const std::vector<Polynomial>*>& + getCellFunctions(void) const; + virtual std::string toString(void) const; protected: @@ -58,5 +75,28 @@ getFunctions(unsigned int closure) const{ return *revBasis; } +inline +const std::vector<const std::vector<Polynomial>*>& +BasisVector::getNodeFunctions(void) const{ + return *node; +} + +inline +const std::vector<std::vector<const std::vector<Polynomial>*>*>& +BasisVector::getEdgeFunctions(void) const{ + return *edge; +} + +inline +const std::vector<std::vector<const std::vector<Polynomial>*>*>& +BasisVector::getFaceFunctions(void) const{ + return *face; +} + +inline +const std::vector<const std::vector<Polynomial>*>& +BasisVector::getCellFunctions(void) const{ + return *cell; +} #endif diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt index a1c28aa15b..3ef4e3e3ac 100644 --- a/FunctionSpace/CMakeLists.txt +++ b/FunctionSpace/CMakeLists.txt @@ -20,6 +20,7 @@ set(SRC HexNodeBasis.cpp HexEdgeBasis.cpp TetNodeBasis.cpp + TetEdgeBasis.cpp FunctionSpace.cpp FunctionSpaceScalar.cpp diff --git a/FunctionSpace/HexEdgeBasis.h b/FunctionSpace/HexEdgeBasis.h index 29a2bc3ea2..6c3a6dc8c6 100644 --- a/FunctionSpace/HexEdgeBasis.h +++ b/FunctionSpace/HexEdgeBasis.h @@ -5,10 +5,10 @@ /** @class HexEdgeBasis - @brief An Edge Basis for Hexahedrons + @brief An Edge Basis for Hexahedra This class can instantiate an Edge-Based Basis - (high or low order) for Hexahedrons.@n + (high or low order) for Hexahedra.@n It uses <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a> @@ -19,7 +19,7 @@ class HexEdgeBasis: public BasisVector{ public: //! @param order The order of the Basis //! - //! Returns a new Edge-Basis for Hexahedrons of the given order + //! Returns a new Edge-Basis for Hexahedra of the given order HexEdgeBasis(const int order); //! Deletes this Basis diff --git a/FunctionSpace/HexNodeBasis.h b/FunctionSpace/HexNodeBasis.h index d6439e4b9b..e3d1603608 100644 --- a/FunctionSpace/HexNodeBasis.h +++ b/FunctionSpace/HexNodeBasis.h @@ -5,10 +5,10 @@ /** @class HexNodeBasis - @brief A Node Basis for Hexahedrons + @brief A Node Basis for Hexahedra This class can instantiate a Node-Based Basis - (high or low order) for Hexahedrons.@n + (high or low order) for Hexahedra.@n It uses <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a> @@ -19,7 +19,7 @@ class HexNodeBasis: public BasisScalar{ public: //! @param order The order of the Basis //! - //! Returns a new Node-Basis for Hexahedrons of the given order + //! Returns a new Node-Basis for Hexahedra of the given order HexNodeBasis(const int order); //! @return Deletes this Basis diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp new file mode 100644 index 0000000000..7c3aac36b2 --- /dev/null +++ b/FunctionSpace/TetEdgeBasis.cpp @@ -0,0 +1,211 @@ +#include "TetEdgeBasis.h" +#include "Legendre.h" + +TetEdgeBasis::TetEdgeBasis(const int order){ + // Set Basis Type // + this->order = order; + + type = 1; + dim = 2; + + nVertex = 0; + nEdge = 3 * (order + 1); + nFace = 0; + nCell = ((order - 1) * order + order - 1); + + size = (order + 1) * (order + 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* u = new Polynomial[orderPlus]; + Polynomial* v = new Polynomial[orderPlus]; + + Polynomial lagrange [3]; + Polynomial lagrangeSum [3]; + Polynomial lagrangeSub [3]; + Polynomial rLagrangeSub [3]; + + // 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 (& revert) // + for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){ + lagrangeSub[i] = lagrange[i] - lagrange[j]; + rLagrangeSub[i] = lagrange[j] - lagrange[i]; + } + + + // Basis & Revert (temporary --- *no* const) // + std::vector<std::vector<Polynomial>*> basis(size); + std::vector<std::vector<Polynomial>*> revBasis(size); + + + // Edge Based (Nedelec) // + int i = 0; + + for(int j = 1; i < 3; j = (j + 1) % 3){ + // Direct + std::vector<Polynomial> tmp = lagrange[j].gradient(); + tmp[0].mul(lagrange[i]); + tmp[1].mul(lagrange[i]); + tmp[2].mul(lagrange[i]); + + basis[i] = new std::vector<Polynomial>(lagrange[i].gradient()); + + basis[i]->at(0).mul(lagrange[j]); + basis[i]->at(1).mul(lagrange[j]); + basis[i]->at(2).mul(lagrange[j]); + + basis[i]->at(0).sub(tmp[0]); + basis[i]->at(1).sub(tmp[1]); + basis[i]->at(2).sub(tmp[2]); + + // Revert + std::vector<Polynomial> tmpR = lagrange[i].gradient(); + tmpR[0].mul(lagrange[j]); + tmpR[1].mul(lagrange[j]); + tmpR[2].mul(lagrange[j]); + + revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient()); + + revBasis[i]->at(0).mul(lagrange[i]); + revBasis[i]->at(1).mul(lagrange[i]); + revBasis[i]->at(2).mul(lagrange[i]); + + revBasis[i]->at(0).sub(tmpR[0]); + revBasis[i]->at(1).sub(tmpR[1]); + revBasis[i]->at(2).sub(tmpR[2]); + + // Next + i++; + } + + // Edge Based (High Order) // + for(int l = 1; l < orderPlus; l++){ + for(int e = 0; e < 3; e++){ + basis[i] = + new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e], + lagrangeSum[e])).gradient()); + + revBasis[i] = + new std::vector<Polynomial>((intLegendre[l].compose(rLagrangeSub[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++){ + std::vector<Polynomial> tmp = v[l2].gradient(); + tmp[0].mul(u[l1]); + tmp[1].mul(u[l1]); + tmp[2].mul(u[l1]); + + basis[i] = new std::vector<Polynomial>(u[l1].gradient()); + + basis[i]->at(0).mul(v[l2]); + basis[i]->at(1).mul(v[l2]); + basis[i]->at(2).mul(v[l2]); + + basis[i]->at(0).add(tmp[0]); + basis[i]->at(1).add(tmp[1]); + basis[i]->at(2).add(tmp[2]); + + revBasis[i] = basis[i]; + + i++; + } + } + + // Cell Based (Type 2) // + for(int l1 = 1; l1 < order; l1++){ + for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ + std::vector<Polynomial> tmp = v[l2].gradient(); + tmp[0].mul(u[l1]); + tmp[1].mul(u[l1]); + tmp[2].mul(u[l1]); + + basis[i] = new std::vector<Polynomial>(u[l1].gradient()); + + basis[i]->at(0).mul(v[l2]); + basis[i]->at(1).mul(v[l2]); + basis[i]->at(2).mul(v[l2]); + + basis[i]->at(0).sub(tmp[0]); + basis[i]->at(1).sub(tmp[1]); + basis[i]->at(2).sub(tmp[2]); + + revBasis[i] = basis[i]; + + i++; + } + } + + // Cell Based (Type 3) // + for(int l = 0; l < orderMinus; l++){ + basis[i] = new std::vector<Polynomial>(*basis[0]); + + basis[i]->at(0).mul(v[l]); + basis[i]->at(1).mul(v[l]); + basis[i]->at(2).mul(v[l]); + + revBasis[i] = basis[i]; + + i++; + } + + // Free Temporary Sapce // + delete[] legendre; + delete[] intLegendre; + delete[] u; + delete[] v; + + + // Set Basis // + this->basis = new std::vector<const std::vector<Polynomial>*> + (basis.begin(), basis.end()); + + this->revBasis = new std::vector<const std::vector<Polynomial>*> + (revBasis.begin(), revBasis.end()); +} + +TetEdgeBasis::~TetEdgeBasis(void){ + for(int i = 0; i < size; i++){ + delete (*basis)[i]; + + if(i >= nVertex && i < nVertex + nEdge) + delete (*revBasis)[i]; + } + + delete basis; + delete revBasis; +} diff --git a/FunctionSpace/TetEdgeBasis.h b/FunctionSpace/TetEdgeBasis.h new file mode 100644 index 0000000000..d8e0062f31 --- /dev/null +++ b/FunctionSpace/TetEdgeBasis.h @@ -0,0 +1,30 @@ +#ifndef _TETEDGEBASIS_H_ +#define _TETEDGEBASIS_H_ + +#include "BasisVector.h" + +/** + @class TetEdgeBasis + @brief An Edge Basis for Tetrahedra + + This class can instantiate an Edge-Based Basis + (high or low order) for Tetrahedra.@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 TetEdgeBasis: public BasisVector{ + public: + //! @param order The order of the Basis + //! + //! Returns a new Edge-Basis for Tetrahedra of the given order + TetEdgeBasis(const int order); + + //! Deletes this Basis + //! + virtual ~TetEdgeBasis(void); +}; + +#endif diff --git a/FunctionSpace/TetNodeBasis.h b/FunctionSpace/TetNodeBasis.h index 2083fddf59..9ade551271 100644 --- a/FunctionSpace/TetNodeBasis.h +++ b/FunctionSpace/TetNodeBasis.h @@ -5,10 +5,10 @@ /** @class TetNodeBasis - @brief A Node Basis for Tetrahedrons + @brief A Node Basis for Tetrahedra This class can instantiate a Node-Based Basis - (high or low order) for Tetrahedrons.@n + (high or low order) for Tetrahedra.@n It uses <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a> @@ -19,7 +19,7 @@ class TetNodeBasis: public BasisScalar{ public: //! @param order The order of the Basis //! - //! Returns a new Node-Basis for Tetrahedrons of the given order + //! Returns a new Node-Basis for Tetrahedra of the given order TetNodeBasis(unsigned int order); //! Deletes this Basis -- GitLab