diff --git a/FunctionSpace/BasisVector.h b/FunctionSpace/BasisVector.h
index 85e3491147a962917d4c2a281711087831f9d55b..6344f5229c6ea5dde4ae0617c137a587384f07c9 100644
--- a/FunctionSpace/BasisVector.h
+++ b/FunctionSpace/BasisVector.h
@@ -21,6 +21,7 @@
class BasisVector: public Basis{
protected:
std::vector<const std::vector<Polynomial>*>* basis;
+ std::vector<const std::vector<Polynomial>*>* revBasis;
public:
//! Deletes this BasisVector
@@ -31,6 +32,12 @@ class BasisVector: public Basis{
//! defining this (vectorial) Basis
const std::vector<const std::vector<Polynomial>*>& getFunctions(void) const;
+ //! @param closure A natural number
+ //! @return Returns the set of @em Polynomial%s
+ //! defining this (scalar) Basis, for the given closure
+ const std::vector<const std::vector<Polynomial>*>&
+ getFunctions(unsigned int closure) const;
+
protected:
//! @internal
//! Instantiate a new BasisVector
@@ -50,4 +57,15 @@ getFunctions(void) const{
return *basis;
}
+inline
+const std::vector<const std::vector<Polynomial>*>& BasisVector::
+getFunctions(unsigned int closure) const{
+ if(!closure)
+ return *basis;
+
+ else
+ return *revBasis;
+}
+
+
#endif
diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index 6673e9f467e695c4d4d7373c6e057c273919ec6b..42effb1d876ab00ad7a15613692c1a66cdedd76f 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -25,7 +25,7 @@ interpolate(const MElement& element,
const_cast<MElement&>(element);
// Get Basis Functions //
- const vector<const vector<Polynomial>*>& fun = getBasis(element).getFunctions();
+ const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
const unsigned int nFun = fun.size();
// Get Reference coordinate //
diff --git a/FunctionSpace/FunctionSpaceVector.cpp b/FunctionSpace/FunctionSpaceVector.cpp
index 821e8505ce8dd95d94e9b51cb42f19143048dfa9..e34b0ba390bcf57930a8d8ae1da7999767ed2a63 100644
--- a/FunctionSpace/FunctionSpaceVector.cpp
+++ b/FunctionSpace/FunctionSpaceVector.cpp
@@ -1,4 +1,41 @@
#include "FunctionSpaceVector.h"
+#include "Exception.h"
+
+using namespace std;
FunctionSpaceVector::~FunctionSpaceVector(void){
+}
+
+ const vector<const vector<Polynomial>*> FunctionSpaceVector::
+getLocalFunctions(const MElement& element) const{
+
+ // Get Closure //
+ map<const MElement*, vector<bool>*>::iterator it =
+ edgeClosure->find(&element);
+
+ if(it == edgeClosure->end())
+ throw Exception("Element not found for closure");
+
+ vector<bool>* closure = it->second;
+ const unsigned int size = closure->size();
+
+ // Get Basis //
+ const vector<const vector<Polynomial>*>& basis =
+ getBasis(element).getFunctions(0);
+
+ const vector<const vector<Polynomial>*>& revBasis =
+ getBasis(element).getFunctions(1);
+
+ // Get Functions //
+ vector<const vector<Polynomial>*> fun(size);
+
+ for(unsigned int i = 0; i < size; i++)
+ if((*closure)[i])
+ fun[i] = basis[i];
+
+ else
+ fun[i] = revBasis[i];
+
+ // Return
+ return fun;
}
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index edb168f1ad6590bc6e408130a3832444714264b0..8c7f5fb3ae4d62b30a8bc6fdc490212403b32bf7 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -29,6 +29,9 @@ class FunctionSpaceVector : public FunctionSpace{
const std::vector<double>& coef,
const fullVector<double>& xyz) const = 0;
+ const std::vector<const std::vector<Polynomial>*>
+ getLocalFunctions(const MElement& element) const;
+
const BasisVector& getBasis(const MElement& element) const;
};
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 8b70bc1fd2cc2f62d49038fb9c25b9244327cf0b..884b04531a815fae9dd3b9c877a5e42fdf50ed47 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -21,14 +21,14 @@ TriEdgeBasis::TriEdgeBasis(const int order){
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
-
- Polynomial* lagrange = new Polynomial[3];
- Polynomial* lagrangeSub = new Polynomial[3];
- Polynomial* lagrangeSum = new Polynomial[3];
-
Polynomial* u = new Polynomial[orderPlus];
Polynomial* v = new Polynomial[orderPlus];
+ Polynomial lagrange [3];
+ Polynomial lagrangeSum [3];
+ Polynomial lagrangeSub [3];
+ Polynomial rLagrangeSub [3];
+
// Classical and Intrated-Scaled Legendre Polynomial //
Legendre::legendre(legendre, order);
Legendre::intScaled(intLegendre, orderPlus);
@@ -45,21 +45,25 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Lagrange Sum //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = lagrange[j] + lagrange[i];
+ lagrangeSum[i] = lagrange[j] + lagrange[i];
- // Lagrange Sub //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSub[i] = lagrange[i] - lagrange[j];
+ // Lagrange Sub (& revert) //
+ for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
+ lagrangeSub[i] = lagrange[i] - lagrange[j];
+ rLagrangeSub[i] = lagrange[j] - lagrange[i];
+ }
- // Basis (temporary --- *no* const) //
- std::vector<std::vector<Polynomial>*> basis(size);
+ // Basis & Revert (temporary --- *no* const) //
+ std::vector<std::vector<Polynomial>*> basis(size);
+ std::vector<std::vector<Polynomial>*> revBasis(size);
// Edge Based (Nedelec) //
int i = 0;
for(int j = 1; i < 3; j = (j + 1) % 3){
+ // Direct
std::vector<Polynomial> tmp = lagrange[j].gradient();
tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]);
@@ -75,6 +79,23 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]);
+ // Revert
+ std::vector<Polynomial> tmpR = lagrange[i].gradient();
+ tmpR[0].mul(lagrange[j]);
+ tmpR[1].mul(lagrange[j]);
+ tmpR[2].mul(lagrange[j]);
+
+ revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
+
+ revBasis[i]->at(0).mul(lagrange[i]);
+ revBasis[i]->at(1).mul(lagrange[i]);
+ revBasis[i]->at(2).mul(lagrange[i]);
+
+ revBasis[i]->at(0).sub(tmpR[0]);
+ revBasis[i]->at(1).sub(tmpR[1]);
+ revBasis[i]->at(2).sub(tmpR[2]);
+
+ // Next
i++;
}
@@ -82,8 +103,12 @@ TriEdgeBasis::TriEdgeBasis(const int order){
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 3; e++){
basis[i] =
- new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e])).gradient());
+ new std::vector<Polynomial>((intLegendre[l].compose(lagrangeSub[e],
+ lagrangeSum[e])).gradient());
+ revBasis[i] =
+ new std::vector<Polynomial>((intLegendre[l].compose(rLagrangeSub[e],
+ lagrangeSum[e])).gradient());
i++;
}
}
@@ -114,6 +139,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(0).add(tmp[0]);
basis[i]->at(1).add(tmp[1]);
basis[i]->at(2).add(tmp[2]);
+
+ revBasis[i] = basis[i];
i++;
}
@@ -137,6 +164,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]);
+ revBasis[i] = basis[i];
+
i++;
}
}
@@ -148,6 +177,8 @@ TriEdgeBasis::TriEdgeBasis(const int order){
basis[i]->at(0).mul(v[l]);
basis[i]->at(1).mul(v[l]);
basis[i]->at(2).mul(v[l]);
+
+ revBasis[i] = basis[i];
i++;
}
@@ -155,21 +186,26 @@ TriEdgeBasis::TriEdgeBasis(const int order){
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
- delete[] lagrange;
- delete[] lagrangeSub;
- delete[] lagrangeSum;
delete[] u;
delete[] v;
// Set Basis //
- this->basis = new std::vector<const std::vector<Polynomial>*>
+ this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end());
+
+ this->revBasis = new std::vector<const std::vector<Polynomial>*>
+ (revBasis.begin(), revBasis.end());
}
TriEdgeBasis::~TriEdgeBasis(void){
- for(int i = 0; i < size; i++)
+ for(int i = 0; i < size; i++){
delete (*basis)[i];
+ if(i >= nVertex && i < nVertex + nEdge)
+ delete (*revBasis)[i];
+ }
+
delete basis;
+ delete revBasis;
}
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index 7a8119120065834a17480d9a4b5e78db44506736..dc3cae96106cb2a01b6812059fc3d9c219b33d48 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -15,7 +15,7 @@ TriNedelecBasis::TriNedelecBasis(void){
size = 3;
// Lagrange //
- Polynomial* lagrange = new Polynomial[3];
+ Polynomial lagrange[3];
lagrange[0] =
Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
@@ -28,10 +28,12 @@ TriNedelecBasis::TriNedelecBasis(void){
// Basis (temporary --- *no* const) //
- std::vector<std::vector<Polynomial>*> basis(size);
+ std::vector<std::vector<Polynomial>*> basis(size);
+ std::vector<std::vector<Polynomial>*> revBasis(size);
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
+ // Direct
std::vector<Polynomial> tmp = lagrange[j].gradient();
tmp[0].mul(lagrange[i]);
tmp[1].mul(lagrange[i]);
@@ -46,15 +48,31 @@ TriNedelecBasis::TriNedelecBasis(void){
basis[i]->at(0).sub(tmp[0]);
basis[i]->at(1).sub(tmp[1]);
basis[i]->at(2).sub(tmp[2]);
- }
- // Free Temporary Sapce //
- delete[] lagrange;
+ // Revert
+ std::vector<Polynomial> tmpR = lagrange[i].gradient();
+ tmpR[0].mul(lagrange[j]);
+ tmpR[1].mul(lagrange[j]);
+ tmpR[2].mul(lagrange[j]);
+
+ revBasis[i] = new std::vector<Polynomial>(lagrange[j].gradient());
+
+ revBasis[i]->at(0).mul(lagrange[i]);
+ revBasis[i]->at(1).mul(lagrange[i]);
+ revBasis[i]->at(2).mul(lagrange[i]);
+
+ revBasis[i]->at(0).sub(tmpR[0]);
+ revBasis[i]->at(1).sub(tmpR[1]);
+ revBasis[i]->at(2).sub(tmpR[2]);
+ }
// Set Basis //
- this->basis = new std::vector<const std::vector<Polynomial>*>
+ this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end());
+
+ this->revBasis = new std::vector<const std::vector<Polynomial>*>
+ (revBasis.begin(), revBasis.end());
}
TriNedelecBasis::~TriNedelecBasis(void){
@@ -62,4 +80,9 @@ TriNedelecBasis::~TriNedelecBasis(void){
delete (*basis)[i];
delete basis;
+
+ for(int i = 0; i < size; i++)
+ delete (*revBasis)[i];
+
+ delete revBasis;
}
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 0be867a6efc811b45edd0e8f0a33e3c7760b4bc2..b6594bd3dead65fa5476af6668a1d96b6be202a4 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -19,10 +19,9 @@ TriNodeBasis::TriNodeBasis(const int order){
Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order];
- Polynomial lagrangeSub[3];
Polynomial lagrangeSum[3];
+ Polynomial lagrangeSub[3];
Polynomial rLagrangeSub[3];
- Polynomial rLagrangeSum[3];
// Classical and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1;
@@ -51,13 +50,9 @@ TriNodeBasis::TriNodeBasis(const int order){
for(int i = 0; i < 3; i++)
(*revBasis)[i] = (*basis)[i];
-
// Lagrange Sum (& revert) //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
+ for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
lagrangeSum[i] = *(*basis)[i] + *(*basis)[j];
- rLagrangeSum[i] = *(*basis)[j] + *(*basis)[i];
- }
-
// Lagrange Sub (& revert) //
for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
@@ -75,7 +70,7 @@ TriNodeBasis::TriNodeBasis(const int order){
intLegendre[l].compose(lagrangeSub[e], lagrangeSum[e]));
(*revBasis)[i] = new Polynomial(
- intLegendre[l].compose(rLagrangeSub[e], rLagrangeSum[e]));
+ intLegendre[l].compose(rLagrangeSub[e], lagrangeSum[e]));
i++;
}