Now Lagrange Basis (for Tri) completly in Basis Framework -- Still need a sexy...

Now Lagrange Basis (for Tri) completly in Basis Framework -- Still need a sexy way to choose between Zaglmayr and Lagrange -- Need to put LagrangeBasis::project elsewhere

FunctionSpace/BasisGenerator.cpp

@@ -12,6 +12,7 @@
 #include "TriNodeBasis.h"
 #include "TriEdgeBasis.h"
 #include "TriNedelecBasis.h"
+#include "TriLagrangeBasis.h"
 
 #include "TetNodeBasis.h"
 #include "TetEdgeBasis.h"
@@ -28,21 +29,57 @@ BasisGenerator::~BasisGenerator(void){
 
 Basis* BasisGenerator::generate(int elementType, 
 				int basisType, 
-				int order){
+				int order,
+				std::string family){
+  
+  if(!family.compare(std::string("zaglmayr")))
+    return generateZaglmayr(elementType, basisType, order);
+
+  else if(!family.compare(std::string("lagrange")))
+    return generateLagrange(elementType, basisType, order);
+
+  else
+    throw Exception("Unknwown Basis Family: %s", family.c_str());
+}
+
+Basis* BasisGenerator::generateZaglmayr(int elementType, 
+					int basisType, 
+					int order){
+  
+  switch(elementType){
+  case TYPE_LIN: return linZaglmayrGen(basisType, order);
+  case TYPE_TRI: return triZaglmayrGen(basisType, order);
+  case TYPE_QUA: return quaZaglmayrGen(basisType, order);
+  case TYPE_TET: return tetZaglmayrGen(basisType, order);
+  case TYPE_HEX: return hexZaglmayrGen(basisType, order);
+
+  default: throw Exception("Unknown Element Type (%d) for Basis Generation", 
+			   elementType);
+  }
+}
+
+Basis* BasisGenerator::generateLagrange(int elementType, 
+					int basisType, 
+					int order){
+  if(basisType != 0)
+    throw 
+      Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
+		basisType);
+
   switch(elementType){
-  case TYPE_LIN: return linGen(basisType, order);
-  case TYPE_TRI: return triGen(basisType, order);
-  case TYPE_QUA: return quaGen(basisType, order);
-  case TYPE_TET: return tetGen(basisType, order);
-  case TYPE_HEX: return hexGen(basisType, order);
+  case TYPE_LIN: throw Exception("Lagrange Basis on Lines not Implemented");
+  case TYPE_TRI: return new TriLagrangeBasis(order);
+  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 Basis Generation", 
 			   elementType);
   }
 }
 
-Basis* BasisGenerator::linGen(int basisType, 
-			      int order){
+Basis* BasisGenerator::linZaglmayrGen(int basisType, 
+				      int order){
   switch(basisType){ 
   case  0: return new LineNodeBasis(order);
   case  1: 
@@ -56,8 +93,8 @@ Basis* BasisGenerator::linGen(int basisType,
   }  
 }
 
-Basis* BasisGenerator::triGen(int basisType, 
-			      int order){
+Basis* BasisGenerator::triZaglmayrGen(int basisType, 
+				      int order){
   switch(basisType){
   case  0: return new TriNodeBasis(order);
   case  1: 
@@ -71,8 +108,8 @@ Basis* BasisGenerator::triGen(int basisType,
   }  
 }
 
-Basis* BasisGenerator::quaGen(int basisType, 
-			      int order){
+Basis* BasisGenerator::quaZaglmayrGen(int basisType, 
+				      int order){
   switch(basisType){
   case  0: return new QuadNodeBasis(order);
   case  1: return new QuadEdgeBasis(order);
@@ -83,8 +120,8 @@ Basis* BasisGenerator::quaGen(int basisType,
   }  
 }
 
-Basis* BasisGenerator::tetGen(int basisType, 
-			      int order){
+Basis* BasisGenerator::tetZaglmayrGen(int basisType, 
+				      int order){
   switch(basisType){
   case  0: return new TetNodeBasis(order);
   case  1: return new TetEdgeBasis(order);
@@ -95,8 +132,8 @@ Basis* BasisGenerator::tetGen(int basisType,
   }  
 }
 
-Basis* BasisGenerator::hexGen(int basisType, 
-			      int order){
+Basis* BasisGenerator::hexZaglmayrGen(int basisType, 
+				      int order){
   switch(basisType){
   case  0: return new HexNodeBasis(order);
   case  1: return new HexEdgeBasis(order);

FunctionSpace/BasisGenerator.h

 #ifndef _BASISGENERATOR_H_
 #define _BASISGENERATOR_H_
 
+#include <string>
 #include "Basis.h"
 
 /**
@@ -20,15 +21,29 @@ class BasisGenerator{
    BasisGenerator(void);
   ~BasisGenerator(void);
 
+  static Basis* generate(int elementType, 
+			 int basisType, 
+			 int order,
+			 std::string family);
+
   static Basis* generate(int elementType, 
 			 int basisType, 
 			 int order);
 
-  static Basis* linGen(int basisType, int order);
-  static Basis* triGen(int basisType, int order);
-  static Basis* quaGen(int basisType, int order);
-  static Basis* tetGen(int basisType, int order);
-  static Basis* hexGen(int basisType, int order);
+  static Basis* linZaglmayrGen(int basisType, int order);
+  static Basis* triZaglmayrGen(int basisType, int order);
+  static Basis* quaZaglmayrGen(int basisType, int order);
+  static Basis* tetZaglmayrGen(int basisType, int order);
+  static Basis* hexZaglmayrGen(int basisType, int order);
+
+ private:
+  static Basis* generateZaglmayr(int elementType, 
+				 int basisType, 
+				 int order);
+
+  static Basis* generateLagrange(int elementType, 
+				 int basisType, 
+				 int order);  
 };
 
 
@@ -50,8 +65,10 @@ class BasisGenerator{
    on which the requested Basis will be created
    @param basisType The Basis type
    @param order The order or the requested Basis
+   @param family A string
 
-   This method will @em instanciate the requested Basis
+   This method will @em instanciate the requested Basis,
+   of the requested family
 
    @return Returns a @em pointer to a newly 
    @em instantiated Basis
@@ -68,9 +85,28 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note Families are:
+   @li @c zaglmayr for 
+   <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>
+   Basis Functions
+   @li @c lagrange for Lagrange's Basis Functions   
+   **
+
+   @fn BasisGenerator::generate
+   @param elementType The type of the element,
+   on which the requested Basis will be created
+   @param basisType The Basis type
+   @param order The order or the requested Basis
+
+   Same as 
+   BasisGenerator::generate(@c elementType, @c basisType, @c order, @c zaglmayr)
+
+   @return Returns a @em pointer to a newly 
+   @em instantiated Basis
    **
 
-   @fn BasisGenerator::linGen
+   @fn BasisGenerator::linZaglmayrGen
    @param basisType The Basis type
    @param order The order or the requested Basis
 
@@ -85,9 +121,11 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note The Basis family will be @c zaglmayr
    **
 
-   @fn BasisGenerator::triGen
+   @fn BasisGenerator::triZaglmayrGen
    @param basisType The Basis type
    @param order The order or the requested Basis
 
@@ -102,9 +140,11 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note The Basis family will be @c zaglmayr
    **
 
-   @fn BasisGenerator::quaGen
+   @fn BasisGenerator::quaZaglmayrGen
    @param basisType The Basis type
    @param order The order or the requested Basis
 
@@ -119,9 +159,11 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note The Basis family will be @c zaglmayr
    **
 
-   @fn BasisGenerator::tetGen
+   @fn BasisGenerator::tetZaglmayrGen
    @param basisType The Basis type
    @param order The order or the requested Basis
 
@@ -136,9 +178,11 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note The Basis family will be @c zaglmayr
    **
 
-   @fn BasisGenerator::hexGen
+   @fn BasisGenerator::hexZaglmayrGen
    @param basisType The Basis type
    @param order The order or the requested Basis
 
@@ -153,7 +197,23 @@ class BasisGenerator{
    @li @c 1 for 1-Form
    @li @c 2 for 2-Form
    @li @c 3 for 3-Form
+
+   @note The Basis family will be @c zaglmayr
    **
  */
 
+//////////////////////
+// Inline Functions //
+//////////////////////
+
+inline Basis* BasisGenerator::generate(int elementType, 
+				       int basisType, 
+				       int order){
+  
+  return BasisGenerator::generate(elementType, 
+				  basisType, 
+				  order,
+				  "zaglmayr");
+}
+
 #endif

FunctionSpace/CMakeLists.txt

@@ -8,6 +8,7 @@ set(SRC
   Legendre.cpp
 
   Basis.cpp
+  LagrangeBasis.cpp
   BasisScalar.cpp
   BasisVector.cpp
   BasisGenerator.cpp
@@ -16,9 +17,6 @@ set(SRC
   CurlBasis.cpp
   DivBasis.cpp
 
-  LagrangeGenerator.cpp
-  LagrangeBasis.cpp
-  TriLagrangeBasis.cpp
 
   LineNodeBasis.cpp
   LineEdgeBasis.cpp
@@ -30,6 +28,7 @@ set(SRC
   TriNodeBasis.cpp
   TriEdgeBasis.cpp
   TriNedelecBasis.cpp
+  TriLagrangeBasis.cpp
 
   HexNodeBasis.cpp
   HexEdgeBasis.cpp

FunctionSpace/FunctionSpace.cpp

@@ -60,7 +60,7 @@ void FunctionSpace::build(const GroupOfElement& goe,
   type  = basisType;
   basis = BasisGenerator::generate(elementType, 
 				   basisType, 
-				   order);
+				   order, "lagrange");
 
   // Number of *Per* Entity functions //
   fPerVertex = basis->getNVertexBased() / nVertex;

FunctionSpace/LagrangeBasis.cpp

@@ -58,3 +58,42 @@ project(const MElement& element,
   // Return ;
   return newCoef;
 }
+
+unsigned int** LagrangeBasis::triEdgeOrder(int order){
+  // Check //
+  if(order <= 1)
+    return NULL;
+
+  // Number of Edge Based Basis //
+  const unsigned int orderMinus = order - 1;
+  const unsigned int size       = 3 * orderMinus;
+  
+  // Alloc //
+  unsigned int i;
+  unsigned int** edgeOrder = new unsigned int*[2];
+  edgeOrder[0]             = new unsigned int[size];
+  edgeOrder[1]             = new unsigned int[size];
+
+  // Direct Numeration //
+  i = 0;
+
+  for(unsigned int j = 0; j < orderMinus; j++){
+    for(unsigned int k = 0; k < size; k += orderMinus){
+      edgeOrder[0][i] = k + j;
+      i++;
+    }
+  }
+
+  // Invert Numeration //
+  i = 0;
+
+ for(unsigned int j = 0; j < orderMinus; j++){
+    for(unsigned int k = 0; k < size; k += orderMinus){
+      edgeOrder[1][i] = k + orderMinus - 1 - j;
+      i++;
+    }
+  }
+
+  // Retrun //
+  return edgeOrder;
+}

FunctionSpace/LagrangeBasis.h

@@ -67,6 +67,16 @@ class LagrangeBasis: public BasisScalar{
   //! Returns a new LagrangeBasis
   //!
   LagrangeBasis(void);
+
+  //! @param order The ordre of the requested LagrangeBasis
+  //! @return Returns a 2D table 
+  //! (of size 2 x (3 * (@c order - 1)))
+  //! with the Edge Triangle Inverted Closure for Lagrange
+  //! @note The first dimension of the returned table
+  //! is the direct closure, and the second dimension 
+  //! the invert closure
+  //! @warning The returned table @em must be @em deleted
+  static unsigned int** triEdgeOrder(int order);
 };
 
 

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 lagrangeOrder 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

@@ -66,27 +66,34 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
   type = 0;
   dim  = 2; 
 
-  nVertex = l->coefficients.size1();
-  nEdge   = 0;
+  nVertex = 3;
+  nEdge   = 3 * (order - 1);
   nFace   = 0;
-  nCell   = 0;
+  nCell   =     (order - 1) * (order - 2) / 2;
 
-  nEdgeClosure = 0;
+  nEdgeClosure = 2;
   nFaceClosure = 0;
 
-  size = nVertex;
- 
+  size = nVertex + nEdge + nFace + nCell;
+
+
+  // Alloc Some Stuff //
+  const int nMonomial = l->monomials.size1();
+  unsigned int** edgeOrder;
+  Polynomial* pEdgeClosureLess = new Polynomial[nEdge];
 
-  // Basis (pure nodal) //
+
+  // Basis //
   node = new vector<Polynomial*>(nVertex);
-  edge = new vector<vector<Polynomial*>*>(0);
+  edge = new vector<vector<Polynomial*>*>(2);
   face = new vector<vector<Polynomial*>*>(0);
-  cell = new vector<Polynomial*>(0);
+  cell = new vector<Polynomial*>(nCell);
 
+  (*edge)[0] = new vector<Polynomial*>(nEdge);
+  (*edge)[1] = new vector<Polynomial*>(nEdge);
 
-  // Instanciate Polynomials //
-  const int nMonomial = l->monomials.size1();
 
+  // Vertex Based //
   for(int i = 0; i < nVertex; i++){
     Polynomial p = Polynomial(0, 0, 0, 0);
     
@@ -98,6 +105,58 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
     
     (*node)[i] = new Polynomial(p);
   }
+
+
+  // Edge Based //
+  // Without Closure
+  for(int i = 0; i < nEdge; i++){
+    int ci              = i + nVertex;
+    pEdgeClosureLess[i] = Polynomial(0, 0, 0, 0);
+    
+    for(int j = 0; j < nMonomial; j++)
+      pEdgeClosureLess[i] = 
+	pEdgeClosureLess[i] + 
+	Polynomial(l->coefficients(ci, j), // Coef 
+		   l->monomials(j, 0),     // powerX
+		   l->monomials(j, 1),     // powerY
+		   0);                     // powerZ
+  }
+
+  // With Closure
+  edgeOrder = triEdgeOrder(order); // Closure Ordering
+
+  for(int i = 0; i < nEdge; i++){
+    (*(*edge)[0])[i] = 
+      new Polynomial(pEdgeClosureLess[edgeOrder[0][i]]);
+    
+    (*(*edge)[1])[i] = 
+      new Polynomial(pEdgeClosureLess[edgeOrder[1][i]]);
+  }
+
+
+  // Cell Based //
+  for(int i = 0; i < nCell; i++){
+    int ci       = i + nVertex + nEdge;
+    Polynomial p = Polynomial(0, 0, 0, 0);
+    
+    for(int j = 0; j < nMonomial; j++)
+      p = p + Polynomial(l->coefficients(ci, j), // Coef 
+			 l->monomials(j, 0),     // powerX
+			 l->monomials(j, 1),     // powerY
+			 0);                     // powerZ
+    
+    (*cell)[i] = new Polynomial(p);
+  }
+
+  
+  // Delete Temporary Space //
+  delete[] pEdgeClosureLess;
+
+  if(edgeOrder){
+    delete[] edgeOrder[0];
+    delete[] edgeOrder[1];
+    delete[] edgeOrder;
+  }
 }
 
 TriLagrangeBasis::~TriLagrangeBasis(void){