Curved Element: Poisson seems OK :-) --- Adaptive View for Scalar Field: OK :-)

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:

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
   

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{

FunctionSpace/LagrangeBasis.cpp

+#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;
+}

FunctionSpace/LagrangeBasis.h

+#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

FunctionSpace/LagrangeGenerator.cpp

+#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);
+  }
+}

FunctionSpace/LagrangeGenerator.h

+#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

FunctionSpace/TriLagrangeBasis.cpp

+#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;
+}

FunctionSpace/TriLagrangeBasis.h

+#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