diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 43ca97949e97bc996a93d64022b1b3aa595ecd21..05bb6fc614da6f776511c50d25f5bcaca0284085 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -1,6 +1,8 @@
#include "HexNodeBasis.h"
#include "Legendre.h"
+using namespace std;
+
HexNodeBasis::HexNodeBasis(const int order){
// Set Basis Type //
this->order = order;
@@ -13,7 +15,7 @@ HexNodeBasis::HexNodeBasis(const int order){
nFace = 6 * (order - 1) * (order - 1);
nCell = (order - 1) * (order - 1) * (order - 1);
- size = (order + 1) * (order + 1) * (order + 1);
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
@@ -25,6 +27,27 @@ HexNodeBasis::HexNodeBasis(const int order){
// Legendre Polynomial //
Legendre::integrated(legendre, order);
+
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][12][2] =
+ {
+ { {0, 1}, {0, 3}, {0, 4}, {1, 2}, {1, 5}, {2, 3},
+ {2, 6}, {3, 7}, {4, 5}, {4, 7}, {5, 6}, {6, 7} },
+
+ { {1, 0}, {3, 0}, {4, 0}, {2, 1}, {5, 1}, {3, 2},
+ {6, 2}, {7, 3}, {5, 4}, {7, 4}, {6, 5}, {7, 6} },
+ };
+
+ const int faceV[6][4] =
+ {
+ {0, 3, 2, 1},
+ {0, 1, 5, 4},
+ {0, 4, 7, 3},
+ {1, 2, 6, 5},
+ {2, 3, 7, 6},
+ {4, 5, 6, 7}
+ };
+
// Lifting //
lifting[0] =
@@ -67,7 +90,7 @@ HexNodeBasis::HexNodeBasis(const int order){
(Polynomial(1, 0, 1, 0)) +
Polynomial(1, 0, 0, 1);
-
+ /*
// Basis //
basis = new std::vector<const Polynomial*>(size);
@@ -194,7 +217,7 @@ HexNodeBasis::HexNodeBasis(const int order){
}
}
-
+ */
// Free Temporary Sapce //
delete[] legendre;
delete[] lifting;
@@ -205,8 +228,38 @@ HexNodeBasis::HexNodeBasis(const int order){
}
HexNodeBasis::~HexNodeBasis(void){
- for(int i = 0; i < size; i++)
- delete (*basis)[i];
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
+
+
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
+ }
+
+ delete edge;
+
+
+ // Face Based //
+ for(int c = 0; c < 6; 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 basis;
+ delete cell;
}
diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 1742794385695c661c2ef87a7f681a20da664b06..f14e14c67b4b09b49fd5c272325f99402052c516 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -1,6 +1,8 @@
#include "QuadEdgeBasis.h"
#include "Legendre.h"
+using namespace std;
+
QuadEdgeBasis::QuadEdgeBasis(const int order){
// Set Basis Type //
this->order = order;
@@ -8,15 +10,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
type = 1;
dim = 2;
- nVertex = 0 ;
- nEdge = 4 * (order + 1) ;
- nFace = 0 ;
+ nVertex = 0;
+ nEdge = 4 * (order + 1);
+ nFace = 0;
nCell = 2 * (order + 1) * order;
- size = 2 * (order + 2) * (order + 1);
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
+
Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus];
@@ -28,13 +31,29 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
Polynomial lagrange[4];
Polynomial lifting[4];
Polynomial lagrangeSum[4];
- Polynomial liftingSub[4];
- Polynomial rLiftingSub[4];
+ Polynomial liftingSub[2][4];
- // Integrated and classical Legendre Polynomial //
+ // Legendre Polynomial //
Legendre::integrated(intLegendre, orderPlus);
Legendre::legendre(legendre, order);
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][4][2] =
+ {
+ { {0, 1}, {1, 2}, {2, 3}, {3, 0} },
+ { {1, 0}, {2, 1}, {3, 2}, {0, 3} }
+ };
+
+ // 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);
+
+
// Lagrange //
lagrange[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
@@ -53,8 +72,10 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
(Polynomial(1, 0, 1, 0));
// Lagrange Sum //
- for(int i = 0, j = 1; i < 4; i++, j = (j + 1) % 4)
- lagrangeSum[i] = lagrange[i] + lagrange[j];
+ for(int e = 0; e < 4; e++)
+ lagrangeSum[e] =
+ lagrange[edgeV[0][e][0]] +
+ lagrange[edgeV[0][e][1]];
// Lifting //
lifting[0] =
@@ -74,59 +95,48 @@ 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];
- rLiftingSub[i] = lifting[i] - lifting[j];
- }
+ for(int e = 0; e < 4; e++){
+ liftingSub[0][e] =
+ lifting[edgeV[0][e][0]] -
+ lifting[edgeV[0][e][1]];
- // Basis (temporary --- *no* const) //
- std::vector<std::vector<Polynomial>*> basis(size);
- std::vector<std::vector<Polynomial>*> revBasis(size);
+ liftingSub[1][e] =
+ lifting[edgeV[1][e][0]] -
+ lifting[edgeV[1][e][1]];
+ }
// Edge Based (Nedelec) //
- int i = 0;
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]);
- basis[i]->at(1).mul(lagrangeSum[e]);
- basis[i]->at(2).mul(lagrangeSum[e]);
-
- basis[i]->at(0).mul(oneHalf);
- basis[i]->at(1).mul(oneHalf);
- basis[i]->at(2).mul(oneHalf);
-
- // Revert
- revBasis[i] = new std::vector<Polynomial>(rLiftingSub[e].gradient());
+ for(int c = 0; c < 2; c++){
+ for(int e = 0; e < 4; e++){
+ (*(*edge)[c])[e] =
+ new vector<Polynomial>(liftingSub[c][e].gradient());
- revBasis[i]->at(0).mul(lagrangeSum[e]);
- revBasis[i]->at(1).mul(lagrangeSum[e]);
- revBasis[i]->at(2).mul(lagrangeSum[e]);
+ (*(*edge)[c])[e]->at(0).mul(lagrangeSum[e]);
+ (*(*edge)[c])[e]->at(1).mul(lagrangeSum[e]);
+ (*(*edge)[c])[e]->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++;
+ (*(*edge)[c])[e]->at(0).mul(oneHalf);
+ (*(*edge)[c])[e]->at(1).mul(oneHalf);
+ (*(*edge)[c])[e]->at(2).mul(oneHalf);
+ }
}
+
// Edge Based (High 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());
-
- revBasis[i] =
- new std::vector<Polynomial>((intLegendre[l].compose(rLiftingSub[e]) *
- lagrangeSum[e]).gradient());
-
- i++;
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
+
+ for(int l = 1; l < orderPlus; l++){
+ for(int e = 0; e < 4; e++){
+ (*(*edge)[c])[i + 4] =
+ new vector<Polynomial>((intLegendre[l].compose(liftingSub[c][e]) *
+ lagrangeSum[e]).gradient());
+
+ i++;
+ }
}
}
@@ -139,6 +149,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
px = px - Polynomial(1, 0, 0, 0);
py = py - Polynomial(1, 0, 0, 0);
+ unsigned int i = 0;
+
for(int l = 0; l < orderPlus; l++){
iLegendreX[l] = intLegendre[l].compose(px);
iLegendreY[l] = intLegendre[l].compose(py);
@@ -149,9 +161,8 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 1) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
- basis[i] = new std::vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
-
- revBasis[i] = basis[i];
+ (*cell)[i] =
+ new vector<Polynomial>((iLegendreX[l1] * iLegendreY[l2]).gradient());
i++;
}
@@ -160,13 +171,11 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 2) //
for(int l1 = 1; l1 < orderPlus; l1++){
for(int l2 = 1; l2 < orderPlus; l2++){
- basis[i] = new std::vector<Polynomial>(3);
-
- basis[i]->at(0) = legendreX[l1] * iLegendreY[l2];
- basis[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
- basis[i]->at(2) = zero;
+ (*cell)[i] = new vector<Polynomial>(3);
- revBasis[i] = basis[i];
+ (*cell)[i]->at(0) = legendreX[l1] * iLegendreY[l2];
+ (*cell)[i]->at(1) = iLegendreX[l1] * legendreY[l2] * -1;
+ (*cell)[i]->at(2) = zero;
i++;
}
@@ -174,18 +183,16 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
// Cell Based (Type 3) //
for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
- basis[i] = new std::vector<Polynomial>(3);
- basis[iPlus] = new std::vector<Polynomial>(3);
+ (*cell)[i] = new vector<Polynomial>(3);
+ (*cell)[iPlus] = new vector<Polynomial>(3);
- basis[i]->at(0) = iLegendreY[l];
- basis[i]->at(1) = zero;
- basis[i]->at(2) = zero;
+ (*cell)[i]->at(0) = iLegendreY[l];
+ (*cell)[i]->at(1) = zero;
+ (*cell)[i]->at(2) = zero;
- basis[iPlus]->at(0) = zero;
- basis[iPlus]->at(1) = iLegendreX[l];
- basis[iPlus]->at(2) = zero;
-
- revBasis[i] = basis[i];
+ (*cell)[iPlus]->at(0) = zero;
+ (*cell)[iPlus]->at(1) = iLegendreX[l];
+ (*cell)[iPlus]->at(2) = zero;
i++;
}
@@ -199,24 +206,34 @@ QuadEdgeBasis::QuadEdgeBasis(const int order){
delete[] iLegendreY;
delete[] legendreX;
delete[] legendreY;
-
-
- // 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++){
- delete (*basis)[i];
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
+
- if(i >= nVertex && i < nVertex + nEdge)
- delete (*revBasis)[i];
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
}
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
+
+
+ // Cell Based //
+ for(int i = 0; i < nCell; i++)
+ delete (*cell)[i];
- delete basis;
- delete revBasis;
+ delete cell;
}
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 72f9b7f74aa46b89a04f871259bc2396ba8aee09..687d3edb771ec4c526798477b5a60a35be637f30 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -1,6 +1,8 @@
#include "QuadNodeBasis.h"
#include "Legendre.h"
+using namespace std;
+
QuadNodeBasis::QuadNodeBasis(const int order){
// Set Basis Type //
this->order = order;
@@ -8,20 +10,38 @@ QuadNodeBasis::QuadNodeBasis(const int order){
type = 0;
dim = 2;
- nVertex = 4 ;
- nEdge = 4 * (order - 1) ;
- nFace = 0 ;
+ nVertex = 4;
+ nEdge = 4 * (order - 1);
+ nFace = 0;
nCell = (order - 1) * (order - 1);
- size = (order + 1) * (order + 1);
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
Polynomial* legendre = new Polynomial[order];
Polynomial lifting[4];
+ Polynomial liftingSub[2][4];
// Legendre Polynomial //
Legendre::integrated(legendre, order);
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][4][2] =
+ {
+ { {0, 1}, {1, 2}, {2, 3}, {3, 0} },
+ { {1, 0}, {2, 1}, {3, 2}, {0, 3} }
+ };
+
+ // 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);
+
+
// Lifting //
lifting[0] =
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
@@ -39,47 +59,48 @@ QuadNodeBasis::QuadNodeBasis(const int order){
(Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
(Polynomial(1, 0, 1, 0));
+ // Lifting Sub //
+ for(int e = 0; e < 4; e++){
+ liftingSub[0][e] =
+ lifting[edgeV[0][e][0]] -
+ lifting[edgeV[0][e][1]];
+
+ liftingSub[1][e] =
+ lifting[edgeV[1][e][0]] -
+ lifting[edgeV[1][e][1]];
+ }
- // Basis //
- basis = new std::vector<const Polynomial*>(size);
- revBasis = new std::vector<const Polynomial*>(size);
-
// Vertex Based (Lagrange) //
- (*basis)[0] =
+ (*node)[0] =
new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
- (*basis)[1] =
+ (*node)[1] =
new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
- (*basis)[2] =
+ (*node)[2] =
new Polynomial((Polynomial(1, 1, 0, 0)) *
(Polynomial(1, 0, 1, 0)));
- (*basis)[3] =
+ (*node)[3] =
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;
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
- for(int l = 1; l < order; l++){
- for(int e1 = 0, e2 = 1; e1 < 4; e1++, e2 = (e2 + 1) % 4){
- (*basis)[i] =
- new Polynomial(legendre[l].compose(lifting[e2] - lifting[e1]) *
- (*(*basis)[e1] + *(*basis)[e2]));
+ for(int l = 1; l < order; l++){
+ for(int e = 0; e < 4; e++){
+ (*(*edge)[c])[i] =
+ new Polynomial(legendre[l].compose(liftingSub[c][e]) *
+ (*(*node)[edgeV[c][e][0]] + *(*node)[edgeV[c][e][1]]));
- (*revBasis)[i] =
- new Polynomial(legendre[l].compose(lifting[e1] - lifting[e2]) *
- (*(*basis)[e1] + *(*basis)[e2]));
- i++;
+ i++;
+ }
}
}
@@ -91,13 +112,13 @@ QuadNodeBasis::QuadNodeBasis(const int order){
px = px - Polynomial(1, 0, 0, 0);
py = py - Polynomial(1, 0, 0, 0);
+ unsigned int i = 0;
+
for(int l1 = 1; l1 < order; l1++){
for(int l2 = 1; l2 < order; l2++){
- (*basis)[i] =
+ (*cell)[i] =
new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
- (*revBasis)[i] = (*basis)[i];
-
i++;
}
}
@@ -108,13 +129,31 @@ QuadNodeBasis::QuadNodeBasis(const int order){
}
QuadNodeBasis::~QuadNodeBasis(void){
- for(int i = 0; i < size; i++){
- delete (*basis)[i];
+ // Vertex Based //
+ for(int i = 0; i < nVertex; i++)
+ delete (*node)[i];
+
+ delete node;
- if(i >= nVertex && i < nVertex + nEdge)
- delete (*revBasis)[i];
+
+ // Edge Based //
+ for(int c = 0; c < 2; c++){
+ for(int i = 0; i < nEdge; i++)
+ delete (*(*edge)[c])[i];
+
+ delete (*edge)[c];
}
+
+ delete edge;
+
+
+ // Face Based //
+ delete face;
+
+
+ // Cell Based //
+ for(int i = 0; i < nCell; i++)
+ delete (*cell)[i];
- delete basis;
- delete revBasis;
+ delete cell;
}
diff --git a/FunctionSpace/TetEdgeBasis.cpp b/FunctionSpace/TetEdgeBasis.cpp
index 757a05cc5f338330b24d9d3a1547ba022cf01ebe..705f36a783ef5463e4ce7d402a84975224ea8ef2 100644
--- a/FunctionSpace/TetEdgeBasis.cpp
+++ b/FunctionSpace/TetEdgeBasis.cpp
@@ -28,7 +28,7 @@ TetEdgeBasis::TetEdgeBasis(int order){
Polynomial* sclLegendre = new Polynomial[orderMinus];
Polynomial* intLegendre = new Polynomial[orderPlus];
- // Classical and Intrated-Scaled Legendre Polynomial //
+ // Legendre Polynomial //
Legendre::legendre(legendre, orderMinusTwo);
Legendre::scaled(sclLegendre, orderMinusTwo);
Legendre::intScaled(intLegendre, orderPlus);
@@ -97,7 +97,8 @@ TetEdgeBasis::TetEdgeBasis(int order){
tmp[1].mul(lagrange[edgeV[c][e][0]]);
tmp[2].mul(lagrange[edgeV[c][e][0]]);
- (*(*edge)[c])[e] = new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
+ (*(*edge)[c])[e] =
+ 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]]);
diff --git a/FunctionSpace/TetNodeBasis.cpp b/FunctionSpace/TetNodeBasis.cpp
index 602900f5a6b9c830f5d5e61d7e94a20823f91e84..c608aa794907fe86b784d7a2de386493fdcf0d57 100644
--- a/FunctionSpace/TetNodeBasis.cpp
+++ b/FunctionSpace/TetNodeBasis.cpp
@@ -18,39 +18,39 @@ TetNodeBasis::TetNodeBasis(int order){
size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
- Polynomial* legendre = new Polynomial[order];
- Polynomial* sclLegendre = new Polynomial[order];
- Polynomial* intLegendre = new Polynomial[order];
-
- // Classical, Scaled and Intrated-Scaled Legendre Polynomial //
const int orderMinus = order - 1;
const int orderMinusTwo = order - 2;
const int orderMinusThree = order - 3;
+ Polynomial* legendre = new Polynomial[order];
+ Polynomial* sclLegendre = new Polynomial[order];
+ Polynomial* intLegendre = new Polynomial[order];
+
+ // Legendre Polynomial //
Legendre::legendre(legendre, orderMinus);
Legendre::scaled(sclLegendre, orderMinus);
Legendre::intScaled(intLegendre, order);
- // Vertices definig Edges //
- const int edgeV[6][2] = {
- {0, 1},
- {1, 2},
- {2, 0},
- {3, 0},
- {3, 2},
- {3, 1}
- };
-
- // Vertices definig Faces //
- const int faceV[4][3] = {
- {0, 2, 1},
- {0, 1, 3},
- {0, 3, 2},
- {3, 1, 2}
- };
-
- // Counter //
- int i;
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][6][2] =
+ {
+ { {0, 1}, {1, 2}, {2, 0},
+ {3, 0}, {3, 2}, {3, 1} },
+
+ { {1, 0}, {2, 1}, {0, 2},
+ {0, 3}, {2, 3}, {1, 3} }
+ };
+
+ // Vertices definig Faces & Permutations //
+ const int faceV[6][4][3] =
+ {
+ { {0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2} },
+ { {2, 0, 1}, {1, 0, 3}, {3, 0, 2}, {1, 3, 2} },
+ { {2, 1, 0}, {1, 3, 0}, {3, 2, 0}, {1, 2, 3} },
+ { {1, 2, 0}, {3, 1, 0}, {2, 3, 0}, {2, 1, 3} },
+ { {1, 0, 2}, {3, 0, 1}, {2, 0, 3}, {2, 3, 1} },
+ { {0, 1, 2}, {0, 3, 1}, {0, 2, 3}, {3, 2, 1} },
+ };
// Basis //
node = new vector<Polynomial*>(nVertex);
@@ -87,97 +87,46 @@ TetNodeBasis::TetNodeBasis(int order){
// Edge Based //
- i = 0;
-
- for(int l = 1; l < order; l++){
- for(int e = 0; e < 6; e++){
- (*(*edge)[0])[i] =
- new Polynomial(intLegendre[l].compose
- (*(*node)[edgeV[e][0]] - *(*node)[edgeV[e][1]],
- *(*node)[edgeV[e][0]] + *(*node)[edgeV[e][1]]));
-
- (*(*edge)[1])[i] =
- new Polynomial(intLegendre[l].compose
- (*(*node)[edgeV[e][1]] - *(*node)[edgeV[e][0]],
- *(*node)[edgeV[e][1]] + *(*node)[edgeV[e][0]]));
-
- i++;
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
+
+ for(int l = 1; l < order; l++){
+ for(int e = 0; e < 6; e++){
+ (*(*edge)[c])[i] =
+ new Polynomial(intLegendre[l].compose
+ (*(*node)[edgeV[c][e][0]] - *(*node)[edgeV[c][e][1]],
+ *(*node)[edgeV[c][e][0]] + *(*node)[edgeV[c][e][1]]));
+
+ i++;
+ }
}
}
// Face Based //
- i = 0;
-
- for(int l1 = 1; l1 < orderMinus; l1++){
- for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
- for(int m = 0; m < 4; m++){
- Polynomial sum =
- *(*node)[faceV[m][0]] +
- *(*node)[faceV[m][1]] +
- *(*node)[faceV[m][2]];
-
- (*(*face)[0])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][0]] - *(*node)[faceV[m][1]],
- *(*node)[faceV[m][0]] + *(*node)[faceV[m][1]]) *
-
- *(*node)[faceV[m][2]] *
-
- sclLegendre[l2].compose
- (*(*node)[faceV[m][2]] * 2 - sum, sum));
-
- (*(*face)[1])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][1]] - *(*node)[faceV[m][0]],
- *(*node)[faceV[m][1]] + *(*node)[faceV[m][0]]) *
-
- *(*node)[faceV[m][2]] *
-
- sclLegendre[l2].compose
- (*(*node)[faceV[m][2]] * 2 - sum, sum));
-
- (*(*face)[2])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][1]] - *(*node)[faceV[m][2]],
- *(*node)[faceV[m][1]] + *(*node)[faceV[m][2]]) *
-
- *(*node)[faceV[m][0]] *
-
- sclLegendre[l2].compose
- (*(*node)[faceV[m][0]] * 2 - sum, sum));
-
- (*(*face)[3])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][2]] - *(*node)[faceV[m][1]],
- *(*node)[faceV[m][2]] + *(*node)[faceV[m][1]]) *
-
- *(*node)[faceV[m][0]] *
-
- sclLegendre[l2].compose
- (*(*node)[faceV[m][0]] * 2 - sum, sum));
-
- (*(*face)[4])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][2]] - *(*node)[faceV[m][0]],
- *(*node)[faceV[m][2]] + *(*node)[faceV[m][0]]) *
-
- *(*node)[faceV[m][1]] *
-
- sclLegendre[l2].compose
- (*(*node)[faceV[m][1]] * 2 - sum, sum));
-
- (*(*face)[5])[i] =
- new Polynomial(intLegendre[l1].compose
- (*(*node)[faceV[m][0]] - *(*node)[faceV[m][2]],
- *(*node)[faceV[m][0]] + *(*node)[faceV[m][2]]) *
+ for(int c = 0; c < 6; c++){
+ unsigned int i = 0;
+
+ for(int l1 = 1; l1 < orderMinus; l1++){
+ for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
+ for(int f = 0; f < 4; f++){
+ Polynomial sum =
+ *(*node)[faceV[c][f][0]] +
+ *(*node)[faceV[c][f][1]] +
+ *(*node)[faceV[c][f][2]];
+
+ (*(*face)[c])[i] =
+ new Polynomial(intLegendre[l1].compose
+ (*(*node)[faceV[c][f][0]] - *(*node)[faceV[c][f][1]],
+ *(*node)[faceV[c][f][0]] + *(*node)[faceV[c][f][1]]) *
- *(*node)[faceV[m][1]] *
+ *(*node)[faceV[c][f][2]] *
- sclLegendre[l2].compose
- (*(*node)[faceV[m][1]] * 2 - sum, sum));
-
- i++;
+ sclLegendre[l2].compose
+ (*(*node)[faceV[c][f][2]] * 2 - sum, sum));
+
+ i++;
+ }
}
}
}
@@ -191,7 +140,7 @@ TetNodeBasis::TetNodeBasis(int order){
Polynomial sub = *(*node)[0] - *(*node)[1];
Polynomial add = *(*node)[0] + *(*node)[1];
- i = 0;
+ unsigned int i = 0;
for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
diff --git a/FunctionSpace/TriEdgeBasis.cpp b/FunctionSpace/TriEdgeBasis.cpp
index 9cac7d26d36bc6a9e3d9a34999000385ede6019d..332bf3fa067f370f01b43f4ad3f3aa14b0828661 100644
--- a/FunctionSpace/TriEdgeBasis.cpp
+++ b/FunctionSpace/TriEdgeBasis.cpp
@@ -15,7 +15,7 @@ TriEdgeBasis::TriEdgeBasis(int order){
nFace = 0;
nCell = ((order - 1) * order + order - 1);
- size = (order + 1) * (order + 2);
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
const int orderPlus = order + 1;
@@ -30,30 +30,16 @@ TriEdgeBasis::TriEdgeBasis(int order){
Polynomial lagrangeSum [3];
Polynomial lagrangeSub [2][3];
- // Classical and Intrated-Scaled Legendre Polynomial //
+ // Legendre Polynomial //
Legendre::legendre(legendre, order);
Legendre::intScaled(intLegendre, orderPlus);
- // 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 i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = lagrange[j] + lagrange[i];
-
- // Lagrange Sub (& revert) //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
- lagrangeSub[0][i] = lagrange[i] - lagrange[j];
- lagrangeSub[1][i] = lagrange[j] - lagrange[i];
- }
-
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][3][2] =
+ {
+ { {0, 1}, {1, 2}, {2, 0} },
+ { {1, 0}, {2, 1}, {0, 2} }
+ };
// Basis //
node = new vector<vector<Polynomial>*>(nVertex);
@@ -63,66 +49,76 @@ TriEdgeBasis::TriEdgeBasis(int order){
(*edge)[0] = new vector<vector<Polynomial>*>(nEdge);
(*edge)[1] = new vector<vector<Polynomial>*>(nEdge);
-
- // Counter
- unsigned int i = 0;
-
- // Edge Based (Nedelec) //
- for(int j = 1; i < 3; j = (j + 1) % 3){
- // Temp
- vector<Polynomial> tmp = lagrange[j].gradient();
-
- // Edge[0]
- tmp[0].mul(lagrange[i]);
- tmp[1].mul(lagrange[i]);
- tmp[2].mul(lagrange[i]);
- (*(*edge)[0])[i] = new vector<Polynomial>(lagrange[i].gradient());
+ // Lagrange //
+ lagrange[0] =
+ Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);
- (*(*edge)[0])[i]->at(0).mul(lagrange[j]);
- (*(*edge)[0])[i]->at(1).mul(lagrange[j]);
- (*(*edge)[0])[i]->at(2).mul(lagrange[j]);
+ lagrange[1] =
+ Polynomial(1, 1, 0, 0);
- (*(*edge)[0])[i]->at(0).sub(tmp[0]);
- (*(*edge)[0])[i]->at(1).sub(tmp[1]);
- (*(*edge)[0])[i]->at(2).sub(tmp[2]);
+ lagrange[2] =
+ Polynomial(1, 0, 1, 0);
- //Edge[1]
- tmp = lagrange[i].gradient();
+ // 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]];
+ }
- tmp[0].mul(lagrange[j]);
- tmp[1].mul(lagrange[j]);
- tmp[2].mul(lagrange[j]);
+
+ // 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]]);
+
+ (*(*edge)[c])[e] =
+ new vector<Polynomial>(lagrange[edgeV[c][e][0]].gradient());
- (*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
+ (*(*edge)[c])[e]->at(0).mul(lagrange[edgeV[c][e][1]]);
+ (*(*edge)[c])[e]->at(1).mul(lagrange[edgeV[c][e][1]]);
+ (*(*edge)[c])[e]->at(2).mul(lagrange[edgeV[c][e][1]]);
+
+ (*(*edge)[c])[e]->at(0).sub(tmp[0]);
+ (*(*edge)[c])[e]->at(1).sub(tmp[1]);
+ (*(*edge)[c])[e]->at(2).sub(tmp[2]);
+ }
+ }
- (*(*edge)[1])[i]->at(0).mul(lagrange[i]);
- (*(*edge)[1])[i]->at(1).mul(lagrange[i]);
- (*(*edge)[1])[i]->at(2).mul(lagrange[i]);
-
- (*(*edge)[1])[i]->at(0).sub(tmp[0]);
- (*(*edge)[1])[i]->at(1).sub(tmp[1]);
- (*(*edge)[1])[i]->at(2).sub(tmp[2]);
- // Next
- i++;
- }
-
// Edge Based (High Order) //
- for(int l = 1; l < orderPlus; l++){
- for(int e = 0; e < 3; e++){
- (*(*edge)[0])[i] =
- new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[0][e],
- lagrangeSum[e])).gradient());
+ for(int c = 0; c < 2; c++){
+ unsigned int i = 0;
- (*(*edge)[1])[i] =
- new vector<Polynomial>((intLegendre[l].compose(lagrangeSub[1][e],
- lagrangeSum[e])).gradient());
- i++;
+ for(int l = 1; l < orderPlus; l++){
+ for(int e = 0; e < 3; e++){
+ (*(*edge)[c])[i + 3] =
+ new vector<Polynomial>
+ ((intLegendre[l].compose(lagrangeSub[c][e],
+ lagrangeSum[e])).gradient());
+
+ i++;
+ }
}
}
-
+
+
// Cell Based (Preliminary) //
Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
@@ -132,7 +128,7 @@ TriEdgeBasis::TriEdgeBasis(int order){
v[l].mul(lagrange[2]);
}
- i = 0;
+ unsigned int i = 0;
// Cell Based (Type 1) //
for(int l1 = 1; l1 < order; l1++){
diff --git a/FunctionSpace/TriNedelecBasis.cpp b/FunctionSpace/TriNedelecBasis.cpp
index c5746c35b0ee722ede79c050a1ab452b3c24fa50..929bdb43b9325736d497f3e6e6686825c9117154 100644
--- a/FunctionSpace/TriNedelecBasis.cpp
+++ b/FunctionSpace/TriNedelecBasis.cpp
@@ -16,6 +16,23 @@ TriNedelecBasis::TriNedelecBasis(void){
size = 3;
+ // Vertices definig Edges & Permutations //
+ const int edgeV[2][3][2] =
+ {
+ { {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);
+
+
// Lagrange //
Polynomial lagrange[3];
@@ -29,53 +46,27 @@ TriNedelecBasis::TriNedelecBasis(void){
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);
-
-
// Nedelec //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
- // Temp
- vector<Polynomial> tmp = lagrange[j].gradient();
-
- // Edge[0]
- tmp[0].mul(lagrange[i]);
- tmp[1].mul(lagrange[i]);
- tmp[2].mul(lagrange[i]);
-
- (*(*edge)[0])[i] = new vector<Polynomial>(lagrange[i].gradient());
-
- (*(*edge)[0])[i]->at(0).mul(lagrange[j]);
- (*(*edge)[0])[i]->at(1).mul(lagrange[j]);
- (*(*edge)[0])[i]->at(2).mul(lagrange[j]);
-
- (*(*edge)[0])[i]->at(0).sub(tmp[0]);
- (*(*edge)[0])[i]->at(1).sub(tmp[1]);
- (*(*edge)[0])[i]->at(2).sub(tmp[2]);
-
- //Edge[1]
- tmp = lagrange[i].gradient();
-
- tmp[0].mul(lagrange[j]);
- tmp[1].mul(lagrange[j]);
- tmp[2].mul(lagrange[j]);
-
- (*(*edge)[1])[i] = new vector<Polynomial>(lagrange[j].gradient());
-
- (*(*edge)[1])[i]->at(0).mul(lagrange[i]);
- (*(*edge)[1])[i]->at(1).mul(lagrange[i]);
- (*(*edge)[1])[i]->at(2).mul(lagrange[i]);
-
- (*(*edge)[1])[i]->at(0).sub(tmp[0]);
- (*(*edge)[1])[i]->at(1).sub(tmp[1]);
- (*(*edge)[1])[i]->at(2).sub(tmp[2]);
- }
+ 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]]);
+
+ (*(*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]);
+ }
+ }
}
TriNedelecBasis::~TriNedelecBasis(void){
diff --git a/FunctionSpace/TriNodeBasis.cpp b/FunctionSpace/TriNodeBasis.cpp
index 4de94cfb7f75f12bf6fc0c7f8ebd6a0f77effb8b..43a11e65808ddf478fb36f37f6b41d39cb362345 100644
--- a/FunctionSpace/TriNodeBasis.cpp
+++ b/FunctionSpace/TriNodeBasis.cpp
@@ -15,21 +15,27 @@ TriNodeBasis::TriNodeBasis(int order){
nFace = 0;
nCell = (order - 1) * (order - 2) / 2;
- size = (order + 1) * (order + 2) / 2;
+ size = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space //
- Polynomial* legendre = new Polynomial[order];
- Polynomial* intLegendre = new Polynomial[order];
+ const int orderMinus = order - 1;
- Polynomial lagrangeSum[3];
- Polynomial lagrangeSub[2][3];
+ Polynomial* legendre = new Polynomial[order];
+ Polynomial* intLegendre = new Polynomial[order];
- // Classical and Intrated-Scaled Legendre Polynomial //
- const int orderMinus = order - 1;
+ 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] =
+ {
+ { {0, 1}, {1, 2}, {2, 0} },
+ { {1, 0}, {2, 1}, {0, 2} }
+ };
// Basis //
node = new vector<Polynomial*>(nVertex);
@@ -55,13 +61,20 @@ TriNodeBasis::TriNodeBasis(int order){
// Lagrange Sum //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3)
- lagrangeSum[i] = *(*node)[i] + *(*node)[j];
+ for(int e = 0; e < 3; e++)
+ lagrangeSum[e] =
+ *(*node)[edgeV[0][e][0]] +
+ *(*node)[edgeV[0][e][1]];
// Lagrange Sub //
- for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
- lagrangeSub[0][i] = *(*node)[j] - *(*node)[i];
- lagrangeSub[1][i] = *(*node)[i] - *(*node)[j];
+ for(int e = 0; e < 3; e++){
+ lagrangeSub[0][e] =
+ *(*node)[edgeV[0][e][0]] -
+ *(*node)[edgeV[0][e][1]];
+
+ lagrangeSub[1][e] =
+ *(*node)[edgeV[1][e][0]] -
+ *(*node)[edgeV[1][e][1]];
}
@@ -96,6 +109,7 @@ TriNodeBasis::TriNodeBasis(int order){
}
}
+
// Free Temporary Sapce //
delete[] legendre;
delete[] intLegendre;