diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 9d53d02e5f16d8b4b0bb2360c086793094c9204d..1742794385695c661c2ef87a7f681a20da664b06 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -25,11 +25,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial* legendreX = new Polynomial[orderPlus];
Polynomial* legendreY = new Polynomial[orderPlus];
- Polynomial* lagrange = new Polynomial[4];
- Polynomial* lagrangeSum = new Polynomial[4];
-
- Polynomial* lifting = new Polynomial[4];
- Polynomial* liftingSub = new Polynomial[4];
+ Polynomial lagrange[4];
+ Polynomial lifting[4];
+ Polynomial lagrangeSum[4];
+ Polynomial liftingSub[4];
+ Polynomial rLiftingSub[4];
// Integrated and classical Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus);
@@ -74,12 +74,14 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
(Polynomial(1, 0, 1, 0));
// Lifting Sub //
- for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
- liftingSub[i] = lifting[j] - lifting[i];
-
+ for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4){
+ liftingSub[i] = lifting[j] - lifting[i];
+ rLiftingSub[i] = lifting[i] - lifting[j];
+ }
// Basis (temporary --- *no* const) //
std::vector<std::vector<Polynomial>*> basis(size);
+ std::vector<std::vector<Polynomial>*> revBasis(size);
// Edge Based (Nedelec) //
@@ -87,6 +89,7 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial oneHalf(0.5, 0, 0, 0);
for(int e = 0; e < 4; e++){
+ // Direct
basis[i] = new std::vector<Polynomial>(liftingSub[e].gradient());
basis[i]->at(0).mul(lagrangeSum[e]);
@@ -97,6 +100,18 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
basis[i]->at(1).mul(oneHalf);
basis[i]->at(2).mul(oneHalf);
+ // Revert
+ revBasis[i] = new std::vector<Polynomial>(rLiftingSub[e].gradient());
+
+ revBasis[i]->at(0).mul(lagrangeSum[e]);
+ revBasis[i]->at(1).mul(lagrangeSum[e]);
+ revBasis[i]->at(2).mul(lagrangeSum[e]);
+
+ revBasis[i]->at(0).mul(oneHalf);
+ revBasis[i]->at(1).mul(oneHalf);
+ revBasis[i]->at(2).mul(oneHalf);
+
+ // Next
i++;
}
@@ -104,7 +119,12 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
for(int l = 1; l < orderPlus; l++){
for(int e = 0; e < 4; e++){
basis[i] =
- new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient());
+ new std::vector<Polynomial>((intLegendre[l].compose(liftingSub[e]) *
+ lagrangeSum[e]).gradient());
+
+ revBasis[i] =
+ new std::vector<Polynomial>((intLegendre[l].compose(rLiftingSub[e]) *
+ lagrangeSum[e]).gradient());
i++;
}
@@ -131,6 +151,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
for(int l2 = 1; l2 < orderPlus; l2++){
basis[i] = new std::vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
+ revBasis[i] = basis[i];
+
i++;
}
}
@@ -144,6 +166,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
basis[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
basis[i]->at(2) = zero;
+ revBasis[i] = basis[i];
+
i++;
}
}
@@ -161,9 +185,12 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
basis[iPlus]->at(1) = iLegendreX[l];
basis[iPlus]->at(2) = zero;
+ revBasis[i] = basis[i];
+
i++;
}
+
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;
@@ -173,21 +200,23 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
delete[] legendreX;
delete[] legendreY;
- delete[] lagrange;
- delete[] lagrangeSum;
-
- delete[] lifting;
- delete[] liftingSub;
-
// Set Basis //
this->basis = new std::vector<const std::vector<Polynomial>*>
(basis.begin(), basis.end());
+
+ this->revBasis = new std::vector<const std::vector<Polynomial>*>
+ (revBasis.begin(), revBasis.end());
}
QuadEdgeBasis::~QuadEdgeBasis(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/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 861ff3eae0e6a0815bdf1ce73e9e9dfbd20592be..72f9b7f74aa46b89a04f871259bc2396ba8aee09 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -17,7 +17,7 @@ QuadNodeBasis::QuadNodeBasis(const int order){
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
- Polynomial* lifting = new Polynomial[4];
+ Polynomial lifting[4];
// Legendre Polynomial //
Legendre::integrated(legendre, order);
@@ -42,7 +42,8 @@ QuadNodeBasis::QuadNodeBasis(const int order){
// Basis //
- basis = new std::vector<const Polynomial*>(size);
+ basis = new std::vector<const Polynomial*>(size);
+ revBasis = new std::vector<const Polynomial*>(size);
// Vertex Based (Lagrange) //
(*basis)[0] =
@@ -61,6 +62,11 @@ QuadNodeBasis::QuadNodeBasis(const int order){
new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)));
+ // Vertex Based (Revert) //
+ for(int i = 0; i < 3; i++)
+ (*revBasis)[i] = (*basis)[i];
+
+
// Edge Based //
int i = 4;
@@ -69,11 +75,15 @@ QuadNodeBasis::QuadNodeBasis(const int order){
(*basis)[i] =
new Polynomial(legendre[l].compose(lifting[e2] - lifting[e1]) *
(*(*basis)[e1] + *(*basis)[e2]));
-
+
+ (*revBasis)[i] =
+ new Polynomial(legendre[l].compose(lifting[e1] - lifting[e2]) *
+ (*(*basis)[e1] + *(*basis)[e2]));
i++;
}
}
+
// Cell Based //
Polynomial px = Polynomial(2, 1, 0, 0);
Polynomial py = Polynomial(2, 0, 1, 0);
@@ -86,18 +96,25 @@ QuadNodeBasis::QuadNodeBasis(const int order){
(*basis)[i] =
new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
+ (*revBasis)[i] = (*basis)[i];
+
i++;
}
}
+
// Free Temporary Sapce //
delete[] legendre;
- delete[] lifting;
}
QuadNodeBasis::~QuadNodeBasis(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/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index b6594bd3dead65fa5476af6668a1d96b6be202a4..2b0ef40f39017bc05f85d532dc2853c8d83ddc75 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -76,6 +76,7 @@ TriNodeBasis::TriNodeBasis(const int order){
}
}
+
// Cell Based //
Polynomial p = *(*basis)[2] * 2 - Polynomial(1, 0, 0, 0);
const int orderMinusTwo = order - 2;
@@ -87,10 +88,12 @@ TriNodeBasis::TriNodeBasis(const int order){
legendre[l2].compose(p) * *(*basis)[2]);
(*revBasis)[i] = (*basis)[i];
+
i++;
}
}
+
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;