diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h index 2f9ce3cf0a424899b54ee2233223e237d7afc7fc..08af96b40d6bb0ce74dee823f9eb93b2f567757c 100644 --- a/FunctionSpace/BasisGenerator.h +++ b/FunctionSpace/BasisGenerator.h @@ -58,6 +58,7 @@ class BasisGenerator{ @note Element types are: @li @c TYPE_TRI for Triangles @li @c TYPE_QUA for Quadrangles + @li @c TYPE_TET for Tetrahedrons @li @c TYPE_HEX for Hexahedrons @note Basis types are: diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt index 19fe8c16a8c651a4ac0e56a9cca3bc8c7947a2ba..18dab32b3c44e9429749548636fd6007b9ed53d8 100644 --- a/FunctionSpace/CMakeLists.txt +++ b/FunctionSpace/CMakeLists.txt @@ -18,11 +18,18 @@ set(SRC QuadNodeBasis.cpp QuadEdgeBasis.cpp + + LagrangeGenerator.cpp + LagrangeBasis.cpp + TriLagrangeBasis.cpp + TriNodeBasis.cpp TriEdgeBasis.cpp TriNedelecBasis.cpp + HexNodeBasis.cpp HexEdgeBasis.cpp + TetNodeBasis.cpp TetEdgeBasis.cpp diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h index 8ee705f4b119e8e96dc55da6787c33d08efbe684..712c4430e38f4962128e2c87bbc56d602225ec91 100644 --- a/FunctionSpace/FunctionSpace.h +++ b/FunctionSpace/FunctionSpace.h @@ -66,8 +66,9 @@ class FunctionSpace{ virtual ~FunctionSpace(void); const GroupOfElement& getSupport(void) const; - unsigned int getType(void) const; unsigned int getOrder(void) const; + unsigned int getType(void) const; + bool isScalar(void) const; unsigned int getNFunctionPerVertex(const MElement& element) const; unsigned int getNFunctionPerEdge(const MElement& element) const; @@ -132,6 +133,11 @@ class FunctionSpace{ FunctionSpace ** + @fn FunctionSpace::getOrder + @return Return the @em order + of this FunctionSpace + ** + @fn FunctionSpace::getType @return Return the @em type of the Basis functions composing @@ -139,9 +145,11 @@ class FunctionSpace{ @see Basis::getType() ** - @fn FunctionSpace::getOrder - @return Return the @em order - of this FunctionSpace + @fn FunctionSpace::isScalar + @return Return @c true if this + FunstionSpace is scalar, and + @c flase otherwise + @see Basis::isScalar() ** @fn FunctionSpace::getNFunctionPerVertex @@ -220,12 +228,16 @@ inline const GroupOfElement& FunctionSpace::getSupport(void) const{ return *goe; } +inline unsigned int FunctionSpace::getOrder(void) const{ + return (unsigned int)(basis->getOrder()); +} + inline unsigned int FunctionSpace::getType(void) const{ return type; } -inline unsigned int FunctionSpace::getOrder(void) const{ - return (unsigned int)(basis->getOrder()); +inline bool FunctionSpace::isScalar(void) const{ + return basis->isScalar(); } inline unsigned int FunctionSpace::getNFunctionPerVertex(const MElement& element) const{ diff --git a/FunctionSpace/LagrangeBasis.cpp b/FunctionSpace/LagrangeBasis.cpp new file mode 100644 index 0000000000000000000000000000000000000000..63b51e1efe49a18193412f0740bb508973948332 --- /dev/null +++ b/FunctionSpace/LagrangeBasis.cpp @@ -0,0 +1,34 @@ +#include "Exception.h" +#include "LagrangeBasis.h" + +using namespace std; + +LagrangeBasis::LagrangeBasis(void){ +} + +LagrangeBasis::~LagrangeBasis(void){ +} + +fullVector<double> LagrangeBasis:: +project(const fullVector<double>& coef, + const std::vector<const Polynomial*>& basis){ + + // Init New Coefs // + const unsigned int size = point->size1(); + fullVector<double> newCoef(size); + + // Interpolation at Lagrange Points // + const unsigned int nCoef = coef.size(); + + for(unsigned int i = 0; i < size; i++){ + newCoef(i) = 0; + + for(unsigned int j = 0; j < nCoef; j++) + newCoef(i) += coef(j) * basis[j]->at((*point)(i, 0), + (*point)(i, 1), + (*point)(i, 2)); + } + + // Return ; + return newCoef; +} diff --git a/FunctionSpace/LagrangeBasis.h b/FunctionSpace/LagrangeBasis.h new file mode 100644 index 0000000000000000000000000000000000000000..0d29cb6c782322d9792899c9e38d9c9755719bc2 --- /dev/null +++ b/FunctionSpace/LagrangeBasis.h @@ -0,0 +1,70 @@ +#ifndef _LAGRANGEBASIS_H_ +#define _LAGRANGEBASIS_H_ + +#include <vector> + +#include "polynomialBasis.h" +#include "BasisScalar.h" + +/** + @interface LagrangeBasis + @brief Interoface for Lagrange Basis + + This is an interface for Lagrange Basis.@n + + These Scalar Basis allow a @em Coefficient Matrix, + and a Monomial Matrix, to be consulted.@n + + A vector from an Other Basis (set of Functions) + can also be projected into a Lagrange Basis.@n + + @todo + Add a method to erase polynomialBasis in polynomialBasis@n + Add a method to get lagrange Point in polynomialBasis + */ + +class LagrangeBasis: public BasisScalar{ + protected: + fullMatrix<double>* point; + const polynomialBasis* l; + + public: + //! Deletes this Basis + //! + virtual ~LagrangeBasis(void); + + //! @return Returns the Coefficient Matrix + const fullMatrix<double>& getCoefficient(void) const; + + //! @return Returns the Monomial Matrix + const fullMatrix<double>& getMonomial(void) const; + + //! @param coef A vector of Real numbers + //! @param basis A vector of Polynomials + //! in a @em Reference Space + //! @return Projects the given Coefficients in this LagrangeBasis@n + fullVector<double> project(const fullVector<double>& coef, + const std::vector<const Polynomial*>& basis); + + protected: + //! Returns a new LagrangeBasis + //! + LagrangeBasis(void); +}; + + +////////////////////// +// Inline Functions // +////////////////////// + +inline const fullMatrix<double>& LagrangeBasis:: +getCoefficient(void) const{ + return l->coefficients; +} + +inline const fullMatrix<double>& LagrangeBasis:: +getMonomial(void) const{ + return l->monomials; +} + +#endif diff --git a/FunctionSpace/LagrangeGenerator.cpp b/FunctionSpace/LagrangeGenerator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..99c8f0c96070ad95b634168e661e3a34544db38d --- /dev/null +++ b/FunctionSpace/LagrangeGenerator.cpp @@ -0,0 +1,22 @@ +#include "Exception.h" +#include "TriLagrangeBasis.h" +#include "LagrangeGenerator.h" + +LagrangeGenerator::LagrangeGenerator(void){ +} + +LagrangeGenerator::~LagrangeGenerator(void){ +} + +LagrangeBasis* LagrangeGenerator::generate(unsigned int elementType, + unsigned int lagrangeOrder){ + switch(elementType){ + case TYPE_TRI: return new TriLagrangeBasis(lagrangeOrder); + case TYPE_QUA: throw Exception("Lagrange Basis on Quads not Implemented"); + case TYPE_TET: throw Exception("Lagrange Basis on Tets not Implemented"); + case TYPE_HEX: throw Exception("Lagrange Basis on Hexs not Implemented"); + + default: throw Exception("Unknown Element Type (%d) for Lagrange Basis", + elementType); + } +} diff --git a/FunctionSpace/LagrangeGenerator.h b/FunctionSpace/LagrangeGenerator.h new file mode 100644 index 0000000000000000000000000000000000000000..dea9a27fef673a89dad0c6a39d38ee180cb53a84 --- /dev/null +++ b/FunctionSpace/LagrangeGenerator.h @@ -0,0 +1,52 @@ +#ifndef _LAGRANGEGENERATOR_H_ +#define _LAGRANGEGENERATOR_H_ + +#include "LagrangeBasis.h" + +/** + @class LagrangeGenerator + @brief Generates a Lagrange Basis + + Generates a Lagrange Basis + */ + +class LagrangeGenerator{ + public: + LagrangeGenerator(void); + ~LagrangeGenerator(void); + + static LagrangeBasis* generate(unsigned int elementType, + unsigned int lagrangeOrder); +}; + +/** + @fn LagrangeGenerator::LagrangeGenerator + Instantiates a new LagrangeGenerator + + @note + A LagrangeGenerator got @em only @em class @em methods, + so it is not required to instanciate it. + ** + + @fn LagrangeGenerator::~LagrangeGenerator + Deletes this LagrangeGenerator + ** + + @fn LagrangeGenerator::generate + @param elementType The type of the element, + on which the requested LagrangeBasis will be created + @param order The order or the requested LagrangeBasis + + This method will @em instanciate the requested LagrangeBasis + + @return Returns a @em pointer to a newly + @em instantiated LagrangeBasis + + @note Element types are: + @li @c TYPE_TRI for Triangles + @li @c TYPE_QUA for Quadrangles + @li @c TYPE_TET for Tetrahedrons + @li @c TYPE_HEX for Hexahedrons +*/ + +#endif diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3c96b885ede5df82be23f3e23c42cce89b926b2b --- /dev/null +++ b/FunctionSpace/TriLagrangeBasis.cpp @@ -0,0 +1,211 @@ +#include "polynomialBasis.h" +#include "Exception.h" +#include "TriLagrangeBasis.h" + +using namespace std; + +TriLagrangeBasis::TriLagrangeBasis(int order){ + // Call polynomialBasis procedure // + int tag; + + switch(order){ + case 1: + tag = MSH_TRI_3; + break; + + case 2: + tag = MSH_TRI_6; + break; + + case 3: + tag = MSH_TRI_10; + break; + + case 4: + tag = MSH_TRI_15; + break; + + case 5: + tag = MSH_TRI_21; + break; + + case 6: + tag = MSH_TRI_28; + break; + + case 7: + tag = MSH_TRI_36; + break; + + case 8: + tag = MSH_TRI_45; + break; + + case 9: + tag = MSH_TRI_55; + break; + + case 10: + tag = MSH_TRI_66; + break; + + default: + throw Exception + ("Can't instanciate an order %d Lagrangian Basis for a Triangle", + order); + + break; + } + + l = polynomialBases::find(tag); + point = new fullMatrix<double>(triPoint(order)); + + // Set Basis Type // + this->order = order; + + type = 0; + dim = 2; + + nVertex = l->coefficients.size1(); + nEdge = 0; + nFace = 0; + nCell = 0; + + nEdgeClosure = 0; + nFaceClosure = 0; + + size = nVertex; + + + // Basis (pure nodal) // + node = new vector<Polynomial*>(nVertex); + edge = new vector<vector<Polynomial*>*>(0); + face = new vector<vector<Polynomial*>*>(0); + cell = new vector<Polynomial*>(0); + + + // Instanciate Polynomials // + const int nMonomial = l->monomials.size1(); + + for(int i = 0; i < nVertex; i++){ + Polynomial p = Polynomial(0, 0, 0, 0); + + for(int j = 0; j < nMonomial; j++) + p = p + Polynomial(l->coefficients(i, j), // Coef + l->monomials(j, 0), // powerX + l->monomials(j, 1), // powerY + 0); // powerZ + + (*node)[i] = new Polynomial(p); + } +} + +TriLagrangeBasis::~TriLagrangeBasis(void){ + // Delete gmsh polynomial Basis // + // --> no method to do so :-( + // --> erased ?? + + // Vertex Based // + for(int i = 0; i < nVertex; i++) + delete (*node)[i]; + + delete node; + + // Edge Based // + for(int c = 0; c < nEdgeClosure; c++){ + for(int i = 0; i < nEdge; i++) + delete (*(*edge)[c])[i]; + + delete (*edge)[c]; + } + + delete edge; + + // Face Based // + for(int c = 0; c < nFaceClosure; c++){ + for(int i = 0; i < nFace; i++) + delete (*(*face)[c])[i]; + + delete (*face)[c]; + } + + delete face; + + // Cell Based // + for(int i = 0; i < nCell; i++) + delete (*cell)[i]; + + delete cell; + + // Delete Lagrange Points // + delete point; +} + +fullMatrix<double> TriLagrangeBasis:: +triPoint(unsigned int order){ + unsigned int nbPoints = (order + 1) * (order + 2) / 2; + fullMatrix<double> point(nbPoints, 3); + + point(0, 0) = 0; + point(0, 1) = 0; + point(0, 2) = 0; + + if(order > 0){ + double dd = 1. / order; + + point(1, 0) = 1; + point(1, 1) = 0; + point(1, 2) = 0; + + point(2, 0) = 0; + point(2, 1) = 1; + point(2, 2) = 0; + + unsigned int index = 3; + + if(order > 1){ + double ksi = 0; + double eta = 0; + + for(unsigned int i = 0; i < order - 1; i++, index++){ + ksi += dd; + + point(index, 0) = ksi; + point(index, 1) = eta; + point(index, 2) = 0; + } + + ksi = 1.; + + for(unsigned int i = 0; i < order - 1; i++, index++){ + ksi -= dd; + eta += dd; + + point(index, 0) = ksi; + point(index, 1) = eta; + point(index, 2) = 0; + } + + eta = 1.; + ksi = 0.; + + for(unsigned int i = 0; i < order - 1; i++, index++){ + eta -= dd; + + point(index, 0) = ksi; + point(index, 1) = eta; + point(index, 2) = 0; + } + + if(order > 2){ + fullMatrix<double> inner = triPoint(order - 3); + + inner.scale(1. - 3. * dd); + inner.add(dd); + point.copy(inner, 0, nbPoints - index, 0, 2, index, 0); + } + } + } + + return point; +} diff --git a/FunctionSpace/TriLagrangeBasis.h b/FunctionSpace/TriLagrangeBasis.h new file mode 100644 index 0000000000000000000000000000000000000000..4ce5848bc80a44a047a199abf6c80eafe5021645 --- /dev/null +++ b/FunctionSpace/TriLagrangeBasis.h @@ -0,0 +1,39 @@ +#ifndef _TRILAGRANGEBASIS_H_ +#define _TRILAGRANGEBASIS_H_ + +#include "LagrangeBasis.h" + +/** + @class TriLagrangeBasis + @brief Lagrange Basis for Triangles + + This class can instantiate a @em Lagrange @em Basis + for a Triangle and for a given Order.@n + + It uses + <a href="http://geuz.org/gmsh/">gmsh</a> Basis.@n + + @todo + Add method to erase polynomialBasis in polynomialBasis + */ + +class TriLagrangeBasis: public LagrangeBasis{ + public: + //! @param odrer A natural number + //! + //! Returns a new for TriLagrangeBasis + //! of the given Order + TriLagrangeBasis(int order); + + //! Deletes this Basis + //! + virtual ~TriLagrangeBasis(void); + + private: + //! @param order A natural number + //! @return Returns Lagrangian Points on a Triangle + //! for the given Order + static fullMatrix<double> triPoint(unsigned int order); +}; + +#endif