Grad for Local Basis of a Function Space + Grad/Curl/Div of a Basis

FunctionSpace/Basis.h

@@ -29,6 +29,9 @@ class Basis{
   int nFace;
   int nCell;
 
+  int nEdgeClosure;
+  int nFaceClosure;
+
   int size;
 
  public:
@@ -79,6 +82,14 @@ class Basis{
   //! @em Based functions of this Basis
   int getNCellBased(void) const;
 
+  //! @return Returns the number of @em edge
+  //! @em closures for this Basis
+  int getNEdgeClosure(void) const;
+
+  //! @return Returns the number of @em face
+  //! @em closures for this Basis
+  int getNFaceClosure(void) const;
+
   //! @return Returns the number of Polynomial%s 
   //! (or Vector%s of Polynomial%s) in the Basis
   int getSize(void) const;
@@ -130,6 +141,14 @@ inline int Basis::getNCellBased(void) const{
   return nCell;
 }
 
+inline int Basis::getNEdgeClosure(void) const{
+  return nEdgeClosure;
+}
+
+inline int Basis::getNFaceClosure(void) const{
+  return nFaceClosure;
+}
+
 inline int Basis::getSize(void) const{
   return size;
 }

FunctionSpace/BasisScalar.cpp

@@ -13,8 +13,21 @@ BasisScalar::~BasisScalar(void){
 string BasisScalar::toString(void) const{
   stringstream stream;
 
-  for(int i = 0; i < size; i++)
-    stream << "f(" << i << ") = " << endl;
+  for(int i = 0; i < nVertex; i++)
+    stream << "f(" << i << ") = " 
+	   << (*node)[i]->toString() << endl;
+
+  for(int i = 0; i < nEdge; i++)
+    stream << "f(" << i + nVertex << ") = " 
+	   << (*(*edge)[0])[i]->toString() << endl;
+
+  for(int i = 0; i < nFace; i++)
+    stream << "f(" << i + nVertex + nEdge << ") = " 
+	   << (*(*face)[0])[i]->toString() << endl;
+
+  for(int i = 0; i < nCell; i++)
+    stream << "f(" << i + nVertex + nEdge + nFace << ") = " 
+	   << (*cell)[i]->toString() << endl;
 
   return stream.str();
 }

FunctionSpace/CMakeLists.txt

@@ -12,6 +12,10 @@ set(SRC
   BasisVector.cpp
   BasisGenerator.cpp
 
+  GradBasis.cpp
+  CurlBasis.cpp
+  DivBasis.cpp
+
   QuadNodeBasis.cpp
   QuadEdgeBasis.cpp
   TriNodeBasis.cpp

FunctionSpace/CurlBasis.cpp

+#include "CurlBasis.h"
+
+using namespace std;
+
+CurlBasis::CurlBasis(const BasisVector& other){
+  // Set Basis Type //
+  order = other.getOrder();
+  
+  type = other.getType();
+  dim  = other.getDim();
+  
+  nVertex = other.getNVertexBased();
+  nEdge   = other.getNEdgeBased();
+  nFace   = other.getNFaceBased();
+  nCell   = other.getNCellBased();
+  
+  nEdgeClosure = other.getNEdgeClosure();
+  nFaceClosure = other.getNFaceClosure();
+
+  size = other.getSize();
+
+  // Alloc Basis //
+  node = new vector<vector<Polynomial>*>(nVertex);
+  edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure);
+  face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure);
+  cell = new vector<vector<Polynomial>*>(nCell);
+
+  for(int i = 0; i < nEdgeClosure; i++)
+    (*edge)[i] = new vector<vector<Polynomial>*>(nEdge);
+
+  for(int i = 0; i < nFaceClosure; i++)
+    (*face)[i] = new vector<vector<Polynomial>*>(nFace);  
+
+  // Node Based //
+  for(int i = 0; i < nVertex; i++)
+    (*node)[i] = 
+      new vector<Polynomial>
+      (Polynomial::curl(other.getNodeFunction(i)));
+  
+  // Edge Based //
+  for(int i = 0; i < nEdgeClosure; i++)
+    for(int j = 0; j < nEdge; j++)
+      (*(*edge)[i])[j] = 
+	new vector<Polynomial>
+	(Polynomial::curl(other.getEdgeFunction(i, j)));
+
+  // Face Based //
+  for(int i = 0; i < nFaceClosure; i++)
+    for(int j = 0; j < nFace; j++)
+      (*(*face)[i])[j] = 
+	new vector<Polynomial>
+	(Polynomial::curl(other.getFaceFunction(i, j)));
+
+  // Cell Based //
+  for(int i = 0; i < nCell; i++)
+    (*cell)[i] = 
+      new vector<Polynomial>
+      (Polynomial::curl(other.getCellFunction(i)));
+}
+
+CurlBasis::~CurlBasis(void){
+  // 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;
+}

FunctionSpace/CurlBasis.h

+#ifndef _CURLBASIS_H_
+#define _CURLBASIS_H_
+
+#include "BasisVector.h"
+
+/**
+   @class CurlBasis
+   @brief The Curl of an other Basis
+ 
+   This class can instantiate a Curl Basis.
+ */
+
+class CurlBasis: public BasisVector{
+ public:
+  //! @param other An Other Basis
+  //!
+  //! Returns a new Basis, which is the curl of
+  //! the given Basis
+  CurlBasis(const BasisVector& other);
+  
+  //! Deletes this Basis
+  //!
+  virtual ~CurlBasis(void);
+};
+
+#endif

FunctionSpace/DivBasis.cpp

+#include "DivBasis.h"
+
+using namespace std;
+
+DivBasis::DivBasis(const BasisVector& other){
+  // Set Basis Type //
+  order = other.getOrder();
+  
+  type = other.getType();
+  dim  = other.getDim();
+  
+  nVertex = other.getNVertexBased();
+  nEdge   = other.getNEdgeBased();
+  nFace   = other.getNFaceBased();
+  nCell   = other.getNCellBased();
+  
+  nEdgeClosure = other.getNEdgeClosure();
+  nFaceClosure = other.getNFaceClosure();
+
+  size = other.getSize();
+
+  // Alloc Basis //
+  node = new vector<Polynomial*>(nVertex);
+  edge = new vector<vector<Polynomial*>*>(2);
+  face = new vector<vector<Polynomial*>*>(6);
+  cell = new vector<Polynomial*>(nCell);
+
+  for(int i = 0; i < nEdgeClosure; i++)
+    (*edge)[i] = new vector<Polynomial*>(nEdge);
+
+  for(int i = 0; i < nFaceClosure; i++)
+    (*face)[i] = new vector<Polynomial*>(nFace);  
+
+  // Node Based //
+  for(int i = 0; i < nVertex; i++)
+    (*node)[i] = 
+      new Polynomial
+      (Polynomial::divergence(other.getNodeFunction(i)));
+  
+  // Edge Based //
+  for(int i = 0; i < nEdgeClosure; i++)
+    for(int j = 0; j < nEdge; j++)
+      (*(*edge)[i])[j] = 
+	new Polynomial
+	(Polynomial::divergence(other.getEdgeFunction(i, j)));
+
+  // Face Based //
+  for(int i = 0; i < nFaceClosure; i++)
+    for(int j = 0; j < nFace; j++)
+      (*(*face)[i])[j] = 
+	new Polynomial
+	(Polynomial::divergence(other.getFaceFunction(i, j)));
+
+  // Cell Based //
+  for(int i = 0; i < nCell; i++)
+    (*cell)[i] = 
+      new Polynomial
+      (Polynomial::divergence(other.getCellFunction(i)));
+}
+
+DivBasis::~DivBasis(void){
+  // 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;
+}

FunctionSpace/DivBasis.h

+#ifndef _DIVBASIS_H_
+#define _DIVBASIS_H_
+
+#include "BasisScalar.h"
+#include "BasisVector.h"
+
+/**
+   @class DivBasis
+   @brief The Divergence of an other Basis
+ 
+   This class can instantiate a Divergence Basis.
+ */
+
+class DivBasis: public BasisScalar{
+ public:
+  //! @param other An Other Basis
+  //!
+  //! Returns a new Basis, which is the divergence of
+  //! the given Basis
+  DivBasis(const BasisVector& other);
+  
+  //! Deletes this Basis
+  //!
+  virtual ~DivBasis(void);
+};
+
+#endif

FunctionSpace/FunctionSpace.cpp

@@ -61,7 +61,7 @@ void FunctionSpace::build(const GroupOfElement& goe,
   basis = BasisGenerator::generate(elementType, 
 				   basisType, 
 				   order);
-  
+
   // Number of *Per* Entity functions //
   fPerVertex = basis->getNVertexBased() / nVertex;
   // NB: fPreVertex = 0 *or* 1
@@ -99,7 +99,7 @@ void FunctionSpace::buildDof(void){
     // Get Dof for this Element
     vector<Dof> myDof = getKeys(*(element[i]));
     unsigned int nDof = myDof.size();
-    
+
     // Create new GroupOfDof
     GroupOfDof* god = new GroupOfDof(nDof, *(element[i])); 
     (*group)[i]     = god;
@@ -221,7 +221,7 @@ const GroupOfDof& FunctionSpace::getGoDFromElement(const MElement& element) cons
     return *(it->second); 
 }
 
-vector<int> FunctionSpace::getEdgeClosure(const MElement& element) const{
+vector<int> FunctionSpace::getEdgeClosure(const MElement& element){
   // Get not const Element (gmsh problem, not mine !) //
   MElement& eelement = const_cast<MElement&>(element);
 
@@ -246,7 +246,7 @@ vector<int> FunctionSpace::getEdgeClosure(const MElement& element) const{
   return edgeClosure;
 }
 
-vector<int> FunctionSpace::getFaceClosure(const MElement& element) const{
+vector<int> FunctionSpace::getFaceClosure(const MElement& element){
   // Get not const Element (gmsh problem, not mine !) //
   MElement& eelement = const_cast<MElement&>(element);
 

FunctionSpace/FunctionSpace.h

@@ -97,8 +97,8 @@ class FunctionSpace{
   void insertDof(Dof& d, GroupOfDof* god);    
 
   // Closure
-  std::vector<int> getEdgeClosure(const MElement& element) const;
-  std::vector<int> getFaceClosure(const MElement& element) const;
+  static std::vector<int> getEdgeClosure(const MElement& element);
+  static std::vector<int> getFaceClosure(const MElement& element);
 };
 
 

FunctionSpace/FunctionSpaceNode.cpp

@@ -10,6 +10,10 @@ FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe,
 				     int order){
   // Build 0Form Basis //
   build(goe, 0, order); 
+
+  // Init BasisScalar //
+  basisScalar = 
+    static_cast<const BasisScalar*>(basis);
 }
     
 FunctionSpaceNode::~FunctionSpaceNode(void){

FunctionSpace/FunctionSpaceScalar.cpp

@@ -3,31 +3,36 @@
 
 using namespace std;
 
+FunctionSpaceScalar::FunctionSpaceScalar(void){
+  hasGrad   = false;
+  gradBasis = NULL;
+}
+
 FunctionSpaceScalar::~FunctionSpaceScalar(void){
+  if(hasGrad)
+    delete gradBasis;
 }
 
 const vector<const Polynomial*> FunctionSpaceScalar::
 getLocalFunctions(const MElement& element) const{
   
   // Get Basis //
-  const BasisScalar&  basis = getBasis(element);
-
-  const unsigned int nFNode = basis.getNVertexBased();
-  const unsigned int nFEdge = basis.getNEdgeBased();
-  const unsigned int nFFace = basis.getNFaceBased();
-  const unsigned int nFCell = basis.getNCellBased();
+  const unsigned int nFNode = basisScalar->getNVertexBased();
+  const unsigned int nFEdge = basisScalar->getNEdgeBased();
+  const unsigned int nFFace = basisScalar->getNFaceBased();
+  const unsigned int nFCell = basisScalar->getNCellBased();
 
   // Get Closure //
   vector<int> edgeClosure = getEdgeClosure(element);
   vector<int> faceClosure = getFaceClosure(element); 
 
   // Get Functions //
-  vector<const Polynomial*> fun(basis.getSize());
+  vector<const Polynomial*> fun(basisScalar->getSize());
   unsigned int i = 0;
   
   // Vertex Based
   for(unsigned int j = 0; j < nFNode; j++){
-    fun[i] = &basis.getNodeFunction(j);
+    fun[i] = &basisScalar->getNodeFunction(j);
     i++;
   }
 
@@ -41,7 +46,7 @@ getLocalFunctions(const MElement& element) const{
   for(unsigned int j = 0; j < nFPerEdge; j++){
     for(unsigned int k = 0; k < nEdge; k++){
       fun[i] = 
-	&basis.getEdgeFunction(edgeClosure[k], fEdge);
+	&basisScalar->getEdgeFunction(edgeClosure[k], fEdge);
       
       fEdge++;
       i++;
@@ -58,7 +63,7 @@ getLocalFunctions(const MElement& element) const{
   for(unsigned int j = 0; j < nFPerFace; j++){
     for(unsigned int k = 0; k < nFace; k++){
       fun[i] = 
-	&basis.getFaceFunction(faceClosure[k], fFace);
+	&basisScalar->getFaceFunction(faceClosure[k], fFace);
       
       fFace++;
       i++;
@@ -67,7 +72,7 @@ getLocalFunctions(const MElement& element) const{
 
   // Cell Based
   for(unsigned int j = 0; j < nFCell; j++){
-    fun[i] = &basis.getCellFunction(j);
+    fun[i] = &basisScalar->getCellFunction(j);
     i++;
   }
 
@@ -75,3 +80,79 @@ getLocalFunctions(const MElement& element) const{
   // Return //
   return fun;
 }
+
+const vector<const vector<Polynomial>*> FunctionSpaceScalar::
+getGradLocalFunctions(const MElement& element) const{
+
+  // Got Grad Basis ? //
+  // --> mutable data 
+  //  --> Just a 'cache memory' 
+  if(!hasGrad){
+    gradBasis = new GradBasis(*basisScalar);
+    hasGrad   = true;
+  }
+
+  // Get Basis //
+  const unsigned int nFNode = gradBasis->getNVertexBased();
+  const unsigned int nFEdge = gradBasis->getNEdgeBased();
+  const unsigned int nFFace = gradBasis->getNFaceBased();
+  const unsigned int nFCell = gradBasis->getNCellBased();
+
+  // Get Closure //
+  vector<int> edgeClosure = getEdgeClosure(element);
+  vector<int> faceClosure = getFaceClosure(element); 
+
+  // Get Functions //
+  vector<const vector<Polynomial>*> fun(gradBasis->getSize());
+  unsigned int i = 0;
+  
+  // Vertex Based
+  for(unsigned int j = 0; j < nFNode; j++){
+    fun[i] = &gradBasis->getNodeFunction(j);
+    i++;
+  }
+
+  // Edge Based
+  // Number of basis function *per* edge
+  //  --> should always be an integer !
+  const unsigned int nEdge     = edgeClosure.size();
+  const unsigned int nFPerEdge = nFEdge / nEdge;
+  unsigned int fEdge = 0;
+
+  for(unsigned int j = 0; j < nFPerEdge; j++){
+    for(unsigned int k = 0; k < nEdge; k++){
+      fun[i] = 
+	&gradBasis->getEdgeFunction(edgeClosure[k], fEdge);
+      
+      fEdge++;
+      i++;
+    }
+  }
+
+  // Face Based
+  // Number of basis function *per* face
+  //  --> should always be an integer !
+  const unsigned int nFace     = faceClosure.size();
+  const unsigned int nFPerFace = nFFace / nFace;
+  unsigned int fFace = 0;
+
+  for(unsigned int j = 0; j < nFPerFace; j++){
+    for(unsigned int k = 0; k < nFace; k++){
+      fun[i] = 
+	&gradBasis->getFaceFunction(faceClosure[k], fFace);
+      
+      fFace++;
+      i++;
+    }
+  }
+
+  // Cell Based
+  for(unsigned int j = 0; j < nFCell; j++){
+    fun[i] = &gradBasis->getCellFunction(j);
+    i++;
+  }
+
+
+  // Return //
+  return fun;  
+}

FunctionSpace/FunctionSpaceScalar.h

@@ -3,6 +3,7 @@
 
 #include "fullMatrix.h"
 #include "BasisScalar.h"
+#include "GradBasis.h"
 #include "FunctionSpace.h"
 
 /**
@@ -21,6 +22,12 @@
 
 
 class FunctionSpaceScalar : public FunctionSpace{
+ protected:
+  const BasisScalar* basisScalar;
+
+  mutable bool       hasGrad;
+  mutable GradBasis* gradBasis;
+
  public:
   virtual ~FunctionSpaceScalar(void);
 
@@ -31,8 +38,14 @@ class FunctionSpaceScalar : public FunctionSpace{
   
   const std::vector<const Polynomial*> 
     getLocalFunctions(const MElement& element) const;
+
+  const std::vector<const std::vector<Polynomial>*>
+    getGradLocalFunctions(const MElement& element) const;
   
   const BasisScalar& getBasis(const MElement& element) const;
+
+ protected:
+  FunctionSpaceScalar(void);
 };
 
 
@@ -77,10 +90,10 @@ class FunctionSpaceScalar : public FunctionSpace{
 //////////////////////
 // Inline Functions //
 //////////////////////
-
+/*
 inline const BasisScalar& FunctionSpaceScalar::
 getBasis(const MElement& element) const{
-  return static_cast<const BasisScalar&>(*basis);
+  return *basisScalar;
 }
-
+*/
 #endif

FunctionSpace/GradBasis.cpp

+#include "GradBasis.h"
+
+using namespace std;
+
+GradBasis::GradBasis(const BasisScalar& other){
+  // Set Basis Type //
+  order = other.getOrder();
+  
+  type = other.getType();
+  dim  = other.getDim();
+  
+  nVertex = other.getNVertexBased();
+  nEdge   = other.getNEdgeBased();
+  nFace   = other.getNFaceBased();
+  nCell   = other.getNCellBased();
+  
+  nEdgeClosure = other.getNEdgeClosure();
+  nFaceClosure = other.getNFaceClosure();
+
+  size = other.getSize();
+
+  // Alloc Basis //
+  node = new vector<vector<Polynomial>*>(nVertex);
+  edge = new vector<vector<vector<Polynomial>*>*>(nEdgeClosure);
+  face = new vector<vector<vector<Polynomial>*>*>(nFaceClosure);
+  cell = new vector<vector<Polynomial>*>(nCell);
+
+  for(int i = 0; i < nEdgeClosure; i++)
+    (*edge)[i] = new vector<vector<Polynomial>*>(nEdge);
+
+  for(int i = 0; i < nFaceClosure; i++)
+    (*face)[i] = new vector<vector<Polynomial>*>(nFace);  
+
+  // Node Based //
+  for(int i = 0; i < nVertex; i++)
+    (*node)[i] = 
+      new vector<Polynomial>
+      (other.getNodeFunction(i).gradient());
+  
+  // Edge Based //
+  for(int i = 0; i < nEdgeClosure; i++)
+    for(int j = 0; j < nEdge; j++)
+      (*(*edge)[i])[j] = 
+	new vector<Polynomial>
+	(other.getEdgeFunction(i, j).gradient());
+
+  // Face Based //
+  for(int i = 0; i < nFaceClosure; i++)
+    for(int j = 0; j < nFace; j++)
+      (*(*face)[i])[j] = 
+	new vector<Polynomial>
+	(other.getFaceFunction(i, j).gradient());
+
+  // Cell Based //
+  for(int i = 0; i < nCell; i++)
+    (*cell)[i] = 
+      new vector<Polynomial>
+      (other.getCellFunction(i).gradient());
+}
+
+GradBasis::~GradBasis(void){
+  // 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;
+}

FunctionSpace/GradBasis.h

+#ifndef _GRADBASIS_H_
+#define _GRADBASIS_H_
+
+#include "BasisVector.h"
+#include "BasisScalar.h"
+
+/**
+   @class GradBasis
+   @brief The Gradient of an other Basis
+ 
+   This class can instantiate a Gradient Basis.
+ */
+
+class GradBasis: public BasisVector{
+ public:
+  //! @param other An Other Basis
+  //!
+  //! Returns a new Basis, which is the gradient of
+  //! the given Basis
+  GradBasis(const BasisScalar& other);
+  
+  //! Deletes this Basis
+  //!
+  virtual ~GradBasis(void);
+};
+
+#endif

FunctionSpace/HexEdgeBasis.cpp

@@ -17,6 +17,9 @@ HexEdgeBasis::HexEdgeBasis(int order){
   nFace   = 12         * order * (order + 1);
   nCell   =  3 * order * order * (order + 1);
 
+  nEdgeClosure = 2;
+  nFaceClosure = 8;
+
   size = 3 * (order + 2) * (order + 2) * (order + 1);
 
   // Alloc Temporary Space //

FunctionSpace/HexNodeBasis.cpp

@@ -16,6 +16,9 @@ HexNodeBasis::HexNodeBasis(int order){
   nFace   =  6 * (order - 1) * (order - 1);
   nCell   =      (order - 1) * (order - 1) * (order - 1);
 
+  nEdgeClosure = 2;
+  nFaceClosure = 8;
+
   size = nVertex + nEdge + nFace + nCell;
 
   // Alloc Temporary Space //

FunctionSpace/Polynomial.cpp

@@ -101,6 +101,47 @@ vector<Polynomial> Polynomial::gradient(void) const{
   return grad;
 }
 
+vector<Polynomial> Polynomial::curl(const vector<Polynomial>& p){
+  vector<Polynomial> rot(3);
+
+  // Partial Derivatives //
+  Polynomial dP0d1 = p[0];
+  Polynomial dP0d2 = p[0];
+  Polynomial dP1d0 = p[1];
+  Polynomial dP1d2 = p[1];
+  Polynomial dP2d0 = p[2];
+  Polynomial dP2d1 = p[2];
+
+  dP0d1.derivative(1);
+  dP0d2.derivative(2);
+  dP1d0.derivative(0);
+  dP1d2.derivative(2);
+  dP2d0.derivative(0);
+  dP2d1.derivative(1);
+
+  // Curl //
+  rot[0] = dP2d1 - dP1d2;
+  rot[1] = dP0d2 - dP2d0;
+  rot[2] = dP1d0 - dP0d1;
+
+  // Return //
+  return rot;
+}
+
+Polynomial Polynomial::divergence(const vector<Polynomial>& p){
+  // Partial Derivatives //
+  Polynomial dP0d0 = p[0];
+  Polynomial dP1d1 = p[1];
+  Polynomial dP2d2 = p[2];
+
+  dP0d0.derivative(0);
+  dP1d1.derivative(1);
+  dP2d2.derivative(2);
+
+  // Return Div //
+  return dP0d0 + dP1d1 + dP2d2;
+}
+
 double Polynomial::at
   (const double x, const double y, const double z) const{
   

FunctionSpace/Polynomial.h

@@ -36,8 +36,10 @@ class Polynomial{
    Polynomial(void);
   ~Polynomial(void);
 
-  void                    derivative(const int dim);
+  void derivative(const int dim);
   std::vector<Polynomial> gradient(void) const;
+  static std::vector<Polynomial> curl(const std::vector<Polynomial>& p);
+  static Polynomial divergence(const std::vector<Polynomial>& p);
 
   double operator()
     (const double x, const double y, const double z) const;  

FunctionSpace/QuadEdgeBasis.cpp

@@ -15,6 +15,9 @@ QuadEdgeBasis::QuadEdgeBasis(int order){
   nFace   = 0;
   nCell   = 2 * (order + 1) * order;
 
+  nEdgeClosure = 2;
+  nFaceClosure = 0;
+
   size = nVertex + nEdge + nFace + nCell;
 
   // Alloc Temporary Space //

FunctionSpace/QuadNodeBasis.cpp

@@ -15,6 +15,9 @@ QuadNodeBasis::QuadNodeBasis(int order){
   nFace   = 0;
   nCell   =     (order - 1) * (order - 1);
 
+  nEdgeClosure = 2;
+  nFaceClosure = 0;
+
   size = nVertex + nEdge + nFace + nCell;
 
   // Alloc Temporary Space //