diff --git a/FunctionSpace/Basis.h b/FunctionSpace/Basis.h
index 73070a1927b772fd40c7b1281017186590a47ca8..6a3badc137b5a33f8a4a3966d898b3e2f2304c73 100644
--- a/FunctionSpace/Basis.h
+++ b/FunctionSpace/Basis.h
@@ -2,6 +2,7 @@
#define _BASIS_H_
#include <string>
+#include "ReferenceSpace.h"
/**
@interface Basis
@@ -18,27 +19,29 @@
class Basis{
protected:
- bool scalar;
+ const ReferenceSpace* refSpace;
+ unsigned int nRefSpace;
+ bool scalar;
- int order;
- int type;
- int dim;
+ unsigned int order;
+ unsigned int type;
+ unsigned int dim;
- int nVertex;
- int nEdge;
- int nFace;
- int nCell;
+ unsigned int nVertex;
+ unsigned int nEdge;
+ unsigned int nFace;
+ unsigned int nCell;
- int nEdgeClosure;
- int nFaceClosure;
-
- int size;
+ unsigned int nFunction;
public:
//! Deletes this Basis
//!
virtual ~Basis(void);
+ //! @return Returns the ReferenceSpace of this Basis
+ const ReferenceSpace& getReferenceSpace(void) const;
+
//! @return Returns:
//! @li @c true, if this is a
//! @em scalar Basis (BasisScalar)
@@ -53,46 +56,39 @@ class Basis{
bool isScalar(void) const;
//! @return Returns the @em polynomial @em order of the Basis
- int getOrder(void) const;
+ unsigned int getOrder(void) const;
//! @return Returns the @em type of the Basis:
//! @li 0 for 0-form
//! @li 1 for 1-form
//! @li 2 for 2-form
//! @li 3 for 3-form
- int getType(void) const;
+ unsigned int getType(void) const;
//! @return Returns the @em dimension
//! (1D, 2D or 3D) of the Basis
- int getDim(void) const;
+ unsigned int getDim(void) const;
//! @return Returns the number of @em Vertex
//! @em Based functions of this Basis
- int getNVertexBased(void) const;
+ unsigned int getNVertexBased(void) const;
//! @return Returns the number of @em Edge
//! @em Based functions of this Basis
- int getNEdgeBased(void) const;
+ unsigned int getNEdgeBased(void) const;
//! @return Returns the number of @em Face
//! @em Based functions of this Basis
- int getNFaceBased(void) const;
+ unsigned int getNFaceBased(void) const;
//! @return Returns the number of @em Cell
//! @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;
+ unsigned int getNCellBased(void) const;
//! @return Returns the number of Polynomial%s
- //! (or Vector%s of Polynomial%s) in the Basis
- int getSize(void) const;
+ //! (or Vector%s of Polynomial%s) Functions
+ //! in this Basis
+ unsigned int getNFunction(void) const;
//! @return Returns the Basis String
virtual std::string toString(void) const = 0;
@@ -109,48 +105,45 @@ class Basis{
// Inline Functions //
//////////////////////
+inline const ReferenceSpace&
+Basis::getReferenceSpace(void) const{
+ return *refSpace;
+}
+
inline bool Basis::isScalar(void) const{
return scalar;
}
-inline int Basis::getOrder(void) const{
+inline unsigned int Basis::getOrder(void) const{
return order;
}
-inline int Basis::getType(void) const{
+inline unsigned int Basis::getType(void) const{
return type;
}
-inline int Basis::getDim(void) const{
+inline unsigned int Basis::getDim(void) const{
return dim;
}
-inline int Basis::getNVertexBased(void) const{
+inline unsigned int Basis::getNVertexBased(void) const{
return nVertex;
}
-inline int Basis::getNEdgeBased(void) const{
+inline unsigned int Basis::getNEdgeBased(void) const{
return nEdge;
}
-inline int Basis::getNFaceBased(void) const{
+inline unsigned int Basis::getNFaceBased(void) const{
return nFace;
}
-inline int Basis::getNCellBased(void) const{
+inline unsigned 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;
+inline unsigned int Basis::getNFunction(void) const{
+ return nFunction;
}
#endif
diff --git a/FunctionSpace/BasisGenerator.cpp b/FunctionSpace/BasisGenerator.cpp
index cfb3efbeb449e4e7a80bacead5d07670ea617ad0..245864cb5960fe76e3482177ed4de18ca63ba7af 100644
--- a/FunctionSpace/BasisGenerator.cpp
+++ b/FunctionSpace/BasisGenerator.cpp
@@ -27,9 +27,9 @@ BasisGenerator::BasisGenerator(void){
BasisGenerator::~BasisGenerator(void){
}
-Basis* BasisGenerator::generate(int elementType,
- int basisType,
- int order,
+Basis* BasisGenerator::generate(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order,
std::string family){
if(!family.compare(std::string("zaglmayr")))
@@ -42,10 +42,9 @@ Basis* BasisGenerator::generate(int elementType,
throw Exception("Unknwown Basis Family: %s", family.c_str());
}
-Basis* BasisGenerator::generateZaglmayr(int elementType,
- int basisType,
- int order){
-
+Basis* BasisGenerator::generateZaglmayr(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order){
switch(elementType){
case TYPE_LIN: return linZaglmayrGen(basisType, order);
case TYPE_TRI: return triZaglmayrGen(basisType, order);
@@ -58,9 +57,9 @@ Basis* BasisGenerator::generateZaglmayr(int elementType,
}
}
-Basis* BasisGenerator::generateLagrange(int elementType,
- int basisType,
- int order){
+Basis* BasisGenerator::generateLagrange(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order){
if(basisType != 0)
throw
Exception("Cannot Have a %d-Form Lagrange Basis (0-Form only)",
@@ -78,8 +77,8 @@ Basis* BasisGenerator::generateLagrange(int elementType,
}
}
-Basis* BasisGenerator::linZaglmayrGen(int basisType,
- int order){
+Basis* BasisGenerator::linZaglmayrGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
case 0: return new LineNodeBasis(order);
case 1:
@@ -93,8 +92,8 @@ Basis* BasisGenerator::linZaglmayrGen(int basisType,
}
}
-Basis* BasisGenerator::triZaglmayrGen(int basisType,
- int order){
+Basis* BasisGenerator::triZaglmayrGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
case 0: return new TriNodeBasis(order);
case 1:
@@ -108,11 +107,11 @@ Basis* BasisGenerator::triZaglmayrGen(int basisType,
}
}
-Basis* BasisGenerator::quaZaglmayrGen(int basisType,
- int order){
+Basis* BasisGenerator::quaZaglmayrGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
- case 0: return new QuadNodeBasis(order);
- case 1: return new QuadEdgeBasis(order);
+ //case 0: return new QuadNodeBasis(order);
+ //case 1: return new QuadEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Quads");
case 3: throw Exception("3-form not implemented on Quads");
@@ -120,11 +119,11 @@ Basis* BasisGenerator::quaZaglmayrGen(int basisType,
}
}
-Basis* BasisGenerator::tetZaglmayrGen(int basisType,
- int order){
+Basis* BasisGenerator::tetZaglmayrGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
- case 0: return new TetNodeBasis(order);
- case 1: return new TetEdgeBasis(order);
+ //case 0: return new TetNodeBasis(order);
+ //case 1: return new TetEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Tetrahedrons");
case 3: throw Exception("3-form not implemented on Tetrahedrons");
@@ -132,11 +131,11 @@ Basis* BasisGenerator::tetZaglmayrGen(int basisType,
}
}
-Basis* BasisGenerator::hexZaglmayrGen(int basisType,
- int order){
+Basis* BasisGenerator::hexZaglmayrGen(unsigned int basisType,
+ unsigned int order){
switch(basisType){
- case 0: return new HexNodeBasis(order);
- case 1: return new HexEdgeBasis(order);
+ //case 0: return new HexNodeBasis(order);
+ //case 1: return new HexEdgeBasis(order);
case 2: throw Exception("2-form not implemented on Hexs");
case 3: throw Exception("3-form not implemented on Hexs");
diff --git a/FunctionSpace/BasisGenerator.h b/FunctionSpace/BasisGenerator.h
index 3858bf284f30bdfc7bdcacced0b912b4820f2fed..2a6dc8c9194540c39881f4e3d07cf1f03372efd0 100644
--- a/FunctionSpace/BasisGenerator.h
+++ b/FunctionSpace/BasisGenerator.h
@@ -21,29 +21,29 @@ class BasisGenerator{
BasisGenerator(void);
~BasisGenerator(void);
- static Basis* generate(int elementType,
- int basisType,
- int order,
+ static Basis* generate(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order,
std::string family);
- static Basis* generate(int elementType,
- int basisType,
- int order);
+ static Basis* generate(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order);
- static Basis* linZaglmayrGen(int basisType, int order);
- static Basis* triZaglmayrGen(int basisType, int order);
- static Basis* quaZaglmayrGen(int basisType, int order);
- static Basis* tetZaglmayrGen(int basisType, int order);
- static Basis* hexZaglmayrGen(int basisType, int order);
+ static Basis* linZaglmayrGen(unsigned int basisType, unsigned int order);
+ static Basis* triZaglmayrGen(unsigned int basisType, unsigned int order);
+ static Basis* quaZaglmayrGen(unsigned int basisType, unsigned int order);
+ static Basis* tetZaglmayrGen(unsigned int basisType, unsigned int order);
+ static Basis* hexZaglmayrGen(unsigned int basisType, unsigned int order);
private:
- static Basis* generateZaglmayr(int elementType,
- int basisType,
- int order);
+ static Basis* generateZaglmayr(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order);
- static Basis* generateLagrange(int elementType,
- int basisType,
- int order);
+ static Basis* generateLagrange(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order);
};
@@ -60,7 +60,7 @@ class BasisGenerator{
Deletes this BasisGenerator
**
- @fn BasisGenerator::generate(int, int, int, std::string)
+ @fn BasisGenerator::generate(unsigned int, unsigned int, unsigned int, std::string)
@param elementType The type of the element,
on which the requested Basis will be created
@param basisType The Basis type
@@ -93,7 +93,7 @@ class BasisGenerator{
@li @c lagrange for Lagrange's Basis Functions
**
- @fn BasisGenerator::generate(int, int, int)
+ @fn BasisGenerator::generate(unsigned int, unsigned int, unsigned int)
@param elementType The type of the element,
on which the requested Basis will be created
@param basisType The Basis type
@@ -206,9 +206,9 @@ class BasisGenerator{
// Inline Functions //
//////////////////////
-inline Basis* BasisGenerator::generate(int elementType,
- int basisType,
- int order){
+inline Basis* BasisGenerator::generate(unsigned int elementType,
+ unsigned int basisType,
+ unsigned int order){
return BasisGenerator::generate(elementType,
basisType,
diff --git a/FunctionSpace/BasisScalar.cpp b/FunctionSpace/BasisScalar.cpp
index ff7c5607d56c4638898570cd18592843d8ccebb2..424075953ce6d77d4d65524321ff83f7de2f45e0 100644
--- a/FunctionSpace/BasisScalar.cpp
+++ b/FunctionSpace/BasisScalar.cpp
@@ -12,32 +12,27 @@ BasisScalar::~BasisScalar(void){
string BasisScalar::toString(void) const{
stringstream stream;
+ unsigned int i = 0;
stream << "Vertex Based:" << endl;
- for(int i = 0; i < nVertex; i++)
- stream << " # " << i + 1 << "\t"
- << "f(" << i + 1 << ") = "
- << (*node)[i]->toString() << endl;
-
- stream << "Edge Based:" << endl;
- for(int i = 0; i < nEdge; i++)
- for(int j = 0; j < nEdgeClosure; j++)
- stream << " # " << i + 1 + nVertex << "\t"
- << "f_" << j << "(" << i + 1 << ") = "
- << (*(*edge)[j])[i]->toString() << endl;
-
- stream << "Face Based:" << endl;
- for(int i = 0; i < nFace; i++)
- for(int j = 0; j < nFaceClosure; j++)
- stream << " # " << i + 1 + nVertex + nEdge << "\t"
- << "f_" << j << "(" << i + 1 << ") = "
- << (*(*face)[j])[i]->toString() << endl;
-
- stream << "Cell Based:" << endl;
- for(int i = 0; i < nCell; i++)
- stream << " # " << i + 1 + nVertex + nEdge + nFace << "\t"
- << "f(" << i + 1 << ") = "
- << (*cell)[i]->toString() << endl;
+ for(; i < nVertex; i++)
+ stream << "f(" << i + 1 << ") = "
+ << (*(*basis)[0])[i]->toString() << endl;
+
+ stream << "Edge Based:" << endl;
+ for(; i < nVertex + nEdge; i++)
+ stream << "f(" << i + 1 << ") = "
+ << (*(*basis)[0])[i]->toString() << endl;
+
+ stream << "Face Based:" << endl;
+ for(; i < nVertex + nEdge + nFace; i++)
+ stream << "f(" << i + 1 << ") = "
+ << (*(*basis)[0])[i]->toString() << endl;
+
+ stream << "Cell Based:" << endl;
+ for(; i < nVertex + nEdge + nFace + nCell; i++)
+ stream << "f(" << i + 1 << ") = "
+ << (*(*basis)[0])[i]->toString() << endl;
return stream.str();
}
diff --git a/FunctionSpace/BasisScalar.h b/FunctionSpace/BasisScalar.h
index 97497043144232c2ea4f4199b76c47efad9d5977..9b4b39f3460af37a010945a85409310d874891e9 100644
--- a/FunctionSpace/BasisScalar.h
+++ b/FunctionSpace/BasisScalar.h
@@ -20,41 +20,32 @@
class BasisScalar: public Basis{
protected:
- std::vector <Polynomial*>* node;
- std::vector<std::vector<Polynomial*>*>* edge;
- std::vector<std::vector<Polynomial*>*>* face;
- std::vector <Polynomial*>* cell;
+ std::vector<std::vector<const Polynomial*>*>* basis;
public:
//! Deletes this BasisScalar
//!
virtual ~BasisScalar(void);
+ //! @param refSpace A natural number
//! @param i A natural number
- //! @return Returns the @c i%th @em Vertex Based
- //! Basis Function
+ //! @return Returns the @c i%th @em
+ //! Basis Function of the @c refSpace%th ReferenceSpace
const Polynomial&
- getNodeFunction(unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the @c i%th @em Edge Based
- //! Basis Function, with the @c closure%th Closure
- const Polynomial&
- getEdgeFunction(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the @c i%th @em Face Based
- //! Basis Function, with the @c closure%th Closure
- const Polynomial&
- getFaceFunction(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @return Returns the @c i%th @em Cell Based
- //! Basis Function
- const Polynomial&
- getCellFunction(unsigned int i) const;
+ getFunction(unsigned int refSpace, unsigned int i) const;
+
+ //! @param refSpace A natural number
+ //! @return Returns the @em all
+ //! Basis Function of the @c refSpace%th ReferenceSpace
+ const std::vector<const Polynomial*>&
+ getFunction(unsigned int refSpace) const;
+
+ //! @param element An Element
+ //! @return Returns the @em all
+ //! Basis Function in the @c given element
+ //! @em ReferenceSpace
+ const std::vector<const Polynomial*>&
+ getFunction(const MElement& element) const;
virtual std::string toString(void) const;
@@ -70,28 +61,22 @@ class BasisScalar: public Basis{
// Inline Function //
//////////////////////
-inline
-const Polynomial&
-BasisScalar::getNodeFunction(unsigned int i) const{
- return *(*node)[i];
-}
-
inline
const Polynomial&
-BasisScalar::getEdgeFunction(unsigned int closure, unsigned int i) const{
- return *(*(*edge)[closure])[i];
+BasisScalar::getFunction(unsigned int refSpace, unsigned int i) const{
+ return *(*(*basis)[refSpace])[i];
}
-inline
-const Polynomial&
-BasisScalar::getFaceFunction(unsigned int closure, unsigned int i) const{
- return *(*(*face)[closure])[i];
+inline
+const std::vector<const Polynomial*>&
+BasisScalar::getFunction(unsigned int refSpace) const{
+ return *(*basis)[refSpace];
}
-inline
-const Polynomial&
-BasisScalar::getCellFunction(unsigned int i) const{
- return *(*cell)[i];
+inline
+const std::vector<const Polynomial*>&
+BasisScalar::getFunction(const MElement& element) const{
+ return *(*basis)[refSpace->getReferenceSpace(element)];
}
#endif
diff --git a/FunctionSpace/BasisVector.cpp b/FunctionSpace/BasisVector.cpp
index 2b2649ebe6f659f4fed7d52a1d1d54facb547f07..bcbbd19c011ee921d4b11345f09e07b0302d6f25 100644
--- a/FunctionSpace/BasisVector.cpp
+++ b/FunctionSpace/BasisVector.cpp
@@ -12,40 +12,35 @@ BasisVector::~BasisVector(void){
string BasisVector::toString(void) const{
stringstream stream;
+ unsigned int i = 0;
stream << "Vertex Based:" << endl;
- for(int i = 0; i < nVertex; i++)
- stream << " # " << i + 1 << "\t"
- << "f(" << i + 1 << ") = " << endl
- << "\t [ " <<(*(*node)[i])[0].toString() << " ]" << endl
- << "\t [ " <<(*(*node)[i])[1].toString() << " ]" << endl
- << "\t [ " <<(*(*node)[i])[2].toString() << " ]" << endl;
-
- stream << "Edge Based:" << endl;
- for(int i = 0; i < nEdge; i++)
- for(int j = 0; j < nEdgeClosure; j++)
- stream << " # " << i + 1 + nVertex << "\t"
- << " f_" << j << "(" << i + 1 << ") = " << endl
- << "\t [ " << (*(*(*edge)[j])[i])[0].toString() << " ]" << endl
- << "\t [ " << (*(*(*edge)[j])[i])[1].toString() << " ]" << endl
- << "\t [ " << (*(*(*edge)[j])[i])[2].toString() << " ]" << endl;
-
- stream << "Face Based:" << endl;
- for(int i = 0; i < nFace; i++)
- for(int j = 0; j < nFaceClosure; j++)
- stream << " # " << i + 1 + nVertex + nEdge << "\t"
- << " f_" << j << "(" << i + 1 << ") = " << endl
- << "\t [ " << (*(*(*face)[j])[i])[0].toString() << " ]" << endl
- << "\t [ " << (*(*(*face)[j])[i])[1].toString() << " ]" << endl
- << "\t [ " << (*(*(*face)[j])[i])[2].toString() << " ]" << endl;
+ for(; i < nVertex; i++)
+ stream << "f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[0].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[1].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[2].toString() << " ]" << endl;
+
+ stream << "Edge Based:" << endl;
+ for(; i < nVertex + nEdge; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[0].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[1].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[2].toString() << " ]" << endl;
+
+ stream << "Face Based:" << endl;
+ for(; i < nVertex + nEdge + nFace; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[0].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[1].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[2].toString() << " ]" << endl;
- stream << "Cell Based:" << endl;
- for(int i = 0; i < nCell; i++)
- stream << " # " << i + 1 + nVertex + nEdge << "\t"
- << " f(" << i + 1 << ") = " << endl
- << "\t [ " << (*(*cell)[i])[0].toString() << " ]" << endl
- << "\t [ " << (*(*cell)[i])[1].toString() << " ]" << endl
- << "\t [ " << (*(*cell)[i])[2].toString() << " ]" << endl;
+ stream << "Cell Based:" << endl;
+ for(; i < nVertex + nEdge + nFace + nCell; i++)
+ stream << " f(" << i + 1 << ") = " << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[0].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[1].toString() << " ]" << endl
+ << "\t[ " << (*(*(*basis)[0])[i])[2].toString() << " ]" << endl;
return stream.str();
}
diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
index d36b0e30e256b2e30b4e18606fe641e5ade4e389..d202662723cb7410aaeef793522a805af55ed91a 100644
--- a/FunctionSpace/BasisVector.h
+++ b/FunctionSpace/BasisVector.h
@@ -20,41 +20,32 @@
class BasisVector: public Basis{
protected:
- 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<std::vector<const std::vector<Polynomial>*>*>* basis;
public:
//! Deletes this BasisVector
//!
virtual ~BasisVector(void);
+ //! @param refSpace A natural number
//! @param i A natural number
- //! @return Returns the @c i%th @em Vertex Based
- //! Basis Function
+ //! @return Returns the @c i%th @em
+ //! Basis Function of the @c refSpace%th ReferenceSpace
const std::vector<Polynomial>&
- getNodeFunction(unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the @c i%th @em Edge Based
- //! Basis Function, with the @c closure%th Closure
- const std::vector<Polynomial>&
- getEdgeFunction(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the @c i%th @em Face Based
- //! Basis Function, with the @c closure%th Closure
- const std::vector<Polynomial>&
- getFaceFunction(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @return Returns the @c i%th @em Cell Based
- //! Basis Function
- const std::vector<Polynomial>&
- getCellFunction(unsigned int i) const;
+ getFunction(unsigned int refSpace, unsigned int i) const;
+
+ //! @param refSpace A natural number
+ //! @return Returns the @em all
+ //! Basis Function of the @c refSpace%th ReferenceSpace
+ const std::vector<const std::vector<Polynomial>*>&
+ getFunction(unsigned int refSpace) const;
+
+ //! @param element An Element
+ //! @return Returns the @em all
+ //! Basis Function in the @c given element
+ //! @em ReferenceSpace
+ const std::vector<const std::vector<Polynomial>*>&
+ getFunction(const MElement& element) const;
virtual std::string toString(void) const;
@@ -70,28 +61,22 @@ class BasisVector: public Basis{
// Inline Functions //
//////////////////////
-inline
-const std::vector<Polynomial>&
-BasisVector::getNodeFunction(unsigned int i) const{
- return *(*node)[i];
-}
-
inline
const std::vector<Polynomial>&
-BasisVector::getEdgeFunction(unsigned int closure, unsigned int i) const{
- return *(*(*edge)[closure])[i];
+BasisVector::getFunction(unsigned int refSpace, unsigned int i) const{
+ return *(*(*basis)[refSpace])[i];
}
-inline
-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::getFunction(unsigned int refSpace) const{
+ return *(*basis)[refSpace];
}
-inline
-const std::vector<Polynomial>&
-BasisVector::getCellFunction(unsigned int i) const{
- return *(*cell)[i];
+inline
+const std::vector<const std::vector<Polynomial>*>&
+BasisVector::getFunction(const MElement& element) const{
+ return *(*basis)[refSpace->getReferenceSpace(element)];
}
#endif
diff --git a/FunctionSpace/CMakeLists.txt b/FunctionSpace/CMakeLists.txt
index c930d90f37a5f151e79f503496fed76f5782128e..1fd77ad570c87390eea96326b9b92e19e8064908 100644
--- a/FunctionSpace/CMakeLists.txt
+++ b/FunctionSpace/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SRC
Legendre.cpp
ReferenceSpace.cpp
+ LineReferenceSpace.cpp
TriReferenceSpace.cpp
TetReferenceSpace.cpp
@@ -22,26 +23,24 @@ set(SRC
DivBasis.cpp
EvaluatedBasis.cpp
- EvaluatedBasisScalar.cpp
- EvaluatedBasisVector.cpp
LineNodeBasis.cpp
LineEdgeBasis.cpp
LineNedelecBasis.cpp
- QuadNodeBasis.cpp
- QuadEdgeBasis.cpp
+# QuadNodeBasis.cpp
+# QuadEdgeBasis.cpp
TriNodeBasis.cpp
TriEdgeBasis.cpp
TriNedelecBasis.cpp
TriLagrangeBasis.cpp
- HexNodeBasis.cpp
- HexEdgeBasis.cpp
+# HexNodeBasis.cpp
+# HexEdgeBasis.cpp
- TetNodeBasis.cpp
- TetEdgeBasis.cpp
+# TetNodeBasis.cpp
+# TetEdgeBasis.cpp
FunctionSpace.cpp
FunctionSpaceScalar.cpp
diff --git a/FunctionSpace/CurlBasis.cpp b/FunctionSpace/CurlBasis.cpp
index bb4d2546f2c6a96266df4b6da73cbe3636815c17..dae7fd5c3c194829558436f7830f45f97e40f7bb 100644
--- a/FunctionSpace/CurlBasis.cpp
+++ b/FunctionSpace/CurlBasis.cpp
@@ -3,94 +3,43 @@
using namespace std;
CurlBasis::CurlBasis(const BasisVector& other){
+ // Reference Space //
+ refSpace = &other.getReferenceSpace();
+ nRefSpace = other.getReferenceSpace().getNReferenceSpace();
+
// Set Basis Type //
order = other.getOrder() - 1;
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();
+ nVertex = other.getNVertexBased();
+ nEdge = other.getNEdgeBased();
+ nFace = other.getNFaceBased();
+ nCell = other.getNCellBased();
+ nFunction = other.getNFunction();
- // 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);
+ // Basis //
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
- for(int i = 0; i < nEdgeClosure; i++)
- (*edge)[i] = new vector<vector<Polynomial>*>(nEdge);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
- 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] =
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int i = 0; i < nFunction; i++)
+ (*(*basis)[s])[i] =
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)));
+ (Polynomial::curl(other.getFunction(s, i)));
}
CurlBasis::~CurlBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
-
- // 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 (*basis)[i];
}
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
-
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/DivBasis.cpp b/FunctionSpace/DivBasis.cpp
index 8a992c8106c1827686025baf9d8ac94e3a342468..34679481960bc42d71c40a99c5ae5117ac59a18d 100644
--- a/FunctionSpace/DivBasis.cpp
+++ b/FunctionSpace/DivBasis.cpp
@@ -3,94 +3,43 @@
using namespace std;
DivBasis::DivBasis(const BasisVector& other){
+ // Reference Space //
+ refSpace = &other.getReferenceSpace();
+ nRefSpace = other.getReferenceSpace().getNReferenceSpace();
+
// Set Basis Type //
order = other.getOrder() - 1;
type = other.getType();
dim = other.getDim();
- nVertex = other.getNVertexBased();
- nEdge = other.getNEdgeBased();
- nFace = other.getNFaceBased();
- nCell = other.getNCellBased();
+ nVertex = other.getNVertexBased();
+ nEdge = other.getNEdgeBased();
+ nFace = other.getNFaceBased();
+ nCell = other.getNCellBased();
+ nFunction = other.getNFunction();
- 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);
+ // Basis //
+ basis = new vector<vector<const Polynomial*>*>(nRefSpace);
- // 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)));
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const Polynomial*>(nFunction);
- // Face Based //
- for(int i = 0; i < nFaceClosure; i++)
- for(int j = 0; j < nFace; j++)
- (*(*face)[i])[j] =
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int i = 0; i < nFunction; i++)
+ (*(*basis)[s])[i] =
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)));
+ (Polynomial::divergence(other.getFunction(s, 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];
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- delete (*edge)[c];
+ delete (*basis)[i];
}
- 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;
+ delete basis;
}
diff --git a/FunctionSpace/EvaluatedBasis.cpp b/FunctionSpace/EvaluatedBasis.cpp
index 3b06376ebb6af29eb34bb6df8c47c7db0a43c1d2..c6e7927210ca1dae4b68120596adea5bbdcdc6ea 100644
--- a/FunctionSpace/EvaluatedBasis.cpp
+++ b/FunctionSpace/EvaluatedBasis.cpp
@@ -1,7 +1,75 @@
+#include "Polynomial.h"
#include "EvaluatedBasis.h"
-EvaluatedBasis::EvaluatedBasis(void){
+using namespace std;
+
+EvaluatedBasis::
+EvaluatedBasis(const BasisScalar& basis,
+ const fullMatrix<double>& point){
+ // Data //
+ refSpace = &basis.getReferenceSpace();
+ nRefSpace = refSpace->getNReferenceSpace();
+ scalar = true;
+ nFunction = basis.getNFunction();
+ nPoint = point.size1();
+
+ // Alloc //
+ eBasis = new vector<fullMatrix<double>*>(nRefSpace);
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ (*eBasis)[i] = new fullMatrix<double>(nFunction, nPoint);
+
+ // Evaluate //
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ for(unsigned int j = 0; j < nFunction; j++)
+ for(unsigned int k = 0; k < nPoint; k++)
+ (*(*eBasis)[i])(j, k) =
+ basis.getFunction(i, j).at(point(k, 0),
+ point(k, 1),
+ point(k, 2));
+}
+
+EvaluatedBasis::
+EvaluatedBasis(const BasisVector& basis,
+ const fullMatrix<double>& point){
+ // Data //
+ refSpace = &basis.getReferenceSpace();
+ nRefSpace = refSpace->getNReferenceSpace();
+ scalar = false;
+ nFunction = basis.getNFunction();
+ nPoint = point.size1();
+
+ // Alloc //
+ eBasis = new vector<fullMatrix<double>*>(nRefSpace);
+
+ // WARNING Each Evaluation Got *3* Component !
+ const unsigned int nPointThree = nPoint * 3;
+
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ (*eBasis)[i] = new fullMatrix<double>(nFunction, nPointThree);
+
+ // Evaluate //
+ fullVector<double> tmp(3);
+
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++){
+ for(unsigned int k = 0; k < nPoint; k++){
+ tmp = Polynomial::at(basis.getFunction(i, j),
+ point(k, 0),
+ point(k, 1),
+ point(k, 2));
+
+ (*(*eBasis)[i])(j, 3 * k) = tmp(0);
+ (*(*eBasis)[i])(j, 3 * k + 1) = tmp(1);
+ (*(*eBasis)[i])(j, 3 * k + 2) = tmp(2);
+ }
+ }
+ }
}
EvaluatedBasis::~EvaluatedBasis(void){
+ for(unsigned int i = 0; i < nRefSpace; i++)
+ delete (*eBasis)[i];
+
+ delete eBasis;
}
diff --git a/FunctionSpace/EvaluatedBasis.h b/FunctionSpace/EvaluatedBasis.h
index 6967d2106560a0517d736ec0132ea39b25a7be5b..a6e4b50de2b6b545fea9f31fe200c26d91bf19a5 100644
--- a/FunctionSpace/EvaluatedBasis.h
+++ b/FunctionSpace/EvaluatedBasis.h
@@ -1,34 +1,52 @@
#ifndef _EVALUATEDBASIS_H_
#define _EVALUATEDBASIS_H_
-/**
- @interface EvaluatedBasis
- @brief Common Interface for Evaluated Basis
+#include <vector>
+#include "fullMatrix.h"
+#include "BasisScalar.h"
+#include "BasisVector.h"
+#include "MElement.h"
+#include "ReferenceSpace.h"
- This class is the @em common @em interface for all EvaluatedBasis.@n
+/**
+ @class EvaluatedBasis
+ @brief A Basis that has been Evaluated
- An EvaluatedBasis is a @em Basis that has been @em Evaluated at
- some points.@n
+ This class is A basis that has been Evaluated at
+ some given points.@n
*/
class EvaluatedBasis{
- protected:
- bool scalar;
-
- unsigned int nVertex;
- unsigned int nEdge;
- unsigned int nFace;
- unsigned int nCell;
-
- unsigned int nEdgeClosure;
- unsigned int nFaceClosure;
+ private:
+ const ReferenceSpace* refSpace;
+ unsigned int nRefSpace;
+ bool scalar;
+ unsigned int nFunction;
unsigned int nPoint;
+ std::vector<fullMatrix<double>*>* eBasis;
+
public:
+ //! @param basis the Basis to Evaluate
+ //! @param point the Evaluation Points
+ //!
+ //! Instanciates the requested EvaluatedBasisScalar
+ //!
+ EvaluatedBasis(const BasisScalar& basis,
+ const fullMatrix<double>& point);
+
+ //! @param basis the Basis to Evaluate
+ //! @param point the Evaluation Points
+ //!
+ //! Instanciates the requested EvaluatedBasisVector
+ //!
+ EvaluatedBasis(const BasisVector& basis,
+ const fullMatrix<double>& point);
+
//! Deletes this EvaluatedBasis
//!
- virtual ~EvaluatedBasis(void);
+ ~EvaluatedBasis(void);
//! @return Returns:
//! @li @c true, if the evaluated
@@ -36,76 +54,53 @@ class EvaluatedBasis{
//! @li @c false, otherwise
bool isScalar(void) const;
- //! @return Returns the number of @em Vertex
- //! @em Based functions of this EvaluatedBasis
- unsigned int getNVertexBased(void) const;
-
- //! @return Returns the number of @em Edge
- //! @em Based functions of this EvaluatedBasis
- unsigned int getNEdgeBased(void) const;
-
- //! @return Returns the number of @em Face
- //! @em Based functions of this EvaluatedBasis
- unsigned int getNFaceBased(void) const;
-
- //! @return Returns the number of @em Cell
- //! @em Based functions of this EvaluatedBasis
- unsigned int getNCellBased(void) const;
-
- //! @return Returns the number of @em edge
- //! @em closures for this EvaluatedBasis
- unsigned int getNEdgeClosure(void) const;
-
- //! @return Returns the number of @em face
- //! @em closures for this EvaluatedBasis
- unsigned int getNFaceClosure(void) const;
+ //! @return Returns the number of
+ //! evaluated @em Functions
+ unsigned int getNFunction(void) const;
//! @return Returns the number of
- //! evaluation Points
+ //! evaluation @em Points
unsigned int getNPoint(void) const;
- protected:
- //! @internal
- //! Instantiate a new EvaluatedBasis
- //!
- //! @endinternal
- EvaluatedBasis(void);
+ //! @param refSpace A natural number
+ //! @return Returns the evaluation of all Basis Functions
+ //! (at every evaluation Points)
+ //! for the @c refSpace%th @em ReferenceSpace
+ const fullMatrix<double>&
+ getEvaluation(unsigned int refSpace) const;
+
+ //! @param element A MElement
+ //! @return Returns the evaluation of all Basis Functions
+ //! (at every evaluation Points)
+ //! in the ReferenceSpace of the given element
+ const fullMatrix<double>&
+ getEvaluation(const MElement& element) const;
};
//////////////////////
-// Inline Functions //
+// Inline Function //
//////////////////////
inline bool EvaluatedBasis::isScalar(void) const{
return scalar;
}
-inline unsigned int EvaluatedBasis::getNVertexBased(void) const{
- return nVertex;
-}
-
-inline unsigned int EvaluatedBasis::getNEdgeBased(void) const{
- return nEdge;
-}
-
-inline unsigned int EvaluatedBasis::getNFaceBased(void) const{
- return nFace;
-}
-
-inline unsigned int EvaluatedBasis::getNCellBased(void) const{
- return nCell;
+inline unsigned int EvaluatedBasis::getNFunction(void) const{
+ return nFunction;
}
-inline unsigned int EvaluatedBasis::getNEdgeClosure(void) const{
- return nEdgeClosure;
+inline unsigned int EvaluatedBasis::getNPoint(void) const{
+ return nPoint;
}
-inline unsigned int EvaluatedBasis::getNFaceClosure(void) const{
- return nFaceClosure;
+inline const fullMatrix<double>&
+EvaluatedBasis::getEvaluation(unsigned int refSpace) const{
+ return *(*eBasis)[refSpace];
}
-inline unsigned int EvaluatedBasis::getNPoint(void) const{
- return nPoint;
+inline const fullMatrix<double>&
+EvaluatedBasis::getEvaluation(const MElement& element) const{
+ return *(*eBasis)[refSpace->getReferenceSpace(element)];
}
#endif
diff --git a/FunctionSpace/EvaluatedBasisScalar.cpp b/FunctionSpace/EvaluatedBasisScalar.cpp
deleted file mode 100644
index 036e819e079198c31068b8a995a65b906c5ddb79..0000000000000000000000000000000000000000
--- a/FunctionSpace/EvaluatedBasisScalar.cpp
+++ /dev/null
@@ -1,124 +0,0 @@
-#include "EvaluatedBasisScalar.h"
-
-using namespace std;
-
-EvaluatedBasisScalar::
-EvaluatedBasisScalar(const BasisScalar& basis,
- const fullMatrix<double>& point){
- // Data //
- scalar = true;
-
- nVertex = basis.getNVertexBased();
- nEdge = basis.getNEdgeBased();
- nFace = basis.getNFaceBased();
- nCell = basis.getNCellBased();
-
- nEdgeClosure = basis.getNEdgeClosure();
- nFaceClosure = basis.getNFaceClosure();
-
- nPoint = point.size1();
-
-
- // Alloc //
- // Node
- node = new vector<vector<double>*>(nVertex);
-
- for(unsigned int j = 0; j < nVertex; j++)
- (*node)[j] = new vector<double>(nPoint);
-
- // Edge
- edge = new vector<vector<vector<double>*>*>(nEdgeClosure);
-
- for(unsigned int i = 0; i < nEdgeClosure; i++){
- (*edge)[i] = new vector<vector<double>*>(nEdge);
-
- for(unsigned int j = 0; j < nEdge; j++)
- (*(*edge)[i])[j] = new vector<double>(nPoint);
- }
-
- // Face
- face = new vector<vector<vector<double>*>*>(nFaceClosure);
-
- for(unsigned int i = 0; i < nFaceClosure; i++){
- (*face)[i] = new vector<vector<double>*>(nFace);
-
- for(unsigned int j = 0; j < nFace; j++)
- (*(*face)[i])[j] = new vector<double>(nPoint);
- }
-
- // Cell
- cell = new vector<vector<double>*>(nCell);
-
- for(unsigned int j = 0; j < nCell; j++)
- (*cell)[j] = new vector<double>(nPoint);
-
-
- // Evaluate //
- // Node
- for(unsigned int j = 0; j < nVertex; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*node)[j])[k] =
- basis.getNodeFunction(j).at(point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Edge
- for(unsigned int i = 0; i < nEdgeClosure; i++)
- for(unsigned int j = 0; j < nEdge; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*(*edge)[i])[j])[k] =
- basis.getEdgeFunction(i, j).at(point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Face
- for(unsigned int i = 0; i < nFaceClosure; i++)
- for(unsigned int j = 0; j < nFace; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*(*face)[i])[j])[k] =
- basis.getFaceFunction(i, j).at(point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Cell
- for(unsigned int j = 0; j < nCell; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*cell)[j])[k] =
- basis.getCellFunction(j).at(point(k, 0),
- point(k, 1),
- point(k, 2));
-}
-
-EvaluatedBasisScalar::~EvaluatedBasisScalar(void){
- // Node //
- for(unsigned int j = 0; j < nVertex; j++)
- delete (*node)[j];
-
- delete node;
-
- // Edge //
- for(unsigned int i = 0; i < nEdgeClosure; i++){
- for(unsigned int j = 0; j < nEdge; j++)
- delete (*(*edge)[i])[j];
-
- delete (*edge)[i];
- }
-
- delete edge;
-
- // Face //
- for(unsigned int i = 0; i < nFaceClosure; i++){
- for(unsigned int j = 0; j < nFace; j++)
- delete (*(*face)[i])[j];
-
- delete (*face)[i];
- }
-
- delete face;
-
- // Cell //
- for(unsigned int j = 0; j < nCell; j++)
- delete (*cell)[j];
-
- delete cell;
-}
diff --git a/FunctionSpace/EvaluatedBasisScalar.h b/FunctionSpace/EvaluatedBasisScalar.h
deleted file mode 100644
index c3a5aff88a2975221673475609c8b442a6faf208..0000000000000000000000000000000000000000
--- a/FunctionSpace/EvaluatedBasisScalar.h
+++ /dev/null
@@ -1,91 +0,0 @@
-#ifndef _EVALUATEDBASISSCALAR_H_
-#define _EVALUATEDBASISSCALAR_H_
-
-#include <vector>
-#include "fullMatrix.h"
-#include "BasisScalar.h"
-#include "EvaluatedBasis.h"
-
-/**
- @class EvaluatedBasisScalar
- @brief A @em Scalar EvaluatedBasis
-
- This class is a @em Scalar EvaluatedBasis.@n
-*/
-
-class EvaluatedBasisScalar: public EvaluatedBasis{
- protected:
- std::vector <std::vector<double>*>* node;
- std::vector<std::vector<std::vector<double>*>*>* edge;
- std::vector<std::vector<std::vector<double>*>*>* face;
- std::vector <std::vector<double>*>* cell;
-
- public:
- //! @param basis the Basis to Evaluate
- //! @param point the Evaluation Points
- //!
- //! Instanciates the requested EvaluatedBasisScalar
- //!
- EvaluatedBasisScalar(const BasisScalar& basis,
- const fullMatrix<double>& point);
-
- //! Deletes this EvaluatedBasisScalar
- //!
- virtual ~EvaluatedBasisScalar(void);
-
- //! @param i A natural number
- //! @return Returns the evaluation of the @c i%th @em Vertex Based
- //! Basis Function (at every evaluation Points)
- const std::vector<double>&
- getNodeEvaluation(unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the evaluation of the @c i%th @em Edge Based
- //! Basis Function (at every evaluation Points), with the @c closure%th Closure
- const std::vector<double>&
- getEdgeEvaluation(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the evaluation of the @c i%th @em Face Based
- //! Basis Function (at every evaluation Points), with the @c closure%th Closure
- const std::vector<double>&
- getFaceEvaluation(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @return Returns the evaluation of the @c i%th @em Cell Based
- //! Basis Function (at every evaluation Points)
- const std::vector<double>&
- getCellEvaluation(unsigned int i) const;
-};
-
-//////////////////////
-// Inline Function //
-//////////////////////
-
-inline
-const std::vector<double>&
-EvaluatedBasisScalar::getNodeEvaluation(unsigned int i) const{
- return *(*node)[i];
-}
-
-inline
-const std::vector<double>&
-EvaluatedBasisScalar::getEdgeEvaluation(unsigned int closure, unsigned int i) const{
- return *(*(*edge)[closure])[i];
-}
-
-inline
-const std::vector<double>&
-EvaluatedBasisScalar::getFaceEvaluation(unsigned int closure, unsigned int i) const{
- return *(*(*face)[closure])[i];
-}
-
-inline
-const std::vector<double>&
-EvaluatedBasisScalar::getCellEvaluation(unsigned int i) const{
- return *(*cell)[i];
-}
-
-#endif
diff --git a/FunctionSpace/EvaluatedBasisVector.cpp b/FunctionSpace/EvaluatedBasisVector.cpp
deleted file mode 100644
index ac437d38aab985dadd62a2c70c571ea280cb8450..0000000000000000000000000000000000000000
--- a/FunctionSpace/EvaluatedBasisVector.cpp
+++ /dev/null
@@ -1,129 +0,0 @@
-#include "Polynomial.h"
-#include "EvaluatedBasisVector.h"
-
-using namespace std;
-
-EvaluatedBasisVector::
-EvaluatedBasisVector(const BasisVector& basis,
- const fullMatrix<double>& point){
- // Data //
- scalar = false;
-
- nVertex = basis.getNVertexBased();
- nEdge = basis.getNEdgeBased();
- nFace = basis.getNFaceBased();
- nCell = basis.getNCellBased();
-
- nEdgeClosure = basis.getNEdgeClosure();
- nFaceClosure = basis.getNFaceClosure();
-
- nPoint = point.size1();
-
-
- // Alloc //
- // Node
- node = new vector<vector<fullVector<double> >*>(nVertex);
-
- for(unsigned int j = 0; j < nVertex; j++)
- (*node)[j] = new vector<fullVector<double> >(nPoint);
-
- // Edge
- edge = new vector<vector<vector<fullVector<double> >*>*>(nEdgeClosure);
-
- for(unsigned int i = 0; i < nEdgeClosure; i++){
- (*edge)[i] = new vector<vector<fullVector<double> >*>(nEdge);
-
- for(unsigned int j = 0; j < nEdge; j++)
- (*(*edge)[i])[j] = new vector<fullVector<double> >(nPoint);
- }
-
- // Face
- face = new vector<vector<vector<fullVector<double> >*>*>(nFaceClosure);
-
- for(unsigned int i = 0; i < nFaceClosure; i++){
- (*face)[i] = new vector<vector<fullVector<double> >*>(nFace);
-
- for(unsigned int j = 0; j < nFace; j++)
- (*(*face)[i])[j] = new vector<fullVector<double> >(nPoint);
- }
-
- // Cell
- cell = new vector<vector<fullVector<double> >*>(nCell);
-
- for(unsigned int j = 0; j < nCell; j++)
- (*cell)[j] = new vector<fullVector<double> >(nPoint);
-
-
- // Evaluate //
- // Node
- for(unsigned int j = 0; j < nVertex; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*node)[j])[k] =
- Polynomial::at(basis.getNodeFunction(j),
- point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Edge
- for(unsigned int i = 0; i < nEdgeClosure; i++)
- for(unsigned int j = 0; j < nEdge; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*(*edge)[i])[j])[k] =
- Polynomial::at(basis.getEdgeFunction(i, j),
- point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Face
- for(unsigned int i = 0; i < nFaceClosure; i++)
- for(unsigned int j = 0; j < nFace; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*(*face)[i])[j])[k] =
- Polynomial::at(basis.getFaceFunction(i, j),
- point(k, 0),
- point(k, 1),
- point(k, 2));
-
- // Cell
- for(unsigned int j = 0; j < nCell; j++)
- for(unsigned int k = 0; k < nPoint; k++)
- (*(*cell)[j])[k] =
- Polynomial::at(basis.getCellFunction(j),
- point(k, 0),
- point(k, 1),
- point(k, 2));
-}
-
-EvaluatedBasisVector::~EvaluatedBasisVector(void){
- // Node //
- for(unsigned int j = 0; j < nVertex; j++)
- delete (*node)[j];
-
- delete node;
-
- // Edge //
- for(unsigned int i = 0; i < nEdgeClosure; i++){
- for(unsigned int j = 0; j < nEdge; j++)
- delete (*(*edge)[i])[j];
-
- delete (*edge)[i];
- }
-
- delete edge;
-
- // Face //
- for(unsigned int i = 0; i < nFaceClosure; i++){
- for(unsigned int j = 0; j < nFace; j++)
- delete (*(*face)[i])[j];
-
- delete (*face)[i];
- }
-
- delete face;
-
- // Cell //
- for(unsigned int j = 0; j < nCell; j++)
- delete (*cell)[j];
-
- delete cell;
-}
diff --git a/FunctionSpace/EvaluatedBasisVector.h b/FunctionSpace/EvaluatedBasisVector.h
deleted file mode 100644
index c5e8c0afaaa53cb2804b95c1321ce8d3ff37e3d8..0000000000000000000000000000000000000000
--- a/FunctionSpace/EvaluatedBasisVector.h
+++ /dev/null
@@ -1,91 +0,0 @@
-#ifndef _EVALUATEDBASISVECTOR_H_
-#define _EVALUATEDBASISVECTOR_H_
-
-#include <vector>
-#include "fullMatrix.h"
-#include "BasisVector.h"
-#include "EvaluatedBasis.h"
-
-/**
- @class EvaluatedBasisVector
- @brief A @em Vectorial EvaluatedBasis
-
- This class is a @em Vectorial EvaluatedBasis.@n
-*/
-
-class EvaluatedBasisVector: public EvaluatedBasis{
- protected:
- std::vector <std::vector<fullVector<double> >*>* node;
- std::vector<std::vector<std::vector<fullVector<double> >*>*>* edge;
- std::vector<std::vector<std::vector<fullVector<double> >*>*>* face;
- std::vector <std::vector<fullVector<double> >*>* cell;
-
- public:
- //! @param basis the Basis to Evaluate
- //! @param point the Evaluation Points
- //!
- //! Instanciates the requested EvaluatedBasisVector
- //!
- EvaluatedBasisVector(const BasisVector& basis,
- const fullMatrix<double>& point);
-
- //! Deletes this EvaluatedBasisVector
- //!
- virtual ~EvaluatedBasisVector(void);
-
- //! @param i A natural number
- //! @return Returns the evaluation of the @c i%th @em Vertex Based
- //! Basis Function (at every evaluation Points)
- const std::vector<fullVector<double> >&
- getNodeEvaluation(unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the evaluation of the @c i%th @em Edge Based
- //! Basis Function (at every evaluation Points), with the @c closure%th Closure
- const std::vector<fullVector<double> >&
- getEdgeEvaluation(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @param closure A natural number
- //! @return Returns the evaluation of the @c i%th @em Face Based
- //! Basis Function (at every evaluation Points), with the @c closure%th Closure
- const std::vector<fullVector<double> >&
- getFaceEvaluation(unsigned int closure, unsigned int i) const;
-
- //! @param i A natural number
- //! @return Returns the evaluation of the @c i%th @em Cell Based
- //! Basis Function (at every evaluation Points)
- const std::vector<fullVector<double> >&
- getCellEvaluation(unsigned int i) const;
-};
-
-//////////////////////
-// Inline Function //
-//////////////////////
-
-inline
-const std::vector<fullVector<double> >&
-EvaluatedBasisVector::getNodeEvaluation(unsigned int i) const{
- return *(*node)[i];
-}
-
-inline
-const std::vector<fullVector<double> >&
-EvaluatedBasisVector::getEdgeEvaluation(unsigned int closure, unsigned int i) const{
- return *(*(*edge)[closure])[i];
-}
-
-inline
-const std::vector<fullVector<double> >&
-EvaluatedBasisVector::getFaceEvaluation(unsigned int closure, unsigned int i) const{
- return *(*(*face)[closure])[i];
-}
-
-inline
-const std::vector<fullVector<double> >&
-EvaluatedBasisVector::getCellEvaluation(unsigned int i) const{
- return *(*cell)[i];
-}
-
-#endif
diff --git a/FunctionSpace/FunctionSpace.cpp b/FunctionSpace/FunctionSpace.cpp
index a0c6d9848f06d85061a54ff09b04ede32c4774fb..79ae8688c628e44c705449498462604353dd5938 100644
--- a/FunctionSpace/FunctionSpace.cpp
+++ b/FunctionSpace/FunctionSpace.cpp
@@ -221,378 +221,6 @@ const GroupOfDof& FunctionSpace::getGoDFromElement(const MElement& element) cons
return *(it->second);
}
-vector<int> FunctionSpace::getEdgeClosure(const MElement& element){
- // 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){
- // 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;
-}
-
-const vector<const Polynomial*>
-FunctionSpace::locBasis(const MElement& element,
- const BasisScalar& basis){
- // Get Basis //
- 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 Closure //
- vector<int> edgeClosure = getEdgeClosure(element);
- vector<int> faceClosure = getFaceClosure(element);
-
- // Get Functions //
- vector<const Polynomial*> fun(basis.getSize());
- unsigned int i = 0;
-
- // 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();
-
- if(nEdge){
- 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();
-
- if(nFace){
- 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++;
- }
-
-
- // Return //
- return fun;
-}
-
-const vector<const vector<Polynomial>*>
-FunctionSpace::locBasis(const MElement& element,
- const BasisVector& basis){
- // Get Basis //
- 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 Closure //
- vector<int> edgeClosure = getEdgeClosure(element);
- vector<int> faceClosure = getFaceClosure(element);
-
- // Get Functions //
- vector<const vector<Polynomial>*> fun(basis.getSize());
- unsigned int i = 0;
-
- // 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();
-
- if(nEdge){
- 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();
-
- if(nFace){
- 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++;
- }
-
-
- // Return //
- return fun;
-}
-
-const vector<const vector<double>*>
-FunctionSpace::locEvalBasis(const MElement& element,
- const EvaluatedBasisScalar& evalBasis){
- // Get Basis //
- const unsigned int nFNode = evalBasis.getNVertexBased();
- const unsigned int nFEdge = evalBasis.getNEdgeBased();
- const unsigned int nFFace = evalBasis.getNFaceBased();
- const unsigned int nFCell = evalBasis.getNCellBased();
-
- // Get Closure //
- vector<int> edgeClosure = getEdgeClosure(element);
- vector<int> faceClosure = getFaceClosure(element);
-
- // Get Functions //
- vector<const vector<double>*> fun(nFNode + nFEdge + nFFace + nFCell);
- unsigned int i = 0;
-
- // Vertex Based
- for(unsigned int j = 0; j < nFNode; j++){
- fun[i] = &evalBasis.getNodeEvaluation(j);
- i++;
- }
-
- // Edge Based
- // Number of basis function *per* edge
- // --> should always be an integer !
- const unsigned int nEdge = edgeClosure.size();
-
- if(nEdge){
- 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] =
- &evalBasis.getEdgeEvaluation(edgeClosure[k], fEdge);
-
- fEdge++;
- i++;
- }
- }
- }
-
- // Face Based
- // Number of basis function *per* face
- // --> should always be an integer !
- const unsigned int nFace = faceClosure.size();
-
- if(nFace){
- 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] =
- &evalBasis.getFaceEvaluation(faceClosure[k], fFace);
-
- fFace++;
- i++;
- }
- }
- }
-
- // Cell Based
- for(unsigned int j = 0; j < nFCell; j++){
- fun[i] = &evalBasis.getCellEvaluation(j);
- i++;
- }
-
-
- // Return //
- return fun;
-}
-
-const vector<const vector<fullVector<double> >*>
-FunctionSpace::locEvalBasis(const MElement& element,
- const EvaluatedBasisVector& evalBasis){
- // Get Basis //
- const unsigned int nFNode = evalBasis.getNVertexBased();
- const unsigned int nFEdge = evalBasis.getNEdgeBased();
- const unsigned int nFFace = evalBasis.getNFaceBased();
- const unsigned int nFCell = evalBasis.getNCellBased();
-
- // Get Closure //
- vector<int> edgeClosure = getEdgeClosure(element);
- vector<int> faceClosure = getFaceClosure(element);
-
- // Get Functions //
- vector<const vector<fullVector<double> >*>
- fun(nFNode + nFEdge + nFFace + nFCell);
-
- unsigned int i = 0;
-
- // Vertex Based
- for(unsigned int j = 0; j < nFNode; j++){
- fun[i] = &evalBasis.getNodeEvaluation(j);
- i++;
- }
-
- // Edge Based
- // Number of basis function *per* edge
- // --> should always be an integer !
- const unsigned int nEdge = edgeClosure.size();
-
- if(nEdge){
- 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] =
- &evalBasis.getEdgeEvaluation(edgeClosure[k], fEdge);
-
- fEdge++;
- i++;
- }
- }
- }
-
- // Face Based
- // Number of basis function *per* face
- // --> should always be an integer !
- const unsigned int nFace = faceClosure.size();
-
- if(nFace){
- 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] =
- &evalBasis.getFaceEvaluation(faceClosure[k], fFace);
-
- fFace++;
- i++;
- }
- }
- }
-
- // Cell Based
- for(unsigned int j = 0; j < nFCell; j++){
- fun[i] = &evalBasis.getCellEvaluation(j);
- i++;
- }
-
-
- // Return //
- return fun;
-}
-
string FunctionSpace::toString(void) const{
return basis->toString();
}
diff --git a/FunctionSpace/FunctionSpace.h b/FunctionSpace/FunctionSpace.h
index 82d7dd64ad41af6d6a28d15e01da329896883f7a..f73d663cffa0bca6f40857fa9978b4609d395673 100644
--- a/FunctionSpace/FunctionSpace.h
+++ b/FunctionSpace/FunctionSpace.h
@@ -8,9 +8,7 @@
#include "Basis.h"
#include "BasisScalar.h"
#include "BasisVector.h"
-
-#include "EvaluatedBasisScalar.h"
-#include "EvaluatedBasisVector.h"
+#include "EvaluatedBasis.h"
#include "Comparators.h"
#include "Dof.h"
@@ -107,29 +105,20 @@ class FunctionSpace{
void buildDof(void);
void insertDof(Dof& d, GroupOfDof* god);
- // Closure
- static std::vector<int> getEdgeClosure(const MElement& element);
- static std::vector<int> getFaceClosure(const MElement& element);
-
// Local Basis
static
- const std::vector<const Polynomial*>
+ const std::vector<const Polynomial*>&
locBasis(const MElement& element,
const BasisScalar& basis);
static
- const std::vector<const std::vector<Polynomial>*>
+ const std::vector<const std::vector<Polynomial>*>&
locBasis(const MElement& element,
const BasisVector& basis);
static
- const std::vector<const std::vector<double>*>
+ const fullMatrix<double>&
locEvalBasis(const MElement& element,
- const EvaluatedBasisScalar& evalBasis);
-
- static
- const std::vector<const std::vector<fullVector<double> >*>
- locEvalBasis(const MElement& element,
- const EvaluatedBasisVector& evalBasis);
+ const EvaluatedBasis& evalBasis);
};
@@ -299,4 +288,25 @@ inline const std::vector<GroupOfDof*>& FunctionSpace::getAllGroups(void) const{
return *group;
}
+inline const std::vector<const Polynomial*>&
+FunctionSpace::locBasis(const MElement& element,
+ const BasisScalar& basis){
+
+ return basis.getFunction(element);
+}
+
+inline const std::vector<const std::vector<Polynomial>*>&
+FunctionSpace::locBasis(const MElement& element,
+ const BasisVector& basis){
+
+ return basis.getFunction(element);
+}
+
+inline const fullMatrix<double>&
+FunctionSpace::locEvalBasis(const MElement& element,
+ const EvaluatedBasis& evalBasis){
+
+ return evalBasis.getEvaluation(element);
+}
+
#endif
diff --git a/FunctionSpace/FunctionSpaceScalar.cpp b/FunctionSpace/FunctionSpaceScalar.cpp
index 84759eff62ad8f89c80a7aef5807a331761b2789..cf828ab28e4b623614f3a422ecc062116c12a894 100644
--- a/FunctionSpace/FunctionSpaceScalar.cpp
+++ b/FunctionSpace/FunctionSpaceScalar.cpp
@@ -30,7 +30,7 @@ preEvaluateLocalFunctions(fullMatrix<double>& points){
delete evalLoc;
// New Struct //
- evalLoc = new EvaluatedBasisScalar(*basisScalar, points);
+ evalLoc = new EvaluatedBasis(*basisScalar, points);
// PreEvaluated //
locPreEvaluated = true;
@@ -51,7 +51,7 @@ preEvaluateGradLocalFunctions(fullMatrix<double>& points){
delete evalGrad;
// New Struct //
- evalGrad = new EvaluatedBasisVector(*gradBasis, points);
+ evalGrad = new EvaluatedBasis(*gradBasis, points);
// PreEvaluated //
gradPreEvaluated = true;
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 795ddbb152120af7bc6a13ca278df06c6e7be06b..979252dbb62fd62e76af1c8a9f78658ab5d75445 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -5,8 +5,7 @@
#include "Exception.h"
#include "BasisScalar.h"
#include "GradBasis.h"
-#include "EvaluatedBasisScalar.h"
-#include "EvaluatedBasisVector.h"
+#include "EvaluatedBasis.h"
#include "FunctionSpace.h"
/**
@@ -36,8 +35,8 @@ class FunctionSpaceScalar : public FunctionSpace{
bool locPreEvaluated;
bool gradPreEvaluated;
- EvaluatedBasisScalar* evalLoc;
- EvaluatedBasisVector* evalGrad;
+ EvaluatedBasis* evalLoc;
+ EvaluatedBasis* evalGrad;
public:
virtual ~FunctionSpaceScalar(void);
@@ -52,19 +51,19 @@ class FunctionSpaceScalar : public FunctionSpace{
const std::vector<double>& coef,
const fullVector<double>& uvw) const = 0;
- const std::vector<const Polynomial*>
+ const std::vector<const Polynomial*>&
getLocalFunctions(const MElement& element) const;
- const std::vector<const std::vector<Polynomial>*>
+ const std::vector<const std::vector<Polynomial>*>&
getGradLocalFunctions(const MElement& element) const;
void preEvaluateLocalFunctions(fullMatrix<double>& points);
void preEvaluateGradLocalFunctions(fullMatrix<double>& points);
- const std::vector<const std::vector<double>*>
+ const fullMatrix<double>&
getEvaluatedLocalFunctions(const MElement& element) const;
- const std::vector<const std::vector<fullVector<double> >*>
+ const fullMatrix<double>&
getEvaluatedGradLocalFunctions(const MElement& element) const;
protected:
@@ -178,12 +177,12 @@ class FunctionSpaceScalar : public FunctionSpace{
// Inline Functions //
//////////////////////
-inline const std::vector<const Polynomial*>
+inline const std::vector<const Polynomial*>&
FunctionSpaceScalar::getLocalFunctions(const MElement& element) const{
return locBasis(element, *basisScalar);
}
-inline const std::vector<const std::vector<Polynomial>*>
+inline const std::vector<const std::vector<Polynomial>*>&
FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
// Got Grad Basis ? //
@@ -197,7 +196,7 @@ FunctionSpaceScalar::getGradLocalFunctions(const MElement& element) const{
return locBasis(element, *gradBasis);
}
-inline const std::vector<const std::vector<double>*>
+inline const fullMatrix<double>&
FunctionSpaceScalar::getEvaluatedLocalFunctions(const MElement& element) const{
if(!locPreEvaluated)
throw Exception("Local Basis Functions not PreEvaluated");
@@ -205,7 +204,7 @@ FunctionSpaceScalar::getEvaluatedLocalFunctions(const MElement& element) const{
return locEvalBasis(element, *evalLoc);
}
-inline const std::vector<const std::vector<fullVector<double> >*>
+inline const fullMatrix<double>&
FunctionSpaceScalar::getEvaluatedGradLocalFunctions(const MElement& element) const{
if(!gradPreEvaluated)
throw Exception("Gradients of Local Basis Functions not PreEvaluated");
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index c58f421b4f534cf34749f6259724c94ab66e801e..fb4dbbd7e91bb5c5a0294b8b2ee4ad0a822bea9f 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -41,7 +41,7 @@ preEvaluateLocalFunctions(fullMatrix<double>& points){
delete evalLoc;
// New Struct //
- evalLoc = new EvaluatedBasisVector(*basisVector, points);
+ evalLoc = new EvaluatedBasis(*basisVector, points);
// PreEvaluated //
locPreEvaluated = true;
@@ -62,7 +62,7 @@ preEvaluateCurlLocalFunctions(fullMatrix<double>& points){
delete evalCurl;
// New Struct //
- evalCurl = new EvaluatedBasisVector(*curlBasis, points);
+ evalCurl = new EvaluatedBasis(*curlBasis, points);
// PreEvaluated //
curlPreEvaluated = true;
@@ -83,7 +83,7 @@ preEvaluateDivLocalFunctions(fullMatrix<double>& points){
delete evalDiv;
// New Struct //
- evalDiv = new EvaluatedBasisScalar(*divBasis, points);
+ evalDiv = new EvaluatedBasis(*divBasis, points);
// PreEvaluated //
divPreEvaluated = true;
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 80250611d942522b7e12203c709a23534a0f25eb..134ff39b3bf023d10c6defff33ba7acb583fa4b8 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -6,8 +6,7 @@
#include "BasisVector.h"
#include "CurlBasis.h"
#include "DivBasis.h"
-#include "EvaluatedBasisScalar.h"
-#include "EvaluatedBasisVector.h"
+#include "EvaluatedBasis.h"
#include "FunctionSpace.h"
/**
@@ -41,9 +40,9 @@ class FunctionSpaceVector : public FunctionSpace{
bool curlPreEvaluated;
bool divPreEvaluated;
- EvaluatedBasisVector* evalLoc;
- EvaluatedBasisVector* evalCurl;
- EvaluatedBasisScalar* evalDiv;
+ EvaluatedBasis* evalLoc;
+ EvaluatedBasis* evalCurl;
+ EvaluatedBasis* evalDiv;
public:
virtual ~FunctionSpaceVector(void);
@@ -58,26 +57,26 @@ class FunctionSpaceVector : public FunctionSpace{
const std::vector<double>& coef,
const fullVector<double>& uvw) const = 0;
- const std::vector<const std::vector<Polynomial>*>
+ const std::vector<const std::vector<Polynomial>*>&
getLocalFunctions(const MElement& element) const;
- const std::vector<const std::vector<Polynomial>*>
+ const std::vector<const std::vector<Polynomial>*>&
getCurlLocalFunctions(const MElement& element) const;
- const std::vector<const Polynomial*>
+ const std::vector<const Polynomial*>&
getDivLocalFunctions(const MElement& element) const;
void preEvaluateLocalFunctions(fullMatrix<double>& points);
void preEvaluateCurlLocalFunctions(fullMatrix<double>& points);
void preEvaluateDivLocalFunctions(fullMatrix<double>& points);
- const std::vector<const std::vector<fullVector<double> >*>
+ const fullMatrix<double>&
getEvaluatedLocalFunctions(const MElement& element) const;
- const std::vector<const std::vector<fullVector<double> >*>
+ const fullMatrix<double>&
getEvaluatedCurlLocalFunctions(const MElement& element) const;
- const std::vector<const std::vector<double>*>
+ const fullMatrix<double>&
getEvaluatedDivLocalFunctions(const MElement& element) const;
protected:
@@ -221,12 +220,12 @@ class FunctionSpaceVector : public FunctionSpace{
// Inline Functions //
//////////////////////
-inline const std::vector<const std::vector<Polynomial>*>
+inline const std::vector<const std::vector<Polynomial>*>&
FunctionSpaceVector::getLocalFunctions(const MElement& element) const{
return locBasis(element, *basisVector);
}
-inline const std::vector<const std::vector<Polynomial>*>
+inline const std::vector<const std::vector<Polynomial>*>&
FunctionSpaceVector::getCurlLocalFunctions(const MElement& element) const{
// Got Curl Basis ? //
@@ -240,7 +239,7 @@ FunctionSpaceVector::getCurlLocalFunctions(const MElement& element) const{
return locBasis(element, *curlBasis);
}
-inline const std::vector<const Polynomial*>
+inline const std::vector<const Polynomial*>&
FunctionSpaceVector::getDivLocalFunctions(const MElement& element) const{
// Got Div Basis ? //
@@ -254,7 +253,7 @@ FunctionSpaceVector::getDivLocalFunctions(const MElement& element) const{
return locBasis(element, *divBasis);
}
-inline const std::vector<const std::vector<fullVector<double> >*>
+inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedLocalFunctions(const MElement& element) const{
if(!locPreEvaluated)
throw Exception("Local Basis Functions not PreEvaluated");
@@ -262,7 +261,7 @@ FunctionSpaceVector::getEvaluatedLocalFunctions(const MElement& element) const{
return locEvalBasis(element, *evalLoc);
}
-inline const std::vector<const std::vector<fullVector<double> >*>
+inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedCurlLocalFunctions(const MElement& element) const{
if(!curlPreEvaluated)
throw Exception("Curls of Local Basis Functions not PreEvaluated");
@@ -270,7 +269,7 @@ FunctionSpaceVector::getEvaluatedCurlLocalFunctions(const MElement& element) con
return locEvalBasis(element, *evalCurl);
}
-inline const std::vector<const std::vector<double>*>
+inline const fullMatrix<double>&
FunctionSpaceVector::getEvaluatedDivLocalFunctions(const MElement& element) const{
if(!divPreEvaluated)
throw Exception("Divergences of Local Basis Functions not PreEvaluated");
diff --git a/FunctionSpace/GradBasis.cpp b/FunctionSpace/GradBasis.cpp
index 96b303bfa9abaec0df9d01f6b1ca0089705dc32f..68830057e63f5724ab384c0594914566fac1992d 100644
--- a/FunctionSpace/GradBasis.cpp
+++ b/FunctionSpace/GradBasis.cpp
@@ -3,94 +3,43 @@
using namespace std;
GradBasis::GradBasis(const BasisScalar& other){
+ // Reference Space //
+ refSpace = &other.getReferenceSpace();
+ nRefSpace = other.getReferenceSpace().getNReferenceSpace();
+
// Set Basis Type //
order = other.getOrder() - 1;
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();
+ nVertex = other.getNVertexBased();
+ nEdge = other.getNEdgeBased();
+ nFace = other.getNFaceBased();
+ nCell = other.getNCellBased();
+ nFunction = other.getNFunction();
- // 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);
+ // Basis //
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
- for(int i = 0; i < nEdgeClosure; i++)
- (*edge)[i] = new vector<vector<Polynomial>*>(nEdge);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
- 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] =
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ for(unsigned int i = 0; i < nFunction; i++)
+ (*(*basis)[s])[i] =
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());
+ (other.getFunction(s, i).gradient());
}
GradBasis::~GradBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
-
- // 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 (*basis)[i];
}
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
-
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/LineEdgeBasis.cpp b/FunctionSpace/LineEdgeBasis.cpp
index e88105404248ca2da641a1bc1e1546f8b9a39ba6..0dcf2c85b80f23aa77347d7be7321d8698d2ddf5 100644
--- a/FunctionSpace/LineEdgeBasis.cpp
+++ b/FunctionSpace/LineEdgeBasis.cpp
@@ -1,32 +1,42 @@
#include "LineEdgeBasis.h"
+#include "LineReferenceSpace.h"
#include "Legendre.h"
using namespace std;
-LineEdgeBasis::LineEdgeBasis(int order){
+LineEdgeBasis::LineEdgeBasis(unsigned int order){
+ // Reference Space //
+ refSpace = new LineReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
+ const vector<const vector<const vector<unsigned int>*>*>&
+ edgeV = refSpace->getAllEdge();
+
// Set Basis Type //
this->order = order;
type = 1;
dim = 1;
- nVertex = 0;
- nEdge = (order + 1);
- nFace = 0;
- nCell = 0;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = nVertex + nEdge + nFace + nCell;
+ nVertex = 0;
+ nEdge = (order + 1);
+ nFace = 0;
+ nCell = 0;
+ nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
- const int orderPlus = order + 1;
+ const unsigned int orderPlus = order + 1;
Polynomial* intLegendre = new Polynomial[orderPlus];
- const Polynomial plusOneHalf(+0.5, 0, 0, 0);
- const Polynomial minusOneHalf(-0.5, 0, 0, 0);
- const Polynomial zero(0, 0, 0, 0);
+ vector<Polynomial> first(3);
+ first[0] = Polynomial(+0.5, 0, 0, 0);
+ first[1] = Polynomial( 0 , 0, 0, 0);
+ first[2] = Polynomial( 0 , 0, 0, 0);
+
+ vector<Polynomial> second(3);
+ second[0] = Polynomial(-0.5, 0, 0, 0);
+ second[1] = Polynomial( 0 , 0, 0, 0);
+ second[2] = Polynomial( 0 , 0, 0, 0);
const Polynomial x[2] = {
Polynomial(+1, 1, 0, 0),
@@ -36,71 +46,44 @@ LineEdgeBasis::LineEdgeBasis(int order){
// Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus);
- // Permutations //
- const int permutation[2] = {0, 1};
-
// 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);
-
- (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
- (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
// Edge Based (Nedelec) //
- (*(*edge)[0])[0] = new vector<Polynomial>(3);
- (*(*edge)[0])[0]->at(0) = plusOneHalf;
- (*(*edge)[0])[0]->at(1) = zero;
- (*(*edge)[0])[0]->at(2) = zero;
-
- (*(*edge)[1])[0] = new vector<Polynomial>(3);
- (*(*edge)[1])[0]->at(0) = minusOneHalf;
- (*(*edge)[1])[0]->at(1) = zero;
- (*(*edge)[1])[0]->at(2) = zero;
-
+ (*(*basis)[0])[0] = new vector<Polynomial>(first);
+ (*(*basis)[1])[0] = new vector<Polynomial>(second);
// Edge Based (High Order) //
- for(int c = 0; c < 2; c++){
- unsigned int i = 0;
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = 1;
- for(int l = 1; l < orderPlus; l++){
- (*(*edge)[c])[i + 1] =
- new vector<Polynomial>((intLegendre[l].compose(x[permutation[c]])).gradient());
+ for(unsigned int l = 1; l < orderPlus; l++){
+ (*(*basis)[s])[i] =
+ new vector<Polynomial>((intLegendre[l].compose
+ (x[(*(*edgeV[s])[0])[0]])).gradient());
i++;
}
}
-
// Free Temporary Space //
delete[] intLegendre;
}
LineEdgeBasis::~LineEdgeBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
+ delete (*basis)[i];
}
-
- delete edge;
-
-
- // Face Based //
- delete face;
-
- // Cell Based //
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/LineEdgeBasis.h b/FunctionSpace/LineEdgeBasis.h
index c5e7e794f46e1983fc334e01860103d2af55e9d8..7986b2a7df77296cbbda8fd47a7456ead65864a8 100644
--- a/FunctionSpace/LineEdgeBasis.h
+++ b/FunctionSpace/LineEdgeBasis.h
@@ -24,7 +24,7 @@ class LineEdgeBasis: public BasisVector{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Lines of the given order
- LineEdgeBasis(int order);
+ LineEdgeBasis(unsigned int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/LineNedelecBasis.cpp b/FunctionSpace/LineNedelecBasis.cpp
index 04b9a79533136ad8cd90ab0edd19be080718dc3f..f0d2f16bcb1b3ba5195f3f1c9fa57a5b0e61bf9b 100644
--- a/FunctionSpace/LineNedelecBasis.cpp
+++ b/FunctionSpace/LineNedelecBasis.cpp
@@ -1,75 +1,59 @@
#include "LineNedelecBasis.h"
+#include "LineReferenceSpace.h"
#include "Legendre.h"
using namespace std;
LineNedelecBasis::LineNedelecBasis(void){
+ // Reference Space //
+ refSpace = new LineReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
// Set Basis Type //
order = 1;
type = 1;
dim = 1;
- nVertex = 0;
- nEdge = 1;
- nFace = 0;
- nCell = 0;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = 1;
+ nVertex = 0;
+ nEdge = 1;
+ nFace = 0;
+ nCell = 0;
+ nFunction = 1;
// Alloc Temporary Space //
- const Polynomial plusOneHalf(+0.5, 0, 0, 0);
- const Polynomial minusOneHalf(-0.5, 0, 0, 0);
- const Polynomial zero(0, 0, 0, 0);
+ vector<Polynomial> first(3);
+ first[0] = Polynomial(+0.5, 0, 0, 0);
+ first[1] = Polynomial( 0 , 0, 0, 0);
+ first[2] = Polynomial( 0 , 0, 0, 0);
+
+ vector<Polynomial> second(3);
+ second[0] = Polynomial(-0.5, 0, 0, 0);
+ second[1] = Polynomial( 0 , 0, 0, 0);
+ second[2] = Polynomial( 0 , 0, 0, 0);
// 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);
-
- (*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
- (*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
// Nedelec //
- (*(*edge)[0])[0] = new vector<Polynomial>(3);
- (*(*edge)[0])[0]->at(0) = plusOneHalf;
- (*(*edge)[0])[0]->at(1) = zero;
- (*(*edge)[0])[0]->at(2) = zero;
-
- (*(*edge)[1])[0] = new vector<Polynomial>(3);
- (*(*edge)[1])[0]->at(0) = minusOneHalf;
- (*(*edge)[1])[0]->at(1) = zero;
- (*(*edge)[1])[0]->at(2) = zero;
+ (*(*basis)[0])[0] = new vector<Polynomial>(first);
+ (*(*basis)[1])[0] = new vector<Polynomial>(second);
}
LineNedelecBasis::~LineNedelecBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
+ delete (*basis)[i];
}
-
- delete edge;
-
-
- // Face Based //
- delete face;
-
- // Cell Based //
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/LineNodeBasis.cpp b/FunctionSpace/LineNodeBasis.cpp
index 4b94c8c5c2456e3810891619c3703ed05482eb7f..ec3bd6e41cd395d028118091a117ecb6ec0347e6 100644
--- a/FunctionSpace/LineNodeBasis.cpp
+++ b/FunctionSpace/LineNodeBasis.cpp
@@ -1,24 +1,28 @@
#include "LineNodeBasis.h"
+#include "LineReferenceSpace.h"
#include "Legendre.h"
using namespace std;
-LineNodeBasis::LineNodeBasis(int order){
+LineNodeBasis::LineNodeBasis(unsigned int order){
+ // Reference Space //
+ refSpace = new LineReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
+ const vector<const vector<const vector<unsigned int>*>*>&
+ edgeV = refSpace->getAllEdge();
+
// Set Basis Type //
this->order = order;
type = 0;
dim = 1;
- nVertex = 2;
- nEdge = (order - 1);
- nFace = 0;
- nCell = 0;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = nVertex + nEdge + nFace + nCell;
+ nVertex = 2;
+ nEdge = (order - 1);
+ nFace = 0;
+ nCell = 0;
+ nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* intLegendre = new Polynomial[order];
@@ -31,69 +35,50 @@ LineNodeBasis::LineNodeBasis(int order){
// Legendre Polynomial //
Legendre::integrated(intLegendre, order);
- // Permutations //
- const int permutation[2] = {0, 1};
-
// Basis //
- node = new vector<Polynomial*>(nVertex);
- edge = new vector<vector<Polynomial*>*>(nEdgeClosure);
- face = new vector<vector<Polynomial*>*>(nFaceClosure);
- cell = new vector<Polynomial*>(nCell);
-
- (*edge)[0] = new vector<Polynomial*>(nEdge);
- (*edge)[1] = new vector<Polynomial*>(nEdge);
+ basis = new vector<vector<const Polynomial*>*>(nRefSpace);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const Polynomial*>(nFunction);
// Vertex Based (Lagrange) //
- (*node)[0] =
- new Polynomial(Polynomial(0.5, 0, 0, 0) -
- Polynomial(0.5, 1, 0, 0));
-
- (*node)[1] =
- new Polynomial(Polynomial(0.5, 0, 0, 0) +
- Polynomial(0.5, 1, 0, 0));
-
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ (*(*basis)[s])[0] =
+ new Polynomial(Polynomial(0.5, 0, 0, 0) -
+ Polynomial(0.5, 1, 0, 0));
+
+ (*(*basis)[s])[1] =
+ new Polynomial(Polynomial(0.5, 0, 0, 0) +
+ Polynomial(0.5, 1, 0, 0));
+ }
// Edge Based //
- for(int c = 0; c < 2; c++){
- unsigned int i = 0;
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nVertex;
- for(int l = 1; l < order; l++){
- (*(*edge)[c])[i] =
- new Polynomial(intLegendre[l].compose(x[permutation[c]]));
+ for(unsigned int l = 1; l < order; l++){
+ (*(*basis)[s])[i] =
+ new Polynomial(intLegendre[l].compose(x[(*(*edgeV[s])[0])[0]]));
i++;
}
}
-
// Free Temporary Sapce //
delete[] intLegendre;
}
LineNodeBasis::~LineNodeBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
+ delete (*basis)[i];
}
-
- delete edge;
-
-
- // Face Based //
- delete face;
-
- // Cell Based //
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/LineNodeBasis.h b/FunctionSpace/LineNodeBasis.h
index 3b07ecaa8b39d260e66e1234672f00f220a24c14..9bc6bbea9fcb378aea7db56e433a1116e35b4393 100644
--- a/FunctionSpace/LineNodeBasis.h
+++ b/FunctionSpace/LineNodeBasis.h
@@ -24,7 +24,7 @@ class LineNodeBasis: public BasisScalar{
//! @param order The order of the Basis
//!
//! Returns a new Node-Basis for Lines of the given order
- LineNodeBasis(int order);
+ LineNodeBasis(unsigned int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/LineReferenceSpace.cpp b/FunctionSpace/LineReferenceSpace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b2be2283b293487998ca9c4064305463ae0ccff
--- /dev/null
+++ b/FunctionSpace/LineReferenceSpace.cpp
@@ -0,0 +1,69 @@
+#include <sstream>
+#include "LineReferenceSpace.h"
+
+using namespace std;
+
+LineReferenceSpace::LineReferenceSpace(void){
+ // Vertex Definition //
+ nVertex = 2;
+
+ // Edge Definition //
+ nEdge = 1;
+ refEdge = new unsigned int*[nEdge];
+
+ for(unsigned int i = 0; i < nEdge; i++)
+ refEdge[i] = new unsigned int[2];
+
+ refEdge[0][0] = 0; refEdge[0][1] = 1;
+
+ // Face Definition //
+ nFace = 0;
+ refFace = NULL;
+
+ // Init All //
+ init();
+}
+
+LineReferenceSpace::~LineReferenceSpace(void){
+ // Delete Ref Edge //
+ for(unsigned int i = 0; i < nEdge; i++)
+ delete[] refEdge[i];
+
+ delete[] refEdge;
+}
+
+string LineReferenceSpace::toLatex(void) const{
+ stringstream stream;
+
+ stream << "\\documentclass{article}" << endl << endl
+
+ << "\\usepackage{longtable}" << endl
+ << "\\usepackage{tikz}" << endl
+ << "\\usetikzlibrary{arrows}" << endl << endl
+
+ << "\\begin{document}" << endl
+ << "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
+ << "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
+
+ << "\\begin{longtable}{c}" << endl << endl;
+
+ for(unsigned int p = 0; p < nPerm; p++){
+ stream << "\\begin{tikzpicture}" << endl
+
+ << "\\node[vertex] (n0) at(0, 0) {$0$};" << endl
+ << "\\node[vertex] (n1) at(3, 0) {$1$};" << endl << endl
+
+ << "\\path[line]"
+ << " (n" << (*(*(*edge)[p])[0])[0] << ")"
+ << " -- "
+ << " (n" << (*(*(*edge)[p])[0])[1] << ");"
+ << endl
+
+ << "\\end{tikzpicture} \\\\ \\\\" << endl << endl;
+ }
+
+ stream << "\\end{longtable}" << endl
+ << "\\end{document}" << endl;
+
+ return stream.str();
+}
diff --git a/FunctionSpace/LineReferenceSpace.h b/FunctionSpace/LineReferenceSpace.h
new file mode 100644
index 0000000000000000000000000000000000000000..16f0cc5c6f9bc830bac19f455eb34f26ed0a2712
--- /dev/null
+++ b/FunctionSpace/LineReferenceSpace.h
@@ -0,0 +1,15 @@
+#ifndef _LINEREFERENCESPACE_H_
+#define _LINEREFERENCESPACE_H_
+
+#include <string>
+#include "ReferenceSpace.h"
+
+class LineReferenceSpace: public ReferenceSpace{
+ public:
+ LineReferenceSpace(void);
+ virtual ~LineReferenceSpace(void);
+
+ std::string toLatex(void) const;
+};
+
+#endif
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index b1069425586158e22a502bce801d23c2b72d8364..7d22185907cbb9e85ba1438050a54dfc6bb5e24d 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -1,132 +1,109 @@
#include "TriEdgeBasis.h"
+#include "TriReferenceSpace.h"
#include "Legendre.h"
using namespace std;
-TriEdgeBasis::TriEdgeBasis(int order){
+TriEdgeBasis::TriEdgeBasis(unsigned int order){
+ // Reference Space //
+ refSpace = new TriReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
+ const vector<const vector<const vector<unsigned int>*>*>&
+ edgeV = refSpace->getAllEdge();
+
// Set Basis Type //
this->order = order;
type = 1;
dim = 2;
- nVertex = 0;
- nEdge = 3 * (order + 1);
- nFace = 0;
- nCell = ((order - 1) * order + order - 1);
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = nVertex + nEdge + nFace + nCell;
+ nVertex = 0;
+ nEdge = 3 * (order + 1);
+ nFace = 0;
+ nCell = ((order - 1) * order + order - 1);
+ nFunction = nVertex + nEdge + nFace + nCell;
- // Alloc Temporary Space //
- const int orderPlus = order + 1;
- const int orderMinus = order - 1;
+ // Alloc Some Space //
+ const unsigned int orderPlus = order + 1;
+ const unsigned int orderMinus = order - 1;
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
Polynomial* u = new Polynomial[orderPlus];
Polynomial* v = new Polynomial[orderPlus];
- Polynomial lagrange [3];
- Polynomial lagrangeSum [3];
- Polynomial lagrangeSub [2][3];
-
// Legendre Polynomial //
Legendre::legendre(legendre, order);
Legendre::intScaled(intLegendre, orderPlus);
- // Vertices definig Edges & Permutations //
- const int edgeV[2][3][2] =
+ // Lagrange //
+ const Polynomial lagrange[3] =
{
- { {0, 1}, {1, 2}, {2, 0} },
- { {1, 0}, {2, 1}, {0, 2} }
- };
+ Polynomial(Polynomial(1, 0, 0, 0) -
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
- // 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);
+ Polynomial(Polynomial(1, 1, 0, 0)),
+ Polynomial(Polynomial(1, 0, 1, 0))
+ };
- // Lagrange //
- lagrange[0] =
- Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
-
- lagrange[1] =
- Polynomial(1, 1, 0, 0);
-
- lagrange[2] =
- Polynomial(1, 0, 1, 0);
-
- // Lagrange Sum //
- for(int e = 0; e < 3; e++)
- lagrangeSum[e] =
- lagrange[edgeV[0][e][0]] +
- lagrange[edgeV[0][e][1]];
-
- // Lagrange Sub //
- for(int e = 0; e < 3; e++){
- lagrangeSub[0][e] =
- lagrange[edgeV[0][e][0]] -
- lagrange[edgeV[0][e][1]];
-
- lagrangeSub[1][e] =
- lagrange[edgeV[1][e][0]] -
- lagrange[edgeV[1][e][1]];
- }
+ // Basis //
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
+
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
-
// Edge Based (Nedelec) //
- for(int c = 0; c < 2; c++){
- for(int e = 0; e < 3; 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]]);
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int e = 0; e < 3; e++){
+ vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
+
+ tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[2].mul(lagrange[(*(*edgeV[s])[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]);
- }
- }
+ tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+ (*(*basis)[s])[e] = new vector<Polynomial>(tmp2);
+ }
+ }
// Edge Based (High Order) //
- for(int c = 0; c < 2; c++){
- unsigned int i = 0;
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = 3;
- for(int l = 1; l < orderPlus; l++){
- for(int e = 0; e < 3; e++){
- (*(*edge)[c])[i + 3] =
+ for(unsigned int l = 1; l < orderPlus; l++){
+ for(unsigned int e = 0; e < 3; e++){
+ (*(*basis)[s])[i] =
new vector<Polynomial>
- ((intLegendre[l].compose(lagrangeSub[c][e],
- lagrangeSum[e])).gradient());
-
+ ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
+ lagrange[(*(*edgeV[s])[e])[1]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[1]] +
+ lagrange[(*(*edgeV[s])[e])[0]])).gradient());
i++;
}
}
}
+ // Cell Based //
- // Cell Based (Preliminary) //
- Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
+ // Preliminaries
+ const Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
- for(int l = 0; l < orderPlus; l++){
- u[l] = intLegendre[l].compose(lagrangeSub[0][0] * -1, lagrangeSum[0]);
+ for(unsigned int l = 0; l < orderPlus; l++){
+ u[l] = intLegendre[l].compose(lagrange[1] - lagrange[0],
+ lagrange[0] + lagrange[1]);
v[l] = legendre[l].compose(p);
v[l].mul(lagrange[2]);
}
@@ -149,65 +126,71 @@ TriEdgeBasis::TriEdgeBasis(int order){
subGradL1L2[1].sub(l1GradL2[1]);
subGradL1L2[2].sub(l1GradL2[2]);
- unsigned int i = 0;
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nEdge;
- // Cell Based (Type 1) //
- for(int l1 = 1; l1 < order; l1++){
- for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
- vector<Polynomial> tmp = v[l2].gradient();
- tmp[0].mul(u[l1]);
- tmp[1].mul(u[l1]);
- tmp[2].mul(u[l1]);
-
- (*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++;
+ // Type 1
+ for(unsigned int l1 = 1; l1 < order; l1++){
+ for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+ vector<Polynomial> tmp1 = v[l2].gradient();
+ vector<Polynomial> tmp2 = u[l1].gradient();
+
+ tmp1[0].mul(u[l1]);
+ tmp1[1].mul(u[l1]);
+ tmp1[2].mul(u[l1]);
+
+ tmp2[0].mul(v[l2]);
+ tmp2[1].mul(v[l2]);
+ tmp2[2].mul(v[l2]);
+
+ tmp2[0].add(tmp1[0]);
+ tmp2[1].add(tmp1[1]);
+ tmp2[2].add(tmp1[2]);
+
+ (*(*basis)[s])[i] = new vector<Polynomial>(tmp2);
+
+ i++;
+ }
}
- }
- // Cell Based (Type 2) //
- for(int l1 = 1; l1 < order; l1++){
- for(int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
- vector<Polynomial> tmp = v[l2].gradient();
- tmp[0].mul(u[l1]);
- tmp[1].mul(u[l1]);
- tmp[2].mul(u[l1]);
-
- (*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).sub(tmp[0]);
- (*cell)[i]->at(1).sub(tmp[1]);
- (*cell)[i]->at(2).sub(tmp[2]);
-
- i++;
+ // Type 2
+ for(unsigned int l1 = 1; l1 < order; l1++){
+ for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){
+ vector<Polynomial> tmp1 = v[l2].gradient();
+ vector<Polynomial> tmp2 = u[l1].gradient();
+
+ tmp1[0].mul(u[l1]);
+ tmp1[1].mul(u[l1]);
+ tmp1[2].mul(u[l1]);
+
+ tmp2[0].mul(v[l2]);
+ tmp2[1].mul(v[l2]);
+ tmp2[2].mul(v[l2]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
+ (*(*basis)[s])[i] = new vector<Polynomial>(tmp2);
+
+ i++;
+ }
}
- }
- // Cell Based (Type 3) //
- for(int l = 0; l < orderMinus; l++){
- vector<Polynomial> subGradL1L2V(subGradL1L2);
-
- subGradL1L2V[0].mul(v[l]);
- subGradL1L2V[1].mul(v[l]);
- subGradL1L2V[2].mul(v[l]);
+ // Type 3
+ for(unsigned int l = 0; l < orderMinus; l++){
+ vector<Polynomial> subGradL1L2V(subGradL1L2);
- (*cell)[i] = new vector<Polynomial>(subGradL1L2V);
-
- i++;
- }
-
+ subGradL1L2V[0].mul(v[l]);
+ subGradL1L2V[1].mul(v[l]);
+ subGradL1L2V[2].mul(v[l]);
+
+ (*(*basis)[s])[i] = new vector<Polynomial>(subGradL1L2V);
+
+ i++;
+ }
+ }
+
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
@@ -216,31 +199,16 @@ TriEdgeBasis::TriEdgeBasis(int order){
}
TriEdgeBasis::~TriEdgeBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
+ delete (*basis)[i];
}
-
- delete edge;
-
-
- // Face Based //
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/TriEdgeBasis.h b/FunctionSpace/TriEdgeBasis.h
index e8876d5f6ef57c5ac17ec9ad72edfaf65ae41129..4f9a1c40b74687147c46534e4baf837ec5415f26 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(int order);
+ TriEdgeBasis(unsigned int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriLagrangeBasis.cpp b/FunctionSpace/TriLagrangeBasis.cpp
index 093647adc921fc717869b2789fd0e5e2dd368e27..cfd97bb7ffa44c52c76863d46d74cad7171344a3 100644
--- a/FunctionSpace/TriLagrangeBasis.cpp
+++ b/FunctionSpace/TriLagrangeBasis.cpp
@@ -1,12 +1,15 @@
#include "BasisFactory.h"
#include "Exception.h"
+#include "TriReferenceSpace.h"
#include "TriLagrangeBasis.h"
using namespace std;
-TriLagrangeBasis::TriLagrangeBasis(int order){
+TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
+ throw Exception("TriLagrangeBasis: Removed");
+ /*
// Call polynomialBasis procedure //
- int tag;
+ unsigned int tag;
switch(order){
case 1:
@@ -60,60 +63,55 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
l = (polynomialBasis*)BasisFactory::create(tag);
point = new fullMatrix<double>(triPoint(order));
+ // Reference Space //
+ refSpace = new TriReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
// Set Basis Type //
this->order = order;
type = 0;
dim = 2;
- nVertex = 3;
- nEdge = 3 * (order - 1);
- nFace = 0;
- nCell = (order - 1) * (order - 2) / 2;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
-
- size = nVertex + nEdge + nFace + nCell;
-
+ nVertex = 3;
+ nEdge = 3 * (order - 1);
+ nFace = 0;
+ nCell = (order - 1) * (order - 2) / 2;
+ nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Some Stuff //
- const int nMonomial = l->monomials.size1();
+ const unsigned int nMonomial = l->monomials.size1();
unsigned int** edgeOrder;
Polynomial* pEdgeClosureLess = new Polynomial[nEdge];
-
-
+
// Basis //
- 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<Polynomial*>(nEdge);
- (*edge)[1] = new vector<Polynomial*>(nEdge);
+ basis = new vector<vector<const Polynomial*>*>(nRefSpace);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const Polynomial*>(nFunction);
// Vertex Based //
- for(int i = 0; i < nVertex; i++){
- Polynomial p = Polynomial(0, 0, 0, 0);
-
- for(int j = 0; j < nMonomial; j++)
- p = p + Polynomial(l->coefficients(i, j), // Coef
- l->monomials(j, 0), // powerX
- l->monomials(j, 1), // powerY
- 0); // powerZ
-
- (*node)[i] = new Polynomial(p);
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int i = 0; i < nVertex; i++){
+ Polynomial p = Polynomial(0, 0, 0, 0);
+
+ for(unsigned int j = 0; j < nMonomial; j++)
+ p = p + Polynomial(l->coefficients(i, j), // Coef
+ l->monomials(j, 0), // powerX
+ l->monomials(j, 1), // powerY
+ 0); // powerZ
+
+ (*(*basis)[s])[i] = new Polynomial(p);
+ }
}
-
// Edge Based //
// Without Closure
- for(int i = 0; i < nEdge; i++){
- int ci = i + nVertex;
+ for(unsigned int i = 0; i < nEdge; i++){
+ unsigned int ci = i + nVertex;
pEdgeClosureLess[i] = Polynomial(0, 0, 0, 0);
- for(int j = 0; j < nMonomial; j++)
+ for(unsigned int j = 0; j < nMonomial; j++)
pEdgeClosureLess[i] =
pEdgeClosureLess[i] +
Polynomial(l->coefficients(ci, j), // Coef
@@ -125,29 +123,27 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
// With Closure
edgeOrder = triEdgeOrder(order); // Closure Ordering
- for(int i = 0; i < nEdge; i++){
- (*(*edge)[0])[i] =
- new Polynomial(pEdgeClosureLess[edgeOrder[0][i]]);
-
- (*(*edge)[1])[i] =
- new Polynomial(pEdgeClosureLess[edgeOrder[1][i]]);
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int i = nVertex; i < nEdge; i++){
+ (*(*basis)[s])[i] =
+ new Polynomial(pEdgeClosureLess[edgeOrder[s][i]]);
+ }
}
-
// Cell Based //
- for(int i = 0; i < nCell; i++){
- int ci = i + nVertex + nEdge;
- Polynomial p = Polynomial(0, 0, 0, 0);
-
- for(int j = 0; j < nMonomial; j++)
- p = p + Polynomial(l->coefficients(ci, j), // Coef
- l->monomials(j, 0), // powerX
- l->monomials(j, 1), // powerY
- 0); // powerZ
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int i = nVertex + nEdge; i < nCell; i++){
+ Polynomial p = Polynomial(0, 0, 0, 0);
- (*cell)[i] = new Polynomial(p);
+ for(unsigned int j = 0; j < nMonomial; j++)
+ p = p + Polynomial(l->coefficients(i, j), // Coef
+ l->monomials(j, 0), // powerX
+ l->monomials(j, 1), // powerY
+ 0); // powerZ
+
+ (*(*basis)[s])[i] = new Polynomial(p);
+ }
}
-
// Delete Temporary Space //
delete[] pEdgeClosureLess;
@@ -157,47 +153,25 @@ TriLagrangeBasis::TriLagrangeBasis(int order){
delete[] edgeOrder[1];
delete[] edgeOrder;
}
+ */
}
TriLagrangeBasis::~TriLagrangeBasis(void){
// Delete gmsh polynomial Basis //
// --> no method to do so :-(
// --> erased ??
+ // ReferenceSpace //
+ delete refSpace;
- // 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;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Face Based //
- for(int c = 0; c < nFaceClosure; c++){
- for(int i = 0; i < nFace; i++)
- delete (*(*face)[c])[i];
-
- delete (*face)[c];
+ delete (*basis)[i];
}
- delete face;
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
-
- delete cell;
-
- // Delete Lagrange Points //
- delete point;
+ delete basis;
}
fullMatrix<double> TriLagrangeBasis::
diff --git a/FunctionSpace/TriLagrangeBasis.h b/FunctionSpace/TriLagrangeBasis.h
index c77554184cd231bcdded5a09569e4f9d7ae08e72..39b1b2f9765935084c10a7d0ed47b0eee730ebfc 100644
--- a/FunctionSpace/TriLagrangeBasis.h
+++ b/FunctionSpace/TriLagrangeBasis.h
@@ -23,7 +23,7 @@ class TriLagrangeBasis: public LagrangeBasis{
//!
//! Returns a new TriLagrangeBasis
//! of the given Order
- TriLagrangeBasis(int order);
+ TriLagrangeBasis(unsigned int order);
//! Deletes this Basis
//!
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 061e6161fcd76edb33a9b13826a9a1e258b27017..dab1f19c863fad9e873aaf76eb64f354e23e7ee0 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -1,103 +1,81 @@
#include "TriNedelecBasis.h"
+#include "TriReferenceSpace.h"
using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
+ // Reference Space //
+ refSpace = new TriReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
+ const vector<const vector<const vector<unsigned int>*>*>&
+ edgeV = refSpace->getAllEdge();
+
// Set Basis Type //
order = 1;
type = 1;
dim = 2;
- nVertex = 0;
- nEdge = 3;
- nFace = 0;
- nCell = 0;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
+ nVertex = 0;
+ nEdge = 3;
+ nFace = 0;
+ nCell = 0;
+ nFunction = 3;
- size = 3;
-
- // Vertices definig Edges & Permutations //
- const int edgeV[2][3][2] =
+ // Lagrange //
+ const Polynomial lagrange[3] =
{
- { {0, 1}, {1, 2}, {2, 0} },
- { {1, 0}, {2, 1}, {0, 2} }
- };
-
- // 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);
+ Polynomial(Polynomial(1, 0, 0, 0) -
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
+ Polynomial(Polynomial(1, 1, 0, 0)),
- // Lagrange //
- Polynomial lagrange[3];
+ Polynomial(Polynomial(1, 0, 1, 0))
+ };
- lagrange[0] =
- Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
-
- lagrange[1] =
- Polynomial(1, 1, 0, 0);
+ // Basis //
+ basis = new vector<vector<const vector<Polynomial>*>*>(nRefSpace);
- lagrange[2] =
- Polynomial(1, 0, 1, 0);
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const vector<Polynomial>*>(nFunction);
+ // Edge Based (Nedelec) //
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ for(unsigned int e = 0; e < 3; e++){
+ vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
+ vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
- // Nedelec //
- for(int c = 0; c < 2; c++){
- for(int e = 0; e < 3; 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]]);
+ tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
+ tmp1[2].mul(lagrange[(*(*edgeV[s])[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]);
+ tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+ tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
+
+ tmp2[0].sub(tmp1[0]);
+ tmp2[1].sub(tmp1[1]);
+ tmp2[2].sub(tmp1[2]);
+
+ (*(*basis)[s])[e] = new vector<Polynomial>(tmp2);
}
- }
+ }
}
TriNedelecBasis::~TriNedelecBasis(void){
- // 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;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
+ delete (*basis)[i];
+ }
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 8d0584ffa2493373310d2cf361fd0ffe905746af..0fab56d2c0cc526d6b5257c05700810af664e954 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -1,117 +1,99 @@
#include "TriNodeBasis.h"
+#include "TriReferenceSpace.h"
#include "Legendre.h"
using namespace std;
-TriNodeBasis::TriNodeBasis(int order){
- // Set Basis Type //
+TriNodeBasis::TriNodeBasis(unsigned int order){
+ // Reference Space //
+ refSpace = new TriReferenceSpace;
+ nRefSpace = refSpace->getNReferenceSpace();
+
+ const vector<const vector<const vector<unsigned int>*>*>&
+ edgeV = refSpace->getAllEdge();
+
+ // Set BasisTwo Type //
this->order = order;
type = 0;
dim = 2;
- nVertex = 3;
- nEdge = 3 * (order - 1);
- nFace = 0;
- nCell = (order - 1) * (order - 2) / 2;
-
- nEdgeClosure = 2;
- nFaceClosure = 0;
+ nVertex = 3;
+ nEdge = 3 * (order - 1);
+ nFace = 0;
+ nCell = (order - 1) * (order - 2) / 2;
+ nFunction = nVertex + nEdge + nFace + nCell;
- size = nVertex + nEdge + nFace + nCell;
-
- // Alloc Temporary Space //
- const int orderMinus = order - 1;
+ // Alloc Some Space //
+ const unsigned int orderMinus = order - 1;
Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
- Polynomial lagrangeSum[3];
- Polynomial lagrangeSub[2][3];
-
// Legendre Polynomial //
Legendre::legendre(legendre, orderMinus);
Legendre::intScaled(intLegendre, order);
- // Vertices definig Edges & Permutations //
- const int edgeV[2][3][2] =
+ // Lagrange Polynomial //
+ const Polynomial lagrange[3] =
{
- { {0, 1}, {1, 2}, {2, 0} },
- { {1, 0}, {2, 1}, {0, 2} }
- };
-
- // Basis //
- 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<Polynomial*>(nEdge);
- (*edge)[1] = new vector<Polynomial*>(nEdge);
-
+ Polynomial(Polynomial(1, 0, 0, 0) -
+ Polynomial(1, 1, 0, 0) -
+ Polynomial(1, 0, 1, 0)),
+
+ Polynomial(Polynomial(1, 1, 0, 0)),
- // Vertex Based (Lagrange) //
- (*node)[0] =
- new Polynomial(Polynomial(1, 0, 0, 0) -
- Polynomial(1, 1, 0, 0) -
- Polynomial(1, 0, 1, 0));
+ Polynomial(Polynomial(1, 0, 1, 0))
+ };
- (*node)[1] =
- new Polynomial(Polynomial(1, 1, 0, 0));
+ // Basis //
+ basis = new vector<vector<const Polynomial*>*>(nRefSpace);
- (*node)[2] =
- new Polynomial(Polynomial(1, 0, 1, 0));
-
+ for(unsigned int s = 0; s < nRefSpace; s++)
+ (*basis)[s] = new vector<const Polynomial*>(nFunction);
- // Lagrange Sum //
- for(int e = 0; e < 3; e++)
- lagrangeSum[e] =
- *(*node)[edgeV[0][e][0]] +
- *(*node)[edgeV[0][e][1]];
-
- // Lagrange Sub //
- for(int e = 0; e < 3; e++){
- lagrangeSub[0][e] =
- *(*node)[edgeV[0][e][1]] -
- *(*node)[edgeV[0][e][0]];
-
- lagrangeSub[1][e] =
- *(*node)[edgeV[1][e][1]] -
- *(*node)[edgeV[1][e][0]];
+ // Vertex Based //
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ (*(*basis)[s])[0] = new Polynomial(lagrange[0]);
+ (*(*basis)[s])[1] = new Polynomial(lagrange[1]);
+ (*(*basis)[s])[2] = new Polynomial(lagrange[2]);
}
-
-
+
// Edge Based //
- for(int c = 0; c < 2; c++){
- unsigned int i = 0;
-
- for(int l = 1; l < order; l++){
- for(int e = 0; e < 3; e++){
- (*(*edge)[c])[i] =
- new Polynomial(intLegendre[l].compose(lagrangeSub[c][e], lagrangeSum[e]));
-
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nVertex;
+
+ for(unsigned int l = 1; l < order; l++){
+ for(unsigned int e = 0; e < 3; e++){
+ (*(*basis)[s])[i] =
+ new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
+ lagrange[(*(*edgeV[s])[e])[0]]
+ ,
+ lagrange[(*(*edgeV[s])[e])[0]] +
+ lagrange[(*(*edgeV[s])[e])[1]]));
i++;
}
}
}
-
-
+
// Cell Based //
- Polynomial p = (*(*node)[2] * 2) - Polynomial(1, 0, 0, 0);
- const int orderMinusTwo = order - 2;
-
- unsigned int i = 0;
+ const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0);
+ const unsigned int orderMinusTwo = order - 2;
- for(int l1 = 1; l1 < orderMinus; l1++){
- for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
- (*cell)[i] =
- new Polynomial(intLegendre[l1].compose(lagrangeSub[0][0], lagrangeSum[0]) *
- legendre[l2].compose(p) * *(*node)[2]);
-
- i++;
+ for(unsigned int s = 0; s < nRefSpace; s++){
+ unsigned int i = nVertex + nEdge;
+
+ for(unsigned int l1 = 1; l1 < orderMinus; l1++){
+ for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
+ (*(*basis)[s])[i] =
+ new Polynomial(intLegendre[l1].compose(lagrange[1] - lagrange[0],
+ lagrange[1] + lagrange[0]) *
+ legendre[l2].compose(p) * lagrange[2]);
+
+ i++;
+ }
}
}
-
// Free Temporary Sapce //
delete[] legendre;
@@ -119,31 +101,16 @@ TriNodeBasis::TriNodeBasis(int order){
}
TriNodeBasis::~TriNodeBasis(void){
- // Vertex Based //
- for(int i = 0; i < nVertex; i++)
- delete (*node)[i];
-
- delete node;
+ // ReferenceSpace //
+ delete refSpace;
+ // Basis //
+ for(unsigned int i = 0; i < nRefSpace; i++){
+ for(unsigned int j = 0; j < nFunction; j++)
+ delete (*(*basis)[i])[j];
- // Edge Based //
- for(int c = 0; c < 2; c++){
- for(int i = 0; i < nEdge; i++)
- delete (*(*edge)[c])[i];
-
- delete (*edge)[c];
+ delete (*basis)[i];
}
-
- delete edge;
-
-
- // Face Based //
- delete face;
-
-
- // Cell Based //
- for(int i = 0; i < nCell; i++)
- delete (*cell)[i];
- delete cell;
+ delete basis;
}
diff --git a/FunctionSpace/TriNodeBasis.h b/FunctionSpace/TriNodeBasis.h
index 018b4174c28d15e83be43c5370e400c0ba2a0da7..215731cb51e0f129854f9b3d91023d44beb567f2 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(int order);
+ TriNodeBasis(unsigned int order);
//! Deletes this Basis
//!