Preliminary for 3D closure -- broken basis (except TriNodeBasis)

FunctionSpace/BasisScalar.h

@@ -20,6 +20,11 @@
 
 class BasisScalar: public Basis{
  protected:
+  std::vector            <const Polynomial*>*   node;
+  std::vector<std::vector<const Polynomial*>*>* edge;
+  std::vector<std::vector<const Polynomial*>*>* face;
+  std::vector            <const Polynomial*>*   cell;
+
   std::vector<const Polynomial*>* basis;
   std::vector<const Polynomial*>* revBasis;
 
@@ -28,15 +33,16 @@ class BasisScalar: public Basis{
   //!
   virtual ~BasisScalar(void);
   
-  //! @return Returns the set of @em Polynomial%s
-  //! defining this (scalar) Basis
-  const std::vector<const Polynomial*>& getFunctions(void) const;
-
   //! @param closure A natural number
   //! @return Returns the set of @em Polynomial%s
   //! 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;
+
   virtual std::string toString(void) const;
 
  protected:
@@ -61,5 +67,28 @@ getFunctions(unsigned int closure) const{
     return *revBasis;
 }
 
+inline
+const std::vector<const Polynomial*>& 
+BasisScalar::getNodeFunctions(void) const{
+  return *node;
+}
+
+inline  
+const std::vector<std::vector<const Polynomial*>*>& 
+BasisScalar::getEdgeFunctions(void) const{
+  return *edge;
+}
+
+inline
+const std::vector<std::vector<const Polynomial*>*>& 
+BasisScalar::getFaceFunctions(void) const{
+  return *face;
+}
+
+inline
+const std::vector<const Polynomial*>& 
+BasisScalar::getCellFunctions(void) const{
+  return *cell;
+}
 
 #endif

FunctionSpace/FunctionSpaceScalar.cpp

@@ -17,21 +17,42 @@ getLocalFunctions(const MElement& element) const{
     throw Exception("Element not found for closure");
 
   vector<bool>* closure   = it->second;
-  const unsigned int size = closure->size();
 
   // Get Basis //
-  const vector<const Polynomial*>&    basis = getBasis(element).getFunctions(0);
-  const vector<const Polynomial*>& revBasis = getBasis(element).getFunctions(1);
+  const vector       <const Polynomial*>&   node = getBasis(element).getNodeFunctions();
+  const vector<vector<const Polynomial*>*>& edge = getBasis(element).getEdgeFunctions();
+  const vector       <const Polynomial*>&   cell = getBasis(element).getCellFunctions();
 
   // Get Functions //
-  vector<const Polynomial*> fun(size);
+  unsigned int i           = 0;
+  const unsigned int nNode = node.size();
+  const unsigned int nEdge = edge[0]->size();
+  const unsigned int nCell = cell.size();
+
+  vector<const Polynomial*> fun(getBasis(element).getSize());
   
-  for(unsigned int i = 0; i < size; i++)
+  // Vertex Based //
+  for(unsigned int j = 0; j < nNode; j++){
+    fun[i] = node[j];
+    i++;
+  }
+
+  // Edge Based //
+  for(unsigned int j = 0; j < nEdge; j++){
     if((*closure)[i])
-      fun[i] = basis[i];
+      fun[i] = (*edge[0])[j];
 
     else
-      fun[i] = revBasis[i];
+      fun[i] = (*edge[1])[j];
+
+    i++;
+  }
+
+  // Vertex Based //
+  for(unsigned int j = 0; j < nCell; j++){
+    fun[i] = cell[j];
+    i++;
+  }
 
   // Return 
   return fun;

FunctionSpace/TriNodeBasis.cpp

 #include "TriNodeBasis.h"
 #include "Legendre.h"
 
+using namespace std;
+
 TriNodeBasis::TriNodeBasis(const int order){
   // Set Basis Type //
   this->order = order;
@@ -20,8 +22,7 @@ TriNodeBasis::TriNodeBasis(const int order){
   Polynomial* intLegendre  = new Polynomial[order];
 
   Polynomial  lagrangeSum[3];
-  Polynomial  lagrangeSub[3];
-  Polynomial rLagrangeSub[3];
+  Polynomial  lagrangeSub[2][3];
 
   // Classical and Intrated-Scaled Legendre Polynomial //
   const int orderMinus = order - 1;
@@ -30,83 +31,102 @@ TriNodeBasis::TriNodeBasis(const int order){
   Legendre::intScaled(intLegendre, order);
  
 
-  // Basis (& revert) //
-     basis = new std::vector<const Polynomial*>(size);
-  revBasis = new std::vector<const Polynomial*>(size);
+  // Basis //
+  node = new vector<const Polynomial*>(nVertex);
+  edge = new vector<vector<const Polynomial*>*>(2);
+  face = new vector<vector<const Polynomial*>*>(0);
+  cell = new vector<const Polynomial*>(nCell);
+
+  (*edge)[0] = new vector<const Polynomial*>(nEdge);
+  (*edge)[1] = new vector<const Polynomial*>(nEdge);
+
 
-  // Vertex Based (Lagrange) // 
-  (*basis)[0] = 
+  // Vertex Based (Lagrange) //
+  (*node)[0] = 
     new Polynomial(Polynomial(1, 0, 0, 0) - 
 		   Polynomial(1, 1, 0, 0) - 
 		   Polynomial(1, 0, 1, 0));
 
-  (*basis)[1] = 
+  (*node)[1] = 
     new Polynomial(Polynomial(1, 1, 0, 0));
 
-  (*basis)[2] = 
+  (*node)[2] = 
     new Polynomial(Polynomial(1, 0, 1, 0));
+ 
 
-  // Vertex Based (revert) //
-  for(int i = 0; i < 3; i++)
-    (*revBasis)[i] = (*basis)[i];  
-
-  // Lagrange Sum (& revert) //
+  // Lagrange Sum //
   for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
-     lagrangeSum[i] = *(*basis)[i] + *(*basis)[j];
+     lagrangeSum[i] = *(*node)[i] + *(*node)[j];
 
-  // Lagrange Sub (& revert) //
+  // Lagrange Sub //
   for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
-     lagrangeSub[i] = *(*basis)[j] - *(*basis)[i];
-    rLagrangeSub[i] = *(*basis)[i] - *(*basis)[j];
+    lagrangeSub[0][i] = *(*node)[j] - *(*node)[i];
+    lagrangeSub[1][i] = *(*node)[i] - *(*node)[j];
   }
 
-  
+ 
   // Edge Based //
-  int i = 3;
-
-  for(int l = 1; l < order; l++){
-    for(int e = 0; e < 3; e++){
-      (*basis)[i] = new Polynomial(
-	intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
-
-      (*revBasis)[i] = new Polynomial(
-	intLegendre[l].compose(rLagrangeSub[e], lagrangeSum[e]));
-      
-      i++;
+  for(int c = 0; c < 2; c++){
+    unsigned int i = 0;
+
+    for(int l = 1; l < order; l++){
+      for(int e = 0; e < 3; e++){
+	(*(*edge)[c])[i] = 
+	  new Polynomial(intLegendre[l].compose(lagrangeSub[c][e], lagrangeSum[e]));
+	
+	i++;
+      }
     }
   }
 
 
   // Cell Based //
-  Polynomial p             = *(*basis)[2] * 2 - Polynomial(1, 0, 0, 0);
+  Polynomial p             = *(*node)[2] * 2 - Polynomial(1, 0, 0, 0);
   const int  orderMinusTwo = order - 2;
+
+  unsigned int i = 0;
   
   for(int l1 = 1; l1 < orderMinus; l1++){
     for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
-      (*basis)[i] = new Polynomial(
-	intLegendre[l1].compose(lagrangeSub[0], lagrangeSum[0]) * 
-	legendre[l2].compose(p) * *(*basis)[2]);
-      
-      (*revBasis)[i] = (*basis)[i];
+      (*cell)[i] = 
+	new Polynomial(intLegendre[l1].compose(lagrangeSub[0][0], lagrangeSum[0]) * 
+		       legendre[l2].compose(p) * *(*node)[2]);
 
       i++;
     }
   }
-
-
+  
   // Free Temporary Sapce //
   delete[] legendre;
   delete[] intLegendre;
 }
 
 TriNodeBasis::~TriNodeBasis(void){
-  for(int i = 0; i < size; i++){
-    delete (*basis)[i];
+  // Vertex Based //
+  for(int i = 0; i < nVertex; i++)
+    delete (*node)[i];
+  
+  delete node;
+
+
+  // Edge Based //
+  for(int c = 0; c < 2; c++){
+    for(int i = 0; i < nEdge; i++)
+      delete (*(*edge)[c])[i];
     
-    if(i >= nVertex && i < nVertex + nEdge)
-      delete (*revBasis)[i];
+    delete (*edge)[c];
   }
+  
+  delete edge;
+
+
+  // Face Based //
+  delete face;
+
+
+  // Cell Based //
+  for(int i = 0; i < nCell; i++)
+    delete (*cell)[i];
 
-  delete basis;
-  delete revBasis;
+  delete cell;
 }