From ced129862c994a104f4d0d01c65447fb0c9ef503 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 6 Nov 2012 17:05:49 +0000
Subject: [PATCH] Curved Element: Poisson seems OK :-) --- Adaptive View for
Scalar Field: OK :-)
---
FunctionSpace/BasisGenerator.h | 1 +
FunctionSpace/CMakeLists.txt | 7 +
FunctionSpace/FunctionSpace.h | 24 +++-
FunctionSpace/LagrangeBasis.cpp | 34 +++++
FunctionSpace/LagrangeBasis.h | 70 +++++++++
FunctionSpace/LagrangeGenerator.cpp | 22 +++
FunctionSpace/LagrangeGenerator.h | 52 +++++++
FunctionSpace/TriLagrangeBasis.cpp | 211 ++++++++++++++++++++++++++++
FunctionSpace/TriLagrangeBasis.h | 39 +++++
9 files changed, 454 insertions(+), 6 deletions(-)
create mode 100644 FunctionSpace/LagrangeBasis.cpp
create mode 100644 FunctionSpace/LagrangeBasis.h
create mode 100644 FunctionSpace/LagrangeGenerator.cpp
create mode 100644 FunctionSpace/LagrangeGenerator.h
create mode 100644 FunctionSpace/TriLagrangeBasis.cpp
create mode 100644 FunctionSpace/TriLagrangeBasis.h
diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h
index 2f9ce3cf0a..08af96b40d 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 19fe8c16a8..18dab32b3c 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 8ee705f4b1..712c4430e3 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 0000000000..63b51e1efe
--- /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 0000000000..0d29cb6c78
--- /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 0000000000..99c8f0c960
--- /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 0000000000..dea9a27fef
--- /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 0000000000..3c96b885ed
--- /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 0000000000..4ce5848bc8
--- /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
--
GitLab