diff --git a/FunctionSpace/HexNodeBasis.cpp b/FunctionSpace/HexNodeBasis.cpp
index 7f742391d5c09aa83a7642615737670c9899f96f..6e197f4dd48538221aa8f04fa60a5f53c6f77580 100644
--- a/FunctionSpace/HexNodeBasis.cpp
+++ b/FunctionSpace/HexNodeBasis.cpp
@@ -119,59 +119,51 @@ HexNodeBasis::HexNodeBasis(size_t order){
basis[s][6] = new Polynomial(lagrange[6]);
basis[s][7] = new Polynomial(lagrange[7]);
}
- /*
+
// Edge Based //
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex;
- for(int l = 1; l < order; l++){
- for(int e = 0; e < 12; e++){
- (*basis)[i] = new Polynomial(
- legendre[l].compose(lifting[edge1[e]] - lifting[edge2[e]]) *
- (*(*basis)[edge1[e]] + *(*basis)[edge2[e]]));
+ for(size_t e = 0; e < 12; e++){
+ for(size_t l = 1; l < order; l++){
+ basis[s][i] =
+ new Polynomial(legendre[l].compose(lifting[edgeIdx[s][e][1]] -
+ lifting[edgeIdx[s][e][0]])
+ *
+ (lagrange[edgeIdx[s][e][0]] +
+ lagrange[edgeIdx[s][e][1]]));
- i++;
+ i++;
+ }
}
}
-
- // Face Based (Preliminary) //
- // Points definig Faces
- const int face1[6] = {0, 3, 2, 1, 5, 4};
- const int face2[6] = {1, 7, 6, 0, 6, 7};
- const int face3[6] = {2, 6, 5, 4, 7, 3};
- const int face4[6] = {3, 2, 1, 5, 4, 0};
-
- // 'Xi' Functions
- for(int f = 0; f < 6; f++)
- xi[f] = lifting[face1[f]] - lifting[face2[f]];
-
- // 'Eta' Functions
- for(int f = 0; f < 6; f++)
- eta[f] = lifting[face1[f]] - lifting[face4[f]];
-
- // 'Lambda' Functions
- for(int f = 0; f < 6; f++)
- lambda[f] =
- *(*basis)[face1[f]] +
- *(*basis)[face2[f]] +
- *(*basis)[face3[f]] +
- *(*basis)[face4[f]];
-
-
// Face Based //
- for(int l1 = 1; l1 < order; l1++){
- for(int l2 = 1; l2 < order; l2++){
- for(int f = 0; f < 6; f++){
- (*basis)[i] = new Polynomial(
- legendre[l1].compose(xi[f]) *
- legendre[l2].compose(eta[f]) *
- lambda[f]);
-
- i++;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex + nEdge;
+
+ for(size_t f = 0; f < 6; f++){
+ for(size_t l1 = 1; l1 < order; l1++){
+ for(size_t l2 = 1; l2 < order; l2++){
+ Polynomial sum =
+ lagrange[faceIdx[s][f][0]] +
+ lagrange[faceIdx[s][f][1]] +
+ lagrange[faceIdx[s][f][2]] +
+ lagrange[faceIdx[s][f][3]];
+
+ basis[s][i] =
+ new Polynomial(legendre[l1].compose(lifting[faceIdx[s][f][0]] -
+ lifting[faceIdx[s][f][1]]) *
+
+ legendre[l2].compose(lifting[faceIdx[s][f][0]] -
+ lifting[faceIdx[s][f][3]]) *
+ sum);
+ i++;
+ }
}
}
}
-
// Cell Based //
Polynomial px = Polynomial(2, 1, 0, 0);
Polynomial py = Polynomial(2, 0, 1, 0);
@@ -181,28 +173,50 @@ HexNodeBasis::HexNodeBasis(size_t order){
py = py - Polynomial(1, 0, 0, 0);
pz = pz - Polynomial(1, 0, 0, 0);
- for(int l1 = 1; l1 < order; l1++){
- for(int l2 = 1; l2 < order; l2++){
- for(int l3 = 1; l3 < order; l3++){
- (*basis)[i] =
- new Polynomial(legendre[l1].compose(px) *
- legendre[l2].compose(py) *
- legendre[l3].compose(pz));
-
- i++;
+ for(size_t s = 0; s < nRefSpace; s++){
+ size_t i = nVertex + nEdge + nFace;
+
+ for(size_t l1 = 1; l1 < order; l1++){
+ for(size_t l2 = 1; l2 < order; l2++){
+ for(size_t l3 = 1; l3 < order; l3++){
+ basis[s][i] =
+ new Polynomial(legendre[l1].compose(px) *
+ legendre[l2].compose(py) *
+ legendre[l3].compose(pz));
+
+ i++;
+ }
}
}
}
+ // Mapping to Gmsh Quad //
+ // x = (u + 1) / 2
+ // y = (v + 1) / 2
+ //
+ // (x, y) = Zaglmayr Ref Quad
+ // (u, v) = Gmsh Ref Quad
+
+ Polynomial mapX(Polynomial(0.5, 1, 0, 0) +
+ Polynomial(0.5, 0, 0, 0));
+
+ Polynomial mapY(Polynomial(0.5, 0, 1, 0) +
+ Polynomial(0.5, 0, 0, 0));
+
+ Polynomial mapZ(Polynomial(0.5, 0, 0, 1) +
+ Polynomial(0.5, 0, 0, 0));
+
+ for(size_t s = 0; s < nRefSpace; s++){
+ for(size_t i = 0; i < nFunction; i++){
+ Polynomial* tmp;
+ tmp = basis[s][i];
+ basis[s][i] = new Polynomial(tmp->compose(mapX, mapY, mapZ));
+ delete tmp;
+ }
+ }
// Free Temporary Sapce //
delete[] legendre;
- delete[] lifting;
-
- delete[] xi;
- delete[] eta;
- delete[] lambda;
-*/
}
HexNodeBasis::~HexNodeBasis(void){
diff --git a/FunctionSpace/Polynomial.cpp b/FunctionSpace/Polynomial.cpp
index e120bf25182f5ce3c78e91dfdd8956b770cbed65..b664eb9bee5902c13e2c5691cb162eb95039c35e 100644
--- a/FunctionSpace/Polynomial.cpp
+++ b/FunctionSpace/Polynomial.cpp
@@ -289,6 +289,17 @@ Polynomial Polynomial::compose(const Polynomial& otherA,
return polynomialFromStack(stk);
}
+Polynomial Polynomial::compose(const Polynomial& otherA,
+ const Polynomial& otherB,
+ const Polynomial& otherC) const{
+ stack<monomial_t> stk;
+
+ for(int i = 0; i < nMon; i++)
+ compose(&mon[i], otherA, otherB, otherC, &stk);
+
+ return polynomialFromStack(stk);
+}
+
void Polynomial::operator=(const Polynomial& other){
if(mon)
delete[] mon;
@@ -603,6 +614,26 @@ void Polynomial::compose(const Polynomial::monomial_t* src,
}
}
+void Polynomial::compose(const Polynomial::monomial_t* src,
+ Polynomial compA, Polynomial compB, Polynomial compC,
+ stack<Polynomial::monomial_t>* stk){
+
+ compA.power(src->power[0]);
+ compB.power(src->power[1]);
+ compC.power(src->power[2]);
+
+ compA.mul(compB);
+ compA.mul(compC);
+ compA.mul(src->coef);
+
+ const int size = compA.nMon;
+
+ for(int i = 0; i < size; i++){
+ if(compA.mon[i].coef != 0)
+ stk->push(compA.mon[i]);
+ }
+}
+
Polynomial Polynomial::
polynomialFromStack(std::stack<Polynomial::monomial_t>& stk){
Polynomial newP;
diff --git a/FunctionSpace/Polynomial.h b/FunctionSpace/Polynomial.h
index b09500d2455e08fedd5a9ff9edd89d8e6e46a988..94951920bfecd6eea054d80a624807d9afbc3be0 100644
--- a/FunctionSpace/Polynomial.h
+++ b/FunctionSpace/Polynomial.h
@@ -68,6 +68,9 @@ class Polynomial{
Polynomial compose(const Polynomial& other) const;
Polynomial compose(const Polynomial& otherA,
const Polynomial& otherB) const;
+ Polynomial compose(const Polynomial& otherA,
+ const Polynomial& otherB,
+ const Polynomial& otherC) const;
void operator=(const Polynomial& other);
@@ -104,6 +107,10 @@ class Polynomial{
Polynomial compA, Polynomial compB,
std::stack<monomial_t>* stk);
+ static void compose(const monomial_t* src,
+ Polynomial compA, Polynomial compB, Polynomial compC,
+ std::stack<monomial_t>* stk);
+
static Polynomial polynomialFromStack(std::stack<monomial_t>& stk);
static monomial_t* copyMonomial(const monomial_t* src, int size);