From 8fc764bef603c7b0974f50ae8efcfc370c1a2507 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Tue, 23 Oct 2012 08:27:49 +0000
Subject: [PATCH] Tet Node Convergence: part2 -- Tet Edge: Part 1
---
FunctionSpace/BasisScalar.cpp | 2 +-
FunctionSpace/BasisScalar.h | 65 +++---
FunctionSpace/BasisVector.cpp | 6 +-
FunctionSpace/BasisVector.h | 64 ++----
FunctionSpace/FunctionSpace.cpp | 74 +++++++
FunctionSpace/FunctionSpace.h | 11 +-
FunctionSpace/FunctionSpaceScalar.cpp | 82 +------
FunctionSpace/FunctionSpaceVector.cpp | 80 +++++--
FunctionSpace/TetEdgeBasis.cpp | 301 ++++++++++++++++++--------
FunctionSpace/TetEdgeBasis.h | 2 +-
FunctionSpace/TetNodeBasis.cpp | 32 +--
FunctionSpace/TetNodeBasis.h | 2 +-
FunctionSpace/TriEdgeBasis.cpp | 194 +++++++++--------
FunctionSpace/TriEdgeBasis.h | 2 +-
FunctionSpace/TriNedelecBasis.cpp | 113 ++++++----
FunctionSpace/TriNodeBasis.cpp | 16 +-
FunctionSpace/TriNodeBasis.h | 2 +-
17 files changed, 610 insertions(+), 438 deletions(-)
diff --git a/FunctionSpace/BasisScalar.cpp b/FunctionSpace/BasisScalar.cpp
index 819ae6a419..3d62551f6b 100644
--- a/FunctionSpace/BasisScalar.cpp
+++ b/FunctionSpace/BasisScalar.cpp
@@ -14,7 +14,7 @@ string BasisScalar::toString(void) const{
stringstream stream;
for(int i = 0; i < size; i++)
- stream << "f(" << i << ") = " << (*basis)[i]->toString() << endl;
+ stream << "f(" << i << ") = " << endl;
return stream.str();
}
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
index aafaef9e33..450b5cc95a 100644
--- a/FunctionSpace/BasisScalar.h
+++ b/FunctionSpace/BasisScalar.h
@@ -20,10 +20,10 @@
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 <Polynomial*>* node;
+ std::vector<std::vector<Polynomial*>*>* edge;
+ std::vector<std::vector<Polynomial*>*>* face;
+ std::vector <Polynomial*>* cell;
std::vector<const Polynomial*>* basis;
std::vector<const Polynomial*>* revBasis;
@@ -33,22 +33,17 @@ class BasisScalar: public Basis{
//!
virtual ~BasisScalar(void);
- //! @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 Polynomial&
+ getNodeFunction(unsigned int i) const;
- const std::vector<std::vector<const Polynomial*>*>&
- getEdgeFunctions(void) const;
+ const Polynomial&
+ getEdgeFunction(unsigned int closure, unsigned int i) const;
- const std::vector<std::vector<const Polynomial*>*>&
- getFaceFunctions(void) const;
+ const Polynomial&
+ getFaceFunction(unsigned int closure, unsigned int i) const;
- const std::vector<const Polynomial*>&
- getCellFunctions(void) const;
+ const Polynomial&
+ getCellFunction(unsigned int i) const;
virtual std::string toString(void) const;
@@ -61,41 +56,31 @@ class BasisScalar: public Basis{
};
//////////////////////
-// Inline Functions //
+// Inline Function //
//////////////////////
inline
-const std::vector<const Polynomial*>& BasisScalar::
-getFunctions(unsigned int closure) const{
- if(!closure)
- return *basis;
-
- else
- return *revBasis;
-}
-
-inline
-const std::vector<const Polynomial*>&
-BasisScalar::getNodeFunctions(void) const{
- return *node;
+const Polynomial&
+BasisScalar::getNodeFunction(unsigned int i) const{
+ return *(*node)[i];
}
inline
-const std::vector<std::vector<const Polynomial*>*>&
-BasisScalar::getEdgeFunctions(void) const{
- return *edge;
+const Polynomial&
+BasisScalar::getEdgeFunction(unsigned int closure, unsigned int i) const{
+ return *(*(*edge)[closure])[i];
}
inline
-const std::vector<std::vector<const Polynomial*>*>&
-BasisScalar::getFaceFunctions(void) const{
- return *face;
+const Polynomial&
+BasisScalar::getFaceFunction(unsigned int closure, unsigned int i) const{
+ return *(*(*face)[closure])[i];
}
inline
-const std::vector<const Polynomial*>&
-BasisScalar::getCellFunctions(void) const{
- return *cell;
+const Polynomial&
+BasisScalar::getCellFunction(unsigned int i) const{
+ return *(*cell)[i];
}
#endif
diff --git a/FunctionSpace/BasisVector.cpp b/FunctionSpace/BasisVector.cpp
index 6eea126b47..b639e1d598 100644
--- a/FunctionSpace/BasisVector.cpp
+++ b/FunctionSpace/BasisVector.cpp
@@ -14,9 +14,9 @@ string BasisVector::toString(void) const{
stringstream stream;
for(int i = 0; i < size; i++)
- stream << "[" << (*basis)[i]->at(0).toString() << "]" << endl
- << "[" << (*basis)[i]->at(1).toString() << "]" << endl
- << "[" << (*basis)[i]->at(2).toString() << "]" << endl
+ stream << "[" << "]" << endl
+ << "[" << "]" << endl
+ << "[" << "]" << endl
<< endl;
return stream.str();
diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
index 6d74ec0dd6..5506d4cea7 100644
--- a/FunctionSpace/BasisVector.h
+++ b/FunctionSpace/BasisVector.h
@@ -20,10 +20,10 @@
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 <std::vector<Polynomial>*>* node;
+ std::vector<std::vector<std::vector<Polynomial>*>*>* edge;
+ std::vector<std::vector<std::vector<Polynomial>*>*>* face;
+ std::vector <std::vector<Polynomial>*>* cell;
std::vector<const std::vector<Polynomial>*>* basis;
std::vector<const std::vector<Polynomial>*>* revBasis;
@@ -33,23 +33,17 @@ class BasisVector: public Basis{
//!
virtual ~BasisVector(void);
- //! @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 std::vector<Polynomial>*>&
- getFunctions(unsigned int closure) const;
-
- const std::vector<const std::vector<Polynomial>*>&
- getNodeFunctions(void) const;
+ const std::vector<Polynomial>&
+ getNodeFunction(unsigned int i) const;
- const std::vector<std::vector<const std::vector<Polynomial>*>*>&
- getEdgeFunctions(void) const;
+ const std::vector<Polynomial>&
+ getEdgeFunction(unsigned int closure, unsigned int i) const;
- const std::vector<std::vector<const std::vector<Polynomial>*>*>&
- getFaceFunctions(void) const;
+ const std::vector<Polynomial>&
+ getFaceFunction(unsigned int closure, unsigned int i) const;
- const std::vector<const std::vector<Polynomial>*>&
- getCellFunctions(void) const;
+ const std::vector<Polynomial>&
+ getCellFunction(unsigned int i) const;
virtual std::string toString(void) const;
@@ -66,37 +60,27 @@ class BasisVector: public Basis{
//////////////////////
inline
-const std::vector<const std::vector<Polynomial>*>& BasisVector::
-getFunctions(unsigned int closure) const{
- if(!closure)
- return *basis;
-
- else
- return *revBasis;
-}
-
-inline
-const std::vector<const std::vector<Polynomial>*>&
-BasisVector::getNodeFunctions(void) const{
- return *node;
+const std::vector<Polynomial>&
+BasisVector::getNodeFunction(unsigned int i) const{
+ return *(*node)[i];
}
inline
-const std::vector<std::vector<const std::vector<Polynomial>*>*>&
-BasisVector::getEdgeFunctions(void) const{
- return *edge;
+const std::vector<Polynomial>&
+BasisVector::getEdgeFunction(unsigned int closure, unsigned int i) const{
+ return *(*(*edge)[closure])[i];
}
inline
-const std::vector<std::vector<const std::vector<Polynomial>*>*>&
-BasisVector::getFaceFunctions(void) const{
- return *face;
+const std::vector<Polynomial>&
+BasisVector::getFaceFunction(unsigned int closure, unsigned int i) const{
+ return *(*(*face)[closure])[i];
}
inline
-const std::vector<const std::vector<Polynomial>*>&
-BasisVector::getCellFunctions(void) const{
- return *cell;
+const std::vector<Polynomial>&
+BasisVector::getCellFunction(unsigned int i) const{
+ return *(*cell)[i];
}
#endif
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index 28647e5dba..68fb467119 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -221,6 +221,80 @@ const GroupOfDof& FunctionSpace::getGoDFromElement(const MElement& element) cons
return *(it->second);
}
+vector<int> FunctionSpace::getEdgeClosure(const MElement& element) const{
+ // Get not const Element (gmsh problem, not mine !) //
+ MElement& eelement = const_cast<MElement&>(element);
+
+ // Get Edge Closure //
+ const unsigned int nEdge = eelement.getNumEdges();
+ vector<int> edgeClosure(nEdge);
+
+ // Look Edges
+ for(unsigned int i = 0; i < nEdge; i++){
+ MEdge edge = eelement.getEdge(i);
+
+ if(edge.getVertex(0)->getNum() <
+ edge.getVertex(1)->getNum())
+
+ edgeClosure[i] = 0;
+
+ else
+ edgeClosure[i] = 1;
+ }
+
+ // Return //
+ return edgeClosure;
+}
+
+vector<int> FunctionSpace::getFaceClosure(const MElement& element) const{
+ // Get not const Element (gmsh problem, not mine !) //
+ MElement& eelement = const_cast<MElement&>(element);
+
+ // Get Face Closure //
+ const unsigned int nFace = eelement.getNumFaces();
+ vector<int> faceClosure(nFace);
+
+ // Look Faces
+ for(unsigned int i = 0; i < nFace; i++){
+ MFace face = eelement.getFace(i);
+
+ unsigned int v[3];
+ v[0] = face.getVertex(0)->getNum();
+ v[1] = face.getVertex(1)->getNum();
+ v[2] = face.getVertex(2)->getNum();
+
+ bool c[3];
+ c[0] = v[0] < v[1];
+ c[1] = v[1] < v[2];
+ c[2] = v[2] < v[0];
+
+ if(c[0])
+ if(c[1])
+ faceClosure[i] = 0;
+
+ else
+ if(c[2])
+ faceClosure[i] = 4;
+
+ else
+ faceClosure[i] = 5;
+
+ else
+ if(c[1])
+ if(c[2])
+ faceClosure[i] = 2;
+
+ else
+ faceClosure[i] = 1;
+
+ else
+ faceClosure[i] = 3;
+ }
+
+ // Return //
+ return faceClosure;
+}
+
string FunctionSpace::toString(void) const{
return basis->toString();
}
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 0ab0758505..4949652788 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -52,9 +52,6 @@ class FunctionSpace{
unsigned int fPerCell;
unsigned int type;
- // Closure //
- std::map<const MElement*, std::vector<bool>*>* edgeClosure;
-
// Dofs //
std::set<const Dof*, DofComparator>* dof;
std::vector<GroupOfDof*>* group;
@@ -89,15 +86,19 @@ class FunctionSpace{
std::string toString(void) const;
protected:
+ // Init
FunctionSpace(void);
void build(const GroupOfElement& goe,
int basisType, int order);
- void buildClosure(void);
+ // Dof
void buildDof(void);
+ void insertDof(Dof& d, GroupOfDof* god);
- void insertDof(Dof& d, GroupOfDof* god);
+ // Closure
+ std::vector<int> getEdgeClosure(const MElement& element) const;
+ std::vector<int> getFaceClosure(const MElement& element) const;
};
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index afd63a4c9e..dc0919a618 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -10,74 +10,16 @@ const vector<const Polynomial*> FunctionSpaceScalar::
getLocalFunctions(const MElement& element) const{
// Get Basis //
- const BasisScalar& basis = getBasis(element);
- const vector <const Polynomial*>& node = basis.getNodeFunctions();
- const vector<vector<const Polynomial*>*>& edge = basis.getEdgeFunctions();
- const vector<vector<const Polynomial*>*>& face = basis.getFaceFunctions();
- const vector <const Polynomial*>& cell = basis.getCellFunctions();
+ 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();
- // Get not const Element (gmsh problem, not mine !) //
- MElement& eelement = const_cast<MElement&>(element);
-
-
- // Get Edge Closure //
- const unsigned int nEdge = eelement.getNumEdges();
- vector<bool> edgeClosure(nEdge);
-
- // Look Edges
- for(unsigned int i = 0; i < nEdge; i++){
- MEdge edge = eelement.getEdge(i);
- edgeClosure[i] =
- edge.getVertex(0)->getNum() < edge.getVertex(1)->getNum();
- }
-
-
- // Get Face Closure //
- const unsigned int nFace = eelement.getNumFaces();
- vector<int> faceClosure(nFace);
-
- // Look Faces
- for(unsigned int i = 0; i < nFace; i++){
- MFace face = eelement.getFace(i);
-
- unsigned int v[3];
- v[0] = face.getVertex(0)->getNum();
- v[1] = face.getVertex(1)->getNum();
- v[2] = face.getVertex(2)->getNum();
-
- bool c[3];
- c[0] = v[0] < v[1];
- c[1] = v[1] < v[2];
- c[2] = v[2] < v[0];
-
- if(c[0])
- if(c[1])
- faceClosure[i] = 0;
-
- else
- if(c[2])
- faceClosure[i] = 4;
-
- else
- faceClosure[i] = 5;
-
- else
- if(c[1])
- if(c[2])
- faceClosure[i] = 2;
-
- else
- faceClosure[i] = 1;
-
- else
- faceClosure[i] = 3;
- }
-
+ // Get Closure //
+ vector<int> edgeClosure = getEdgeClosure(element);
+ vector<int> faceClosure = getFaceClosure(element);
// Get Functions //
vector<const Polynomial*> fun(basis.getSize());
@@ -85,23 +27,21 @@ getLocalFunctions(const MElement& element) const{
// Vertex Based
for(unsigned int j = 0; j < nFNode; j++){
- fun[i] = node[j];
+ fun[i] = &basis.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++){
- if(edgeClosure[k])
- fun[i] = (*edge[0])[fEdge];
-
- else
- fun[i] = (*edge[1])[fEdge];
+ fun[i] =
+ &basis.getEdgeFunction(edgeClosure[k], fEdge);
fEdge++;
i++;
@@ -111,12 +51,14 @@ getLocalFunctions(const MElement& element) const{
// 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] = (*face[faceClosure[k]])[fFace];
+ fun[i] =
+ &basis.getFaceFunction(faceClosure[k], fFace);
fFace++;
i++;
@@ -125,7 +67,7 @@ getLocalFunctions(const MElement& element) const{
// Cell Based
for(unsigned int j = 0; j < nFCell; j++){
- fun[i] = cell[j];
+ fun[i] = &basis.getCellFunction(j);
i++;
}
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index e34b0ba390..c326944054 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -6,36 +6,72 @@ using namespace std;
FunctionSpaceVector::~FunctionSpaceVector(void){
}
- const vector<const vector<Polynomial>*> FunctionSpaceVector::
+const vector<const vector<Polynomial>*> FunctionSpaceVector::
getLocalFunctions(const MElement& element) const{
- // Get Closure //
- map<const MElement*, vector<bool>*>::iterator it =
- edgeClosure->find(&element);
-
- if(it == edgeClosure->end())
- throw Exception("Element not found for closure");
+ // Get Basis //
+ const BasisVector& basis = getBasis(element);
- vector<bool>* closure = it->second;
- const unsigned int size = closure->size();
+ const unsigned int nFNode = basis.getNVertexBased();
+ const unsigned int nFEdge = basis.getNEdgeBased();
+ const unsigned int nFFace = basis.getNFaceBased();
+ const unsigned int nFCell = basis.getNCellBased();
- // Get Basis //
- const vector<const vector<Polynomial>*>& basis =
- getBasis(element).getFunctions(0);
-
- const vector<const vector<Polynomial>*>& revBasis =
- getBasis(element).getFunctions(1);
+ // Get Closure //
+ vector<int> edgeClosure = getEdgeClosure(element);
+ vector<int> faceClosure = getFaceClosure(element);
// Get Functions //
- vector<const vector<Polynomial>*> fun(size);
+ vector<const vector<Polynomial>*> fun(basis.getSize());
+ unsigned int i = 0;
- for(unsigned int i = 0; i < size; i++)
- if((*closure)[i])
- fun[i] = basis[i];
+ // Vertex Based
+ for(unsigned int j = 0; j < nFNode; j++){
+ fun[i] = &basis.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] =
+ &basis.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] =
+ &basis.getFaceFunction(faceClosure[k], fFace);
+
+ fFace++;
+ i++;
+ }
+ }
+
+ // Cell Based
+ for(unsigned int j = 0; j < nFCell; j++){
+ fun[i] = &basis.getCellFunction(j);
+ i++;
+ }
- else
- fun[i] = revBasis[i];
- // Return
+ // Return //
return fun;
}
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 7c3aac36b2..a155460801 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -1,38 +1,76 @@
#include "TetEdgeBasis.h"
#include "Legendre.h"
-TetEdgeBasis::TetEdgeBasis(const int order){
+using namespace std;
+
+TetEdgeBasis::TetEdgeBasis(int order){
// Set Basis Type //
this->order = order;
type = 1;
- dim = 2;
+ dim = 3;
nVertex = 0;
- nEdge = 3 * (order + 1);
- nFace = 0;
- nCell = ((order - 1) * order + order - 1);
+ nEdge = 6 * (order + 1);
+ nFace = 4 * (order + 1) * (order - 1);
+ nCell = (order + 1) * (order - 1) * (order - 2) / 2;
- size = (order + 1) * (order + 2);
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
const int orderMinus = order - 1;
+ const int orderMinusTwo = order - 2;
- Polynomial* legendre = new Polynomial[orderPlus];
- Polynomial* intLegendre = new Polynomial[orderPlus];
- Polynomial* u = new Polynomial[orderPlus];
- Polynomial* v = new Polynomial[orderPlus];
+ Polynomial lagrange[4];
- Polynomial lagrange [3];
- Polynomial lagrangeSum [3];
- Polynomial lagrangeSub [3];
- Polynomial rLagrangeSub [3];
+ Polynomial* legendre = new Polynomial[orderMinus];
+ Polynomial* sclLegendre = new Polynomial[orderMinus];
+ Polynomial* intLegendre = new Polynomial[orderPlus];
// Classical and Intrated-Scaled Legendre Polynomial //
- Legendre::legendre(legendre, order);
+ Legendre::legendre(legendre, orderMinusTwo);
+ Legendre::scaled(sclLegendre, orderMinusTwo);
Legendre::intScaled(intLegendre, orderPlus);
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][6][2] =
+ {
+ { {0, 1}, {1, 2}, {2, 0},
+ {3, 0}, {3, 2}, {3, 1} },
+
+ { {1, 0}, {2, 1}, {0, 2},
+ {0, 3}, {2, 3}, {1, 3} }
+ };
+
+ // Vertices definig Faces & Permutations //
+ const int faceV[6][4][3] =
+ {
+ { {0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2} },
+ { {2, 0, 1}, {1, 0, 3}, {3, 0, 2}, {1, 3, 2} },
+ { {2, 1, 0}, {1, 3, 0}, {3, 2, 0}, {1, 2, 3} },
+ { {1, 2, 0}, {3, 1, 0}, {2, 3, 0}, {2, 1, 3} },
+ { {1, 0, 2}, {3, 0, 1}, {2, 0, 3}, {2, 3, 1} },
+ { {0, 1, 2}, {0, 3, 1}, {0, 2, 3}, {3, 2, 1} },
+ };
+
+ // Basis //
+ node = new vector<vector<Polynomial>*>(nVertex);
+ edge = new vector<vector<vector<Polynomial>*>*>(2);
+ face = new vector<vector<vector<Polynomial>*>*>(6);
+ cell = new vector<vector<Polynomial>*>(nCell);
+
+ (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
+ (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+
+ (*face)[0] = new vector<vector<Polynomial>*>(nFace);
+ (*face)[1] = new vector<vector<Polynomial>*>(nFace);
+ (*face)[2] = new vector<vector<Polynomial>*>(nFace);
+ (*face)[3] = new vector<vector<Polynomial>*>(nFace);
+ (*face)[4] = new vector<vector<Polynomial>*>(nFace);
+ (*face)[5] = new vector<vector<Polynomial>*>(nFace);
+
+
// Lagrange //
lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
@@ -43,76 +81,135 @@ TetEdgeBasis::TetEdgeBasis(const int order){
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);
-
+ lagrange[3] =
+ Polynomial(1, 0, 0, 1);
+
// 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++;
+ for(int c = 0; c < 2; c++){
+ for(int e = 0; e < 6; e++){
+ vector<Polynomial> tmp = lagrange[edgeV[c][e][1]].gradient();
+
+ tmp[0].mul(lagrange[edgeV[c][e][0]]);
+ tmp[1].mul(lagrange[edgeV[c][e][0]]);
+ tmp[2].mul(lagrange[edgeV[c][e][0]]);
+
+ (*(*edge)[c])[e] = new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
+
+ (*(*edge)[c])[e]->at(0).mul(lagrange[edgeV[c][e][1]]);
+ (*(*edge)[c])[e]->at(1).mul(lagrange[edgeV[c][e][1]]);
+ (*(*edge)[c])[e]->at(2).mul(lagrange[edgeV[c][e][1]]);
+
+ (*(*edge)[c])[e]->at(0).sub(tmp[0]);
+ (*(*edge)[c])[e]->at(1).sub(tmp[1]);
+ (*(*edge)[c])[e]->at(2).sub(tmp[2]);
+ }
}
+
// 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++;
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
+
+ for(int l = 1; l < orderPlus; l++){
+ for(int e = 0; e < 6; e++){
+ (*(*edge)[c])[i + 6] =
+ new vector<Polynomial>
+ (intLegendre[l].compose
+ (lagrange[edgeV[c][e][0]] - lagrange[edgeV[c][e][1]],
+ lagrange[edgeV[c][e][0]] + lagrange[edgeV[c][e][1]]).gradient());
+
+ i++;
+ }
}
}
-
+
+
+ // Face Based //
+ for(int c = 0; c < 6; c++){
+ unsigned int i = 0;
+
+ for(int l1 = 1; l1 < order; l1++){
+ for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+ for(int f = 0; f < 4; f++){
+ // Preliminary Type 1
+ Polynomial sum =
+ lagrange[faceV[c][f][0]] +
+ lagrange[faceV[c][f][1]] +
+ lagrange[faceV[c][f][2]];
+
+ Polynomial u = intLegendre[l1].compose
+ (lagrange[faceV[c][f][0]] - lagrange[faceV[c][f][1]],
+ lagrange[faceV[c][f][0]] + lagrange[faceV[c][f][1]]);
+
+ Polynomial v = lagrange[faceV[c][f][2]] * sclLegendre[l2].compose
+ (lagrange[faceV[c][f][2]] * 2 - sum, sum);
+
+ // Preliminary Type 2
+ vector<Polynomial> gradU = u.gradient();
+ vector<Polynomial> gradV = v.gradient();
+
+ vector<Polynomial> vGradU(gradU);
+ vGradU[0].mul(v);
+ vGradU[1].mul(v);
+ vGradU[2].mul(v);
+
+ vector<Polynomial> uGradV(gradU);
+ uGradV[0].mul(u);
+ uGradV[1].mul(u);
+ uGradV[2].mul(u);
+
+ vector<Polynomial> subGradUV(vGradU);
+ subGradUV[0].sub(uGradV[0]);
+ subGradUV[1].sub(uGradV[1]);
+ subGradUV[2].sub(uGradV[2]);
+
+ // Preliminary Type 3
+ vector<Polynomial> gradL1 = lagrange[faceV[c][f][0]].gradient();
+ vector<Polynomial> gradL2 = lagrange[faceV[c][f][1]].gradient();
+
+ vector<Polynomial> l2GradL1(gradL1);
+ l2GradL1[0].mul(l2);
+ l2GradL1[1].mul(l2);
+ l2GradL1[2].mul(l2);
+
+ vector<Polynomial> l1GradL2(gradL2);
+ l1GradL2[0].mul(l1);
+ l1GradL2[1].mul(l1);
+ l1GradL2[2].mul(l1);
+
+ vector<Polynomial> subGradL1L2V(l2GradL1);
+ subGradL1L2V[0].sub(l1GradL2[0]);
+ subGradL1L2V[1].sub(l1GradL2[1]);
+ subGradL1L2V[2].sub(l1GradL2[2]);
+
+ subGradL1L2V[0].mul(v);
+ subGradL1L2V[1].mul(v);
+ subGradL1L2V[2].mul(v);
+
+
+ // Type 1
+ (*(*face)[c])[i] =
+ new vector<Polynomial>((u * v).gradient());
+
+ i++;
+
+ // Type 2
+ (*(*face)[c])[i] =
+ new vector<Polynomial>(subGradUV);
+
+ i++;
+
+ // Type 3
+ (*(*face)[c])[i] =
+ new vector<Polynomial>(subGradL1L2V);
+
+ i++;
+ }
+ }
+ }
+ }
+ /*
// Cell Based (Preliminary) //
Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
@@ -182,30 +279,46 @@ TetEdgeBasis::TetEdgeBasis(const int order){
i++;
}
-
+ */
// Free Temporary Sapce //
delete[] legendre;
+ delete[] sclLegendre;
delete[] intLegendre;
- delete[] u;
- delete[] v;
+}
+TetEdgeBasis::~TetEdgeBasis(void){
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
- // 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());
-}
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
+ }
+
+ delete edge;
-TetEdgeBasis::~TetEdgeBasis(void){
- for(int i = 0; i < size; i++){
- delete (*basis)[i];
- if(i >= nVertex && i < nVertex + nEdge)
- delete (*revBasis)[i];
+ // Face Based //
+ for(int c = 0; c < 6; c++){
+ for(int i = 0; i < nFace; i++)
+ delete (*(*face)[c])[i];
+
+ delete (*face)[c];
}
- delete basis;
- delete revBasis;
+ delete face;
+
+
+ // Cell Based //
+ for(int i = 0; i < nCell; i++)
+ delete (*cell)[i];
+
+ delete cell;
}
diff --git a/FunctionSpace/TetEdgeBasis.h b/FunctionSpace/TetEdgeBasis.h
index d8e0062f31..65d78d2afe 100644
--- a/FunctionSpace/TetEdgeBasis.h
+++ b/FunctionSpace/TetEdgeBasis.h
@@ -20,7 +20,7 @@ class TetEdgeBasis: public BasisVector{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Tetrahedra of the given order
- TetEdgeBasis(const int order);
+ TetEdgeBasis(int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index efd411af37..602900f5a6 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -3,7 +3,7 @@
using namespace std;
-TetNodeBasis::TetNodeBasis(unsigned int order){
+TetNodeBasis::TetNodeBasis(int order){
// Set Basis Type //
this->order = order;
@@ -15,7 +15,7 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
nFace = 2 * (order - 1) * (order - 2);
nCell = (order - 1) * (order - 2) * (order - 3) / 6;
- size = nVertex + nEdge + nFace + nCell;
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
@@ -50,23 +50,23 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
};
// Counter //
- int i;
+ int i;
// Basis //
- node = new vector<const Polynomial*>(nVertex);
- edge = new vector<vector<const Polynomial*>*>(2);
- face = new vector<vector<const Polynomial*>*>(6);
- cell = new vector<const Polynomial*>(nCell);
+ node = new vector<Polynomial*>(nVertex);
+ edge = new vector<vector<Polynomial*>*>(2);
+ face = new vector<vector<Polynomial*>*>(6);
+ cell = new vector<Polynomial*>(nCell);
- (*edge)[0] = new vector<const Polynomial*>(nEdge);
- (*edge)[1] = new vector<const Polynomial*>(nEdge);
+ (*edge)[0] = new vector<Polynomial*>(nEdge);
+ (*edge)[1] = new vector<Polynomial*>(nEdge);
- (*face)[0] = new vector<const Polynomial*>(nFace);
- (*face)[1] = new vector<const Polynomial*>(nFace);
- (*face)[2] = new vector<const Polynomial*>(nFace);
- (*face)[3] = new vector<const Polynomial*>(nFace);
- (*face)[4] = new vector<const Polynomial*>(nFace);
- (*face)[5] = new vector<const Polynomial*>(nFace);
+ (*face)[0] = new vector<Polynomial*>(nFace);
+ (*face)[1] = new vector<Polynomial*>(nFace);
+ (*face)[2] = new vector<Polynomial*>(nFace);
+ (*face)[3] = new vector<Polynomial*>(nFace);
+ (*face)[4] = new vector<Polynomial*>(nFace);
+ (*face)[5] = new vector<Polynomial*>(nFace);
// Vertex Based (Lagrange) //
@@ -89,7 +89,7 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
// Edge Based //
i = 0;
- for(unsigned int l = 1; l < order; l++){
+ for(int l = 1; l < order; l++){
for(int e = 0; e < 6; e++){
(*(*edge)[0])[i] =
new Polynomial(intLegendre[l].compose
diff --git a/FunctionSpace/TetNodeBasis.h b/FunctionSpace/TetNodeBasis.h
index 9ade551271..792e5cc999 100644
--- a/FunctionSpace/TetNodeBasis.h
+++ b/FunctionSpace/TetNodeBasis.h
@@ -20,7 +20,7 @@ class TetNodeBasis: public BasisScalar{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Tetrahedra of the given order
- TetNodeBasis(unsigned int order);
+ TetNodeBasis(int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 884b04531a..9cac7d26d3 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -1,7 +1,9 @@
#include "TriEdgeBasis.h"
#include "Legendre.h"
-TriEdgeBasis::TriEdgeBasis(const int order){
+using namespace std;
+
+TriEdgeBasis::TriEdgeBasis(int order){
// Set Basis Type //
this->order = order;
@@ -24,10 +26,9 @@ TriEdgeBasis::TriEdgeBasis(const int order){
Polynomial* u = new Polynomial[orderPlus];
Polynomial* v = new Polynomial[orderPlus];
- Polynomial lagrange [3];
- Polynomial lagrangeSum [3];
- Polynomial lagrangeSub [3];
- Polynomial rLagrangeSub [3];
+ Polynomial lagrange [3];
+ Polynomial lagrangeSum [3];
+ Polynomial lagrangeSub [2][3];
// Classical and Intrated-Scaled Legendre Polynomial //
Legendre::legendre(legendre, order);
@@ -45,167 +46,182 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = lagrange[j] + lagrange[i];
+ 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];
+ lagrangeSub[0][i] = lagrange[i] - lagrange[j];
+ lagrangeSub[1][i] = lagrange[j] - lagrange[i];
}
- // Basis & Revert (temporary --- *no* const) //
- std::vector<std::vector<Polynomial>*> basis(size);
- std::vector<std::vector<Polynomial>*> revBasis(size);
-
+ // Basis //
+ node = new vector<vector<Polynomial>*>(nVertex);
+ edge = new vector<vector<vector<Polynomial>*>*>(2);
+ face = new vector<vector<vector<Polynomial>*>*>(0);
+ cell = new vector<vector<Polynomial>*>(nCell);
+
+ (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
+ (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+
+ // Counter
+ unsigned int i = 0;
+
// Edge Based (Nedelec) //
- int i = 0;
-
for(int j = 1; i < 3; j = (j + 1) % 3){
- // Direct
- std::vector<Polynomial> tmp = lagrange[j].gradient();
+ // Temp
+ vector<Polynomial> tmp = lagrange[j].gradient();
+
+ // Edge[0]
tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]);
tmp[2].mul(lagrange[i]);
- basis[i] = new std::vector<Polynomial>(lagrange[i].gradient());
+ (*(*edge)[0])[i] = new 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]);
+ (*(*edge)[0])[i]->at(0).mul(lagrange[j]);
+ (*(*edge)[0])[i]->at(1).mul(lagrange[j]);
+ (*(*edge)[0])[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]);
+ (*(*edge)[0])[i]->at(0).sub(tmp[0]);
+ (*(*edge)[0])[i]->at(1).sub(tmp[1]);
+ (*(*edge)[0])[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]);
+ //Edge[1]
+ tmp = lagrange[i].gradient();
- revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
+ tmp[0].mul(lagrange[j]);
+ tmp[1].mul(lagrange[j]);
+ tmp[2].mul(lagrange[j]);
- 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]);
+ (*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
+
+ (*(*edge)[1])[i]->at(0).mul(lagrange[i]);
+ (*(*edge)[1])[i]->at(1).mul(lagrange[i]);
+ (*(*edge)[1])[i]->at(2).mul(lagrange[i]);
+
+ (*(*edge)[1])[i]->at(0).sub(tmp[0]);
+ (*(*edge)[1])[i]->at(1).sub(tmp[1]);
+ (*(*edge)[1])[i]->at(2).sub(tmp[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());
+ (*(*edge)[0])[i] =
+ new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[0][e],
+ lagrangeSum[e])).gradient());
- revBasis[i] =
- new std::vector<Polynomial>((intLegendre[l].compose(rLagrangeSub[e],
- lagrangeSum[e])).gradient());
+ (*(*edge)[1])[i] =
+ new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[1][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]);
+ u[l] = intLegendre[l].compose(lagrangeSub[0][0] * -1, lagrangeSum[0]);
v[l] = legendre[l].compose(p);
v[l].mul(lagrange[2]);
}
+ i = 0;
+
// 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();
+ 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];
+ (*cell)[i] = new vector<Polynomial>(u[l1].gradient());
+ (*cell)[i]->at(0).mul(v[l2]);
+ (*cell)[i]->at(1).mul(v[l2]);
+ (*cell)[i]->at(2).mul(v[l2]);
+
+ (*cell)[i]->at(0).add(tmp[0]);
+ (*cell)[i]->at(1).add(tmp[1]);
+ (*cell)[i]->at(2).add(tmp[2]);
+
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();
+ 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());
+ (*cell)[i] = new 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]);
+ (*cell)[i]->at(0).mul(v[l2]);
+ (*cell)[i]->at(1).mul(v[l2]);
+ (*cell)[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]);
+ (*cell)[i]->at(0).sub(tmp[0]);
+ (*cell)[i]->at(1).sub(tmp[1]);
+ (*cell)[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]);
+ (*cell)[i] = new vector<Polynomial>(*(*cell)[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];
+ (*cell)[i]->at(0).mul(v[l]);
+ (*cell)[i]->at(1).mul(v[l]);
+ (*cell)[i]->at(2).mul(v[l]);
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());
}
TriEdgeBasis::~TriEdgeBasis(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;
+
- if(i >= nVertex && i < nVertex + nEdge)
- delete (*revBasis)[i];
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[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;
}
diff --git a/FunctionSpace/TriEdgeBasis.h b/FunctionSpace/TriEdgeBasis.h
index 48b5f0be4c..e8876d5f6e 100644
--- a/FunctionSpace/TriEdgeBasis.h
+++ b/FunctionSpace/TriEdgeBasis.h
@@ -20,7 +20,7 @@ class TriEdgeBasis: public BasisVector{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Triangles of the given order
- TriEdgeBasis(const int order);
+ TriEdgeBasis(int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index dc3cae9610..c5746c35b0 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -1,5 +1,7 @@
#include "TriNedelecBasis.h"
+using namespace std;
+
TriNedelecBasis::TriNedelecBasis(void){
// Set Basis Type //
order = 1;
@@ -27,62 +29,81 @@ TriNedelecBasis::TriNedelecBasis(void){
Polynomial(1, 0, 1, 0);
- // Basis (temporary --- *no* const) //
- std::vector<std::vector<Polynomial>*> basis(size);
- std::vector<std::vector<Polynomial>*> revBasis(size);
-
+ // Basis //
+ node = new vector<vector<Polynomial>*>(nVertex);
+ edge = new vector<vector<vector<Polynomial>*>*>(2);
+ face = new vector<vector<vector<Polynomial>*>*>(0);
+ cell = new vector<vector<Polynomial>*>(nCell);
+
+ (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
+ (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+
+ // Nedelec //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
- // Direct
- std::vector<Polynomial> tmp = lagrange[j].gradient();
+ // Temp
+ vector<Polynomial> tmp = lagrange[j].gradient();
+
+ // Edge[0]
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]);
+
+ (*(*edge)[0])[i] = new vector<Polynomial>(lagrange[i].gradient());
+
+ (*(*edge)[0])[i]->at(0).mul(lagrange[j]);
+ (*(*edge)[0])[i]->at(1).mul(lagrange[j]);
+ (*(*edge)[0])[i]->at(2).mul(lagrange[j]);
+
+ (*(*edge)[0])[i]->at(0).sub(tmp[0]);
+ (*(*edge)[0])[i]->at(1).sub(tmp[1]);
+ (*(*edge)[0])[i]->at(2).sub(tmp[2]);
+
+ //Edge[1]
+ tmp = lagrange[i].gradient();
+
+ tmp[0].mul(lagrange[j]);
+ tmp[1].mul(lagrange[j]);
+ tmp[2].mul(lagrange[j]);
+
+ (*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
+
+ (*(*edge)[1])[i]->at(0).mul(lagrange[i]);
+ (*(*edge)[1])[i]->at(1).mul(lagrange[i]);
+ (*(*edge)[1])[i]->at(2).mul(lagrange[i]);
+
+ (*(*edge)[1])[i]->at(0).sub(tmp[0]);
+ (*(*edge)[1])[i]->at(1).sub(tmp[1]);
+ (*(*edge)[1])[i]->at(2).sub(tmp[2]);
}
-
-
- // 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());
}
TriNedelecBasis::~TriNedelecBasis(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];
+
+ delete (*edge)[c];
+ }
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
- delete basis;
- for(int i = 0; i < size; i++)
- delete (*revBasis)[i];
+ // Cell Based //
+ for(int i = 0; i < nCell; i++)
+ delete (*cell)[i];
- delete revBasis;
+ delete cell;
}
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 13260e0f7e..4de94cfb7f 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -3,7 +3,7 @@
using namespace std;
-TriNodeBasis::TriNodeBasis(const int order){
+TriNodeBasis::TriNodeBasis(int order){
// Set Basis Type //
this->order = order;
@@ -32,13 +32,13 @@ TriNodeBasis::TriNodeBasis(const int order){
// 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);
+ node = new vector<Polynomial*>(nVertex);
+ edge = new vector<vector<Polynomial*>*>(2);
+ face = new vector<vector<Polynomial*>*>(0);
+ cell = new vector<Polynomial*>(nCell);
- (*edge)[0] = new vector<const Polynomial*>(nEdge);
- (*edge)[1] = new vector<const Polynomial*>(nEdge);
+ (*edge)[0] = new vector<Polynomial*>(nEdge);
+ (*edge)[1] = new vector<Polynomial*>(nEdge);
// Vertex Based (Lagrange) //
@@ -56,7 +56,7 @@ TriNodeBasis::TriNodeBasis(const int order){
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = *(*node)[i] + *(*node)[j];
+ lagrangeSum[i] = *(*node)[i] + *(*node)[j];
// Lagrange Sub //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
diff --git a/FunctionSpace/TriNodeBasis.h b/FunctionSpace/TriNodeBasis.h
index 8746838d8c..018b4174c2 100644
--- a/FunctionSpace/TriNodeBasis.h
+++ b/FunctionSpace/TriNodeBasis.h
@@ -20,7 +20,7 @@ class TriNodeBasis: public BasisScalar{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Triangles of the given order
- TriNodeBasis(const int order);
+ TriNodeBasis(int order);
//! Deletes this Basis
//!
--
GitLab