Preliminaries for New Interfaces of Vector Basis and TetEdgeBasis

FunctionSpace/BasisScalar.h

@@ -38,10 +38,17 @@ class BasisScalar: public Basis{
   //! 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;
+  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;
 

FunctionSpace/BasisVector.h

@@ -20,6 +20,11 @@
 
 class BasisVector: public Basis{
  protected:
+  std::vector            <const std::vector<Polynomial>*>*   node;
+  std::vector<std::vector<const std::vector<Polynomial>*>*>* edge;
+  std::vector<std::vector<const std::vector<Polynomial>*>*>* face;
+  std::vector            <const std::vector<Polynomial>*>*   cell;
+
   std::vector<const std::vector<Polynomial>*>* basis;
   std::vector<const std::vector<Polynomial>*>* revBasis;
 
@@ -34,6 +39,18 @@ class BasisVector: public Basis{
   const std::vector<const std::vector<Polynomial>*>& 
     getFunctions(unsigned int closure) const;
 
+  const std::vector<const std::vector<Polynomial>*>&
+    getNodeFunctions(void) const;
+  
+  const std::vector<std::vector<const std::vector<Polynomial>*>*>&
+    getEdgeFunctions(void) const;
+  
+  const std::vector<std::vector<const std::vector<Polynomial>*>*>&
+    getFaceFunctions(void) const;
+ 
+  const std::vector<const std::vector<Polynomial>*>&
+    getCellFunctions(void) const;
+
   virtual std::string toString(void) const;
 
  protected:
@@ -58,5 +75,28 @@ getFunctions(unsigned int closure) const{
     return *revBasis;
 }
 
+inline
+const std::vector<const std::vector<Polynomial>*>& 
+BasisVector::getNodeFunctions(void) const{
+  return *node;
+}
+
+inline  
+const std::vector<std::vector<const std::vector<Polynomial>*>*>& 
+BasisVector::getEdgeFunctions(void) const{
+  return *edge;
+}
+
+inline
+const std::vector<std::vector<const std::vector<Polynomial>*>*>& 
+BasisVector::getFaceFunctions(void) const{
+  return *face;
+}
+
+inline
+const std::vector<const std::vector<Polynomial>*>& 
+BasisVector::getCellFunctions(void) const{
+  return *cell;
+}
 
 #endif

FunctionSpace/CMakeLists.txt

@@ -20,6 +20,7 @@ set(SRC
   HexNodeBasis.cpp
   HexEdgeBasis.cpp
   TetNodeBasis.cpp
+  TetEdgeBasis.cpp
   
   FunctionSpace.cpp
   FunctionSpaceScalar.cpp

FunctionSpace/HexEdgeBasis.h

@@ -5,10 +5,10 @@
 
 /**
    @class HexEdgeBasis
-   @brief An Edge Basis for Hexahedrons
+   @brief An Edge Basis for Hexahedra
  
    This class can instantiate an Edge-Based Basis 
-   (high or low order) for Hexahedrons.@n
+   (high or low order) for Hexahedra.@n
    
    It uses 
    <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
@@ -19,7 +19,7 @@ class HexEdgeBasis: public BasisVector{
  public:
   //! @param order The order of the Basis
   //!
-  //! Returns a new Edge-Basis for Hexahedrons of the given order
+  //! Returns a new Edge-Basis for Hexahedra of the given order
   HexEdgeBasis(const int order);
   
   //! Deletes this Basis

FunctionSpace/HexNodeBasis.h

@@ -5,10 +5,10 @@
 
 /**
    @class HexNodeBasis
-   @brief A Node Basis for Hexahedrons
+   @brief A Node Basis for Hexahedra
  
    This class can instantiate a Node-Based Basis 
-   (high or low order) for Hexahedrons.@n
+   (high or low order) for Hexahedra.@n
    
    It uses 
    <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
@@ -19,7 +19,7 @@ class HexNodeBasis: public BasisScalar{
  public:
   //! @param order The order of the Basis
   //!
-  //! Returns a new Node-Basis for Hexahedrons of the given order
+  //! Returns a new Node-Basis for Hexahedra of the given order
   HexNodeBasis(const int order);
 
   //! @return Deletes this Basis

FunctionSpace/TetEdgeBasis.cpp

+#include "TetEdgeBasis.h"
+#include "Legendre.h"
+
+TetEdgeBasis::TetEdgeBasis(const int order){
+  // Set Basis Type //
+  this->order = order;
+ 
+  type = 1;
+  dim  = 2;
+
+  nVertex = 0;
+  nEdge   = 3 * (order + 1);
+  nFace   = 0;
+  nCell   = ((order - 1) * order + order - 1);
+
+  size = (order + 1) * (order + 2);
+
+  // Alloc Temporary Space //
+  const int orderPlus     = order + 1;
+  const int orderMinus    = order - 1;
+
+  Polynomial* legendre    = new Polynomial[orderPlus];
+  Polynomial* intLegendre = new Polynomial[orderPlus];
+  Polynomial* u           = new Polynomial[orderPlus];
+  Polynomial* v           = new Polynomial[orderPlus];
+
+  Polynomial  lagrange    [3];
+  Polynomial  lagrangeSum [3];
+  Polynomial  lagrangeSub [3];
+  Polynomial rLagrangeSub [3];
+
+  // Classical and Intrated-Scaled Legendre Polynomial //
+  Legendre::legendre(legendre, order);
+  Legendre::intScaled(intLegendre, orderPlus);
+
+  // Lagrange //
+  lagrange[0] = 
+    Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
+
+  lagrange[1] = 
+    Polynomial(1, 1, 0, 0);
+
+  lagrange[2] = 
+    Polynomial(1, 0, 1, 0);
+
+  // Lagrange Sum //
+  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
+     lagrangeSum[i] = lagrange[j] + lagrange[i];
+
+  // Lagrange Sub (& revert) //
+  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
+      lagrangeSub[i] = lagrange[i] - lagrange[j];
+     rLagrangeSub[i] = lagrange[j] - lagrange[i];
+  }
+
+
+  // Basis & Revert (temporary --- *no* const) //
+  std::vector<std::vector<Polynomial>*>    basis(size);
+  std::vector<std::vector<Polynomial>*> revBasis(size);
+
+
+  // Edge Based (Nedelec) //
+  int i = 0;
+
+  for(int j = 1; i < 3; j = (j + 1) % 3){
+    // Direct
+    std::vector<Polynomial> tmp = lagrange[j].gradient();
+    tmp[0].mul(lagrange[i]);
+    tmp[1].mul(lagrange[i]);
+    tmp[2].mul(lagrange[i]);
+
+    basis[i] = new std::vector<Polynomial>(lagrange[i].gradient());
+
+    basis[i]->at(0).mul(lagrange[j]);
+    basis[i]->at(1).mul(lagrange[j]);
+    basis[i]->at(2).mul(lagrange[j]);      
+
+    basis[i]->at(0).sub(tmp[0]);
+    basis[i]->at(1).sub(tmp[1]);
+    basis[i]->at(2).sub(tmp[2]);
+
+    // Revert 
+    std::vector<Polynomial> tmpR = lagrange[i].gradient();
+    tmpR[0].mul(lagrange[j]);
+    tmpR[1].mul(lagrange[j]);
+    tmpR[2].mul(lagrange[j]);
+
+    revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
+
+    revBasis[i]->at(0).mul(lagrange[i]);
+    revBasis[i]->at(1).mul(lagrange[i]);
+    revBasis[i]->at(2).mul(lagrange[i]);
+   
+    revBasis[i]->at(0).sub(tmpR[0]);
+    revBasis[i]->at(1).sub(tmpR[1]);
+    revBasis[i]->at(2).sub(tmpR[2]);
+
+    // Next
+    i++;
+  }
+
+  // Edge Based (High Order) //
+  for(int l = 1; l < orderPlus; l++){
+    for(int e = 0; e < 3; e++){
+      basis[i] = 
+	new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e],  
+							    lagrangeSum[e])).gradient());
+
+      revBasis[i] = 
+	new std::vector<Polynomial>((intLegendre[l].compose(rLagrangeSub[e], 
+							    lagrangeSum[e])).gradient());
+      i++;
+    }
+  }
+
+  // Cell Based (Preliminary) //
+  Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
+  
+  for(int l = 0; l < orderPlus; l++){
+    u[l] = intLegendre[l].compose(lagrangeSub[0] * -1, lagrangeSum[0]);
+    v[l] = legendre[l].compose(p);
+    v[l].mul(lagrange[2]);
+  }
+
+  // Cell Based (Type 1) //
+  for(int l1 = 1; l1 < order; l1++){
+    for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+      std::vector<Polynomial> tmp = v[l2].gradient();
+      tmp[0].mul(u[l1]);
+      tmp[1].mul(u[l1]);
+      tmp[2].mul(u[l1]);
+
+      basis[i] = new std::vector<Polynomial>(u[l1].gradient());
+      
+      basis[i]->at(0).mul(v[l2]);
+      basis[i]->at(1).mul(v[l2]);
+      basis[i]->at(2).mul(v[l2]);
+
+      basis[i]->at(0).add(tmp[0]);
+      basis[i]->at(1).add(tmp[1]);
+      basis[i]->at(2).add(tmp[2]);
+
+      revBasis[i] = basis[i];
+      
+      i++;
+    }
+  }
+
+  // Cell Based (Type 2) //
+  for(int l1 = 1; l1 < order; l1++){
+    for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+      std::vector<Polynomial> tmp = v[l2].gradient();
+      tmp[0].mul(u[l1]);
+      tmp[1].mul(u[l1]);
+      tmp[2].mul(u[l1]);
+
+      basis[i] = new std::vector<Polynomial>(u[l1].gradient());
+
+      basis[i]->at(0).mul(v[l2]);
+      basis[i]->at(1).mul(v[l2]);
+      basis[i]->at(2).mul(v[l2]);
+
+      basis[i]->at(0).sub(tmp[0]);
+      basis[i]->at(1).sub(tmp[1]);
+      basis[i]->at(2).sub(tmp[2]);
+ 
+      revBasis[i] = basis[i];
+
+      i++;
+    }
+  }
+
+  // Cell Based (Type 3) //
+  for(int l = 0; l < orderMinus; l++){
+    basis[i] = new std::vector<Polynomial>(*basis[0]);
+
+    basis[i]->at(0).mul(v[l]);
+    basis[i]->at(1).mul(v[l]);
+    basis[i]->at(2).mul(v[l]);
+
+    revBasis[i] = basis[i];
+    
+    i++;
+  }
+
+  // Free Temporary Sapce //
+  delete[] legendre;
+  delete[] intLegendre;
+  delete[] u;
+  delete[] v;
+
+
+  // Set Basis //
+  this->basis    = new std::vector<const std::vector<Polynomial>*>
+    (basis.begin(), basis.end());
+
+  this->revBasis = new std::vector<const std::vector<Polynomial>*>
+    (revBasis.begin(), revBasis.end());
+}
+
+TetEdgeBasis::~TetEdgeBasis(void){
+  for(int i = 0; i < size; i++){
+    delete (*basis)[i];
+
+    if(i >= nVertex && i < nVertex + nEdge)
+      delete (*revBasis)[i];
+  }
+
+  delete basis;
+  delete revBasis;
+}

FunctionSpace/TetEdgeBasis.h

+#ifndef _TETEDGEBASIS_H_
+#define _TETEDGEBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+   @class TetEdgeBasis
+   @brief An Edge Basis for Tetrahedra
+ 
+   This class can instantiate an Edge-Based Basis 
+   (high or low order) for Tetrahedra.@n
+   
+   It uses 
+   <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
+   Basis for @em high @em order Polynomial%s generation.@n
+*/
+
+class TetEdgeBasis: public BasisVector{
+ public:
+  //! @param order The order of the Basis
+  //!
+  //! Returns a new Edge-Basis for Tetrahedra of the given order
+  TetEdgeBasis(const int order);
+  
+  //! Deletes this Basis
+  //!
+  virtual ~TetEdgeBasis(void);
+};
+
+#endif

FunctionSpace/TetNodeBasis.h

@@ -5,10 +5,10 @@
 
 /**
    @class TetNodeBasis
-   @brief A Node Basis for Tetrahedrons
+   @brief A Node Basis for Tetrahedra
  
    This class can instantiate a Node-Based Basis 
-   (high or low order) for Tetrahedrons.@n
+   (high or low order) for Tetrahedra.@n
    
    It uses 
    <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
@@ -19,7 +19,7 @@ class TetNodeBasis: public BasisScalar{
  public:
   //! @param order The order of the Basis
   //!
-  //! Returns a new Node-Basis for Tetrahedrons of the given order
+  //! Returns a new Node-Basis for Tetrahedra of the given order
   TetNodeBasis(unsigned int order);
   
   //! Deletes this Basis