diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 4c272d9988e2ebc93bc9b7cfc75ab6cfc614f4ee..76b0a5f48e89d2533431b27c0db16d6112a4d9bd 100644
--- a/FunctionSpace/Basis.h
+++ b/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;
}
diff --git a/FunctionSpace/BasisScalar.cpp b/FunctionSpace/BasisScalar.cpp
index 3d62551f6b5764496ac88b71b134d275d96f574c..1e792ed9365c7a2c5314cb0db87800cc5253bb2e 100644
--- a/FunctionSpace/BasisScalar.cpp
+++ b/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();
}
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index 3ef4e3e3ac5a78885976bbeae7886cbb3ae69c4e..19fe8c16a8c651a4ac0e56a9cca3bc8c7947a2ba 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/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
diff --git a/FunctionSpace/CurlBasis.cpp b/FunctionSpace/CurlBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc55ef79a18dbbdb8e215c171a71e73dcfac36c8
--- /dev/null
+++ b/FunctionSpace/CurlBasis.cpp
@@ -0,0 +1,96 @@
+#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;
+}
diff --git a/FunctionSpace/CurlBasis.h b/FunctionSpace/CurlBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..ffedacf3fe8dac79e6a6c7384a90d33155598011
--- /dev/null
+++ b/FunctionSpace/CurlBasis.h
@@ -0,0 +1,26 @@
+#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
diff --git a/FunctionSpace/DivBasis.cpp b/FunctionSpace/DivBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..15a4a387201b154df72fd2f2478e0c4d99681e93
--- /dev/null
+++ b/FunctionSpace/DivBasis.cpp
@@ -0,0 +1,96 @@
+#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;
+}
diff --git a/FunctionSpace/DivBasis.h b/FunctionSpace/DivBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ebf45e9bd99b2f6b604441904cd6e41708ce64f
--- /dev/null
+++ b/FunctionSpace/DivBasis.h
@@ -0,0 +1,27 @@
+#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
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 68fb4671196e69f69045a63fe62e3419ebd9103e..868209ecad7f281037bdabf15bd84c20f3413c1f 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/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);
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 4949652788f064fc64411acaeed51eebf826b26b..18f0cf387354a07a783b9aea50f2c53d16564da0 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/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);
};
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index cb41efbd271b93430e03cbf154e177d2b6ca15bc..10928c769c632edda0f450eaa29eaa0fa5e1769a 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/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){
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index dc0919a618ab06ddd3cf957b1ab3b7d3eac7162f..8480ee7543e27dae8b6d5e1fd61d8005b2d45ddc 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/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;
+}
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 7f1b536cb382db8f2b4e386919195484367fb168..8b0999e0b70746fa8924bde330f76829de4aab2a 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/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
diff --git a/FunctionSpace/GradBasis.cpp b/FunctionSpace/GradBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c4829e3ba5e307effe93b551d89d4ad5c354ec88
--- /dev/null
+++ b/FunctionSpace/GradBasis.cpp
@@ -0,0 +1,96 @@
+#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;
+}
diff --git a/FunctionSpace/GradBasis.h b/FunctionSpace/GradBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..34b32974c1bc336408136dcc7304d1fce5a42135
--- /dev/null
+++ b/FunctionSpace/GradBasis.h
@@ -0,0 +1,27 @@
+#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
diff --git a/FunctionSpace/HexEdgeBasis.cpp b/FunctionSpace/HexEdgeBasis.cpp
index f456994021c31ff7ef520335d332e0d6a5060758..7ebbbdee4101e8be34135ecf1fff54d4a65fb1b4 100644
--- a/FunctionSpace/HexEdgeBasis.cpp
+++ b/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 //
diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index cd3c99b97b34f61f0934bc57053ad4bd0262a1f4..5cc68e7f0cc0401270948fc6a9edcb035dee1375 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/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 //
diff --git a/FunctionSpace/Polynomial.cpp b/FunctionSpace/Polynomial.cpp
index 1edda763e19697aaccf4cebdedc9331709790d4a..fe952cc856924a8b60542b38924673aa2378b4b0 100644
--- a/FunctionSpace/Polynomial.cpp
+++ b/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{
diff --git a/FunctionSpace/Polynomial.h b/FunctionSpace/Polynomial.h
index fe5571e04bbd655987f775081750aedd87905b58..1e30fabe2e0cf1ce99a76cbf145b325befc78fa4 100644
--- a/FunctionSpace/Polynomial.h
+++ b/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;
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 3b76a66b6b4f1e45ed26b6617d226b2962171cea..56e6570e3baaf322d7851ed2f5f322943141e6d8 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/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 //
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 2d1e2e233d727ff0cfc57ad3eeca7860f11fbf24..68b0fc26a06f9e7f5e7fed6f261e1957a8e79a4e 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/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 //
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 705f36a783ef5463e4ce7d402a84975224ea8ef2..bc98182939468c6998fcecfdb1db3d31d10c5d21 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -15,6 +15,9 @@ TetEdgeBasis::TetEdgeBasis(int order){
nFace = 4 * (order + 1) * (order - 1);
nCell = (order + 1) * (order - 1) * (order - 2) / 2;
+ nEdgeClosure = 2;
+ nFaceClosure = 6;
+
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
@@ -343,7 +346,7 @@ TetEdgeBasis::TetEdgeBasis(int order){
}
TetEdgeBasis::~TetEdgeBasis(void){
- // Vertex Based //
+ // Vertex Based //
for(int i = 0; i < nVertex; i++)
delete (*node)[i];
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index c608aa794907fe86b784d7a2de386493fdcf0d57..05125b2d5d27a78df91d9e4da4a4de008e4cb744 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -15,6 +15,9 @@ TetNodeBasis::TetNodeBasis(int order){
nFace = 2 * (order - 1) * (order - 2);
nCell = (order - 1) * (order - 2) * (order - 3) / 6;
+ nEdgeClosure = 2;
+ nFaceClosure = 6;
+
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 332bf3fa067f370f01b43f4ad3f3aa14b0828661..eb02ba7dc87eb62c0da62a51a034944862543000 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -15,6 +15,9 @@ TriEdgeBasis::TriEdgeBasis(int order){
nFace = 0;
nCell = ((order - 1) * order + order - 1);
+ nEdgeClosure = 2;
+ nFaceClosure = 0;
+
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 929bdb43b9325736d497f3e6e6686825c9117154..061e6161fcd76edb33a9b13826a9a1e258b27017 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -14,6 +14,9 @@ TriNedelecBasis::TriNedelecBasis(void){
nFace = 0;
nCell = 0;
+ nEdgeClosure = 2;
+ nFaceClosure = 0;
+
size = 3;
// Vertices definig Edges & Permutations //
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 43a11e65808ddf478fb36f37f6b41d39cb362345..1cf10b7f7b29d6e38ad10b089c46d7e1a1e6a354 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -15,6 +15,9 @@ TriNodeBasis::TriNodeBasis(int order){
nFace = 0;
nCell = (order - 1) * (order - 2) / 2;
+ nEdgeClosure = 2;
+ nFaceClosure = 0;
+
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //