From 85b08109c38a36013e070582f7ce45213926cade Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 25 Oct 2012 13:26:25 +0000
Subject: [PATCH] Grad for Local Basis of a Function Space + Grad/Curl/Div of a
Basis
---
FunctionSpace/Basis.h | 19 +++++
FunctionSpace/BasisScalar.cpp | 17 ++++-
FunctionSpace/CMakeLists.txt | 4 +
FunctionSpace/CurlBasis.cpp | 96 ++++++++++++++++++++++++
FunctionSpace/CurlBasis.h | 26 +++++++
FunctionSpace/DivBasis.cpp | 96 ++++++++++++++++++++++++
FunctionSpace/DivBasis.h | 27 +++++++
FunctionSpace/FunctionSpace.cpp | 8 +-
FunctionSpace/FunctionSpace.h | 4 +-
FunctionSpace/FunctionSpaceNode.cpp | 4 +
FunctionSpace/FunctionSpaceScalar.cpp | 103 +++++++++++++++++++++++---
FunctionSpace/FunctionSpaceScalar.h | 19 ++++-
FunctionSpace/GradBasis.cpp | 96 ++++++++++++++++++++++++
FunctionSpace/GradBasis.h | 27 +++++++
FunctionSpace/HexEdgeBasis.cpp | 3 +
FunctionSpace/HexNodeBasis.cpp | 3 +
FunctionSpace/Polynomial.cpp | 41 ++++++++++
FunctionSpace/Polynomial.h | 4 +-
FunctionSpace/QuadEdgeBasis.cpp | 3 +
FunctionSpace/QuadNodeBasis.cpp | 3 +
FunctionSpace/TetEdgeBasis.cpp | 5 +-
FunctionSpace/TetNodeBasis.cpp | 3 +
FunctionSpace/TriEdgeBasis.cpp | 3 +
FunctionSpace/TriNedelecBasis.cpp | 3 +
FunctionSpace/TriNodeBasis.cpp | 3 +
25 files changed, 596 insertions(+), 24 deletions(-)
create mode 100644 FunctionSpace/CurlBasis.cpp
create mode 100644 FunctionSpace/CurlBasis.h
create mode 100644 FunctionSpace/DivBasis.cpp
create mode 100644 FunctionSpace/DivBasis.h
create mode 100644 FunctionSpace/GradBasis.cpp
create mode 100644 FunctionSpace/GradBasis.h
diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 4c272d9988..76b0a5f48e 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 3d62551f6b..1e792ed936 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 3ef4e3e3ac..19fe8c16a8 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 0000000000..fc55ef79a1
--- /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 0000000000..ffedacf3fe
--- /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 0000000000..15a4a38720
--- /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 0000000000..2ebf45e9bd
--- /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 68fb467119..868209ecad 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 4949652788..18f0cf3873 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 cb41efbd27..10928c769c 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 dc0919a618..8480ee7543 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 7f1b536cb3..8b0999e0b7 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 0000000000..c4829e3ba5
--- /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 0000000000..34b32974c1
--- /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 f456994021..7ebbbdee41 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 cd3c99b97b..5cc68e7f0c 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 1edda763e1..fe952cc856 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 fe5571e04b..1e30fabe2e 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 3b76a66b6b..56e6570e3b 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 2d1e2e233d..68b0fc26a0 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 705f36a783..bc98182939 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 c608aa7949..05125b2d5d 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 332bf3fa06..eb02ba7dc8 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 929bdb43b9..061e6161fc 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 43a11e6580..1cf10b7f7b 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 //
--
GitLab