HexEdgeBasis seems ok

FunctionSpace/HexEdgeBasis.cpp

+#include <vector>
 #include "HexEdgeBasis.h"
 #include "Legendre.h"
 
+using namespace std;
+
 HexEdgeBasis::HexEdgeBasis(const int order){
   // Set Basis Type //
   this->order = order;
@@ -16,16 +19,27 @@ HexEdgeBasis::HexEdgeBasis(const int order){
   Polynomial* legendre      = new Polynomial[orderPlus];
   Polynomial* intLegendre   = new Polynomial[orderPlus];
 
+  Polynomial* xi            = new Polynomial[6];
+  Polynomial* eta           = new Polynomial[6];
+  Polynomial* lambda        = new Polynomial[6];
+
+  vector<Polynomial>* grXi  = new vector<Polynomial>[6];
+  vector<Polynomial>* grEta = new vector<Polynomial>[6];
+
+  Polynomial** iLegendreXi  = new Polynomial*[orderPlus];
+  Polynomial** iLegendreEta = new Polynomial*[orderPlus];
+  Polynomial** legendreXi   = new Polynomial*[orderPlus];
+  Polynomial** legendreEta  = new Polynomial*[orderPlus];
+
   Polynomial* iLegendreX    = new Polynomial[orderPlus];
   Polynomial* iLegendreY    = new Polynomial[orderPlus];
-  Polynomial* legendreX   = new Polynomial[orderPlus];
-  Polynomial* legendreY   = new Polynomial[orderPlus];
+  Polynomial* iLegendreZ    = new Polynomial[orderPlus];
 
   Polynomial* lagrange      = new Polynomial[8];
-  Polynomial* lagrangeSum = new Polynomial[8];
+  Polynomial* lagrangeSum   = new Polynomial[12];
 
   Polynomial* lifting       = new Polynomial[8];
-  Polynomial* liftingSub  = new Polynomial[8];
+  Polynomial* liftingSub    = new Polynomial[12];
 
 
   // Integrated and classical Legendre Polynomial //
@@ -33,6 +47,11 @@ HexEdgeBasis::HexEdgeBasis(const int order){
   Legendre::legendre(legendre, order);
 
 
+  // Points definig Edges //
+  int edge1[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3};
+  int edge2[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7};
+
+
   // Lagrange // 
   lagrange[0] = 
     (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
@@ -75,8 +94,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
      Polynomial(1, 0, 0, 1);
 
   // Lagrange Sum //
-  for(int i = 0, j = 1; i < 8; i++, j = (j + 1) % 8)
-    lagrangeSum[i] = lagrange[i] + lagrange[j];
+  for(int i = 0; i < 12; i++)
+    lagrangeSum[i] = lagrange[edge1[i]] + lagrange[edge2[i]];
 
     
   // Lifting //
@@ -121,8 +140,8 @@ HexEdgeBasis::HexEdgeBasis(const int order){
      Polynomial(1, 0, 0, 1);
 
   // Lifting Sub //
-  for(int i = 0, j = 1; i < 8; i++, j = (j + 1) % 8)
-    liftingSub[i] = lifting[i] - lifting[j];
+  for(int i = 0; i < 12; i++)
+    liftingSub[i] = lifting[edge1[i]] - lifting[edge2[i]];
 
 
   // Basis //
@@ -136,7 +155,7 @@ HexEdgeBasis::HexEdgeBasis(const int order){
   int i = 0; // Function Counter
   Polynomial oneHalf(0.5, 0, 0, 0);
 
-  for(int e = 0; e < 8; e++){
+  for(int e = 0; e < 12; e++){
     (*basis)[i] = 
       (liftingSub[e]).gradient();
     
@@ -154,7 +173,7 @@ HexEdgeBasis::HexEdgeBasis(const int order){
 
   // Edge Based (High Order) //
   for(int l = 1; l < orderPlus; l++){
-    for(int e = 0; e < 8; e++){
+    for(int e = 0; e < 12; e++){
       (*basis)[i] = 
 	(intLegendre[l].compose(liftingSub[e]) * lagrangeSum[e]).gradient();
      
@@ -162,70 +181,259 @@ HexEdgeBasis::HexEdgeBasis(const int order){
     }
   }
 
-  /*  
+
+  // Face Based (Preliminary) //
+  // Points definig Faces
+  int face1[6] = {0, 3, 2, 1, 5, 4};
+  int face2[6] = {1, 7, 6, 0, 6, 7};
+  int face3[6] = {2, 6, 5, 4, 7, 3};
+  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] = 
+      lagrange[face1[f]] +
+      lagrange[face2[f]] +
+      lagrange[face3[f]] +
+      lagrange[face4[f]];
+  
+  // Gradients
+  for(int f = 0; f < 6; f++){
+    grXi[f]  = xi[f].gradient();
+    grEta[f] = eta[f].gradient();
+  }
+
+  // Compositions
+  for(int l = 0; l < orderPlus; l++){
+    iLegendreXi[l]  = new Polynomial[6];
+    iLegendreEta[l] = new Polynomial[6]; 
+ 
+    legendreXi[l]  = new Polynomial[6];   
+    legendreEta[l] = new Polynomial[6];
+
+    for(int f = 0; f < 6; f++){
+      iLegendreXi[l][f]  = intLegendre[l].compose(xi[f]);
+      iLegendreEta[l][f] = intLegendre[l].compose(eta[f]);
+      
+      legendreXi[l][f]   = legendre[l].compose(xi[f]);
+      legendreEta[l][f]  = legendre[l].compose(eta[f]);
+    }
+  }
+ 
+
+  // Face Based (Type 1) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l2 = 1; l2 < orderPlus; l2++){
+      for(int f = 0; f < 6; f++){
+	(*basis)[i] = (iLegendreXi[l1][f]   * 
+		       iLegendreEta[l2][f]  *
+		       lambda[f]).gradient();
+
+	 i++;
+      }
+    }
+  }
+
+
+  // Face Based (Type 2) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l2 = 1; l2 < orderPlus; l2++){
+      for(int f = 0; f < 6; f++){
+	Polynomial tmp1 = 
+	  legendreXi[l1][f]   * 
+	  iLegendreEta[l2][f];
+		       
+	Polynomial tmp2 =
+	  iLegendreXi[l1][f] *
+	  legendreEta[l2][f];
+
+	vector<Polynomial> gr1 = grXi[f];
+	gr1[0].mul(tmp1);
+	gr1[1].mul(tmp1);
+	gr1[2].mul(tmp1);
+
+	vector<Polynomial> gr2 = grEta[f];
+	gr2[0].mul(tmp2);
+	gr2[1].mul(tmp2);
+	gr2[2].mul(tmp2);	
+	
+	(*basis)[i][0] = (gr1[0] - gr2[0]) * lambda[f];
+	(*basis)[i][1] = (gr1[1] - gr2[1]) * lambda[f];
+	(*basis)[i][2] = (gr1[2] - gr2[2]) * lambda[f];
+
+	i++;
+      }
+    }
+  }
+
+
+  // Face Based (Type 3 -- Xi) //
+  for(int l = 1; l < orderPlus; l++){
+    for(int f = 0; f < 6; f++){
+      Polynomial tmp = iLegendreEta[l][f] * lambda[f];
+      
+      (*basis)[i] = grXi[f];
+      
+      (*basis)[i][0].mul(tmp);
+      (*basis)[i][1].mul(tmp);
+      (*basis)[i][2].mul(tmp);
+      
+      i++;   
+    }
+  }
+
+
+  // Face Based (Type 3 -- Eta) //
+  for(int l = 1; l < orderPlus; l++){
+    for(int f = 0; f < 6; f++){
+      Polynomial tmp = iLegendreXi[l][f] * lambda[f];
+      
+      (*basis)[i] = grEta[f];
+      
+      (*basis)[i][0].mul(tmp);
+      (*basis)[i][1].mul(tmp);
+      (*basis)[i][2].mul(tmp);
+      
+      i++;   
+    }
+  }
+
+  
   // Cell Based (Preliminary) //
   Polynomial px   = Polynomial(2, 1, 0, 0);
   Polynomial py   = Polynomial(2, 0, 1, 0);
+  Polynomial pz   = Polynomial(2, 0, 0, 1);
   Polynomial zero = Polynomial(0, 0, 0, 0);
  
   px = px - Polynomial(1, 0, 0, 0);
   py = py - Polynomial(1, 0, 0, 0);
+  pz = pz - Polynomial(1, 0, 0, 0);
 
   for(int l = 0; l < orderPlus; l++){
     iLegendreX[l] = intLegendre[l].compose(px);
     iLegendreY[l] = intLegendre[l].compose(py);
-     legendreX[l] =    legendre[l].compose(px);
-     legendreY[l] =    legendre[l].compose(py);
+    iLegendreZ[l] = intLegendre[l].compose(pz);
   }
 
+
   // Cell Based (Type 1) //
+  int cellStart  = i; // Type one cell base counter
+  int cellNumber = 0;
+
   for(int l1 = 1; l1 < orderPlus; l1++){
     for(int l2 = 1; l2 < orderPlus; l2++){
-      (*basis)[i] = (iLegendreX[l1] * iLegendreY[l2]).gradient();
+      for(int l3 = 1; l3 < orderPlus; l3++){
+	(*basis)[i] = (iLegendreX[l1] * 
+		       iLegendreY[l2] *
+		       iLegendreZ[l3]).gradient();
 
 	i++;
+	cellNumber++;
+      }
     }
   }
   
-  // Cell Based (Type 2) //
-  for(int l1 = 1; l1 < orderPlus; l1++){
+
+  // Cell Based (Type 2 -- First Part) //
+  for(int j = 0; j < cellNumber; j++){
+    int off = j + cellStart;
+ 
+    (*basis)[i][0] = (*basis)[off][0];
+    (*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0); 
+    (*basis)[i][2] = (*basis)[off][2];
+
+    i++;
+  }
+
+
+  // Cell Based (Type 2 -- Second Part) //
+  for(int j = 0; j < cellNumber; j++){
+    int off = j + cellStart;
+
+    (*basis)[i][0] = (*basis)[off][0];
+    (*basis)[i][1] = (*basis)[off][1] * Polynomial(-1, 0, 0, 0); 
+    (*basis)[i][2] = (*basis)[off][2] * Polynomial(-1, 0, 0, 0);
+
+    i++;
+  }
+
+
+  // Cell Based (Type 3 -- First Part) //
   for(int l2 = 1; l2 < orderPlus; l2++){
-      (*basis)[i][0] =  legendreX[l1] * iLegendreY[l2];
-      (*basis)[i][1] = iLegendreX[l1] *  legendreY[l2] * -1;
+    for(int l3 = 1; l3 < orderPlus; l3++){ 
+      (*basis)[i][0] = iLegendreY[l2] * iLegendreZ[l3];
+      (*basis)[i][1] = zero;
       (*basis)[i][2] = zero;
       
       i++;
     }
   }
 
-  // Cell Based (Type 3) //
-  for(int l = 1, iPlus = i + order; l < orderPlus; l++, iPlus++){
-    (*basis)[i][0] = iLegendreY[l];
-    (*basis)[i][1] = zero;
+
+  // Cell Based (Type 3 -- Second Part) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l3 = 1; l3 < orderPlus; l3++){      
+      (*basis)[i][0] = zero;
+      (*basis)[i][1] = iLegendreX[l1] * iLegendreZ[l3];
       (*basis)[i][2] = zero;
       
-    (*basis)[iPlus][0] = zero;
-    (*basis)[iPlus][1] = iLegendreX[l];
-    (*basis)[iPlus][2] = zero;
+      i++;
+    }
+  }  
+
+
+  // Cell Based (Type 3 -- Thrid Part) //
+  for(int l1 = 1; l1 < orderPlus; l1++){
+    for(int l2 = 1; l2 < orderPlus; l2++){
+      (*basis)[i][0] = zero;
+      (*basis)[i][1] = zero;
+      (*basis)[i][2] = iLegendreX[l1] * iLegendreY[l2];
       
       i++;
     }
-  */
+  }  
+
 
   // Free Temporary Sapce //
   delete[] legendre;
   delete[] intLegendre;
 
-  delete[] iLegendreX;
-  delete[] iLegendreY;
-  delete[] legendreX;
-  delete[] legendreY;
-
   delete[] lagrange;
   delete[] lagrangeSum;
 
   delete[] lifting;
   delete[] liftingSub;
+
+  delete[] xi;
+  delete[] eta;
+  delete[] lambda;
+
+  delete[] grXi;
+  delete[] grEta;
+
+  for(int l = 0; l < orderPlus; l++){
+    delete[] iLegendreXi[l];
+    delete[] iLegendreEta[l]; 
+    delete[] legendreXi[l];   
+    delete[] legendreEta[l];
+  }
+
+  delete[] iLegendreXi;
+  delete[] iLegendreEta; 
+  delete[] legendreXi;   
+  delete[] legendreEta;
+
+  delete[] iLegendreX;
+  delete[] iLegendreY;
+  delete[] iLegendreZ;
 }
 
 HexEdgeBasis::~HexEdgeBasis(void){
@@ -236,7 +444,7 @@ HexEdgeBasis::~HexEdgeBasis(void){
 #include <cstdio>
 int main(void){
   const int P = 1;
-  const double d = 0.05;
+  const double d = 0.1;
   const char x[3] = {'X', 'Y', 'Z'};
 
   HexEdgeBasis b(P);
@@ -260,7 +468,7 @@ int main(void){
   printf("p = zeros(%d, 3);\n", b.getSize());
   printf("\n");
 
-  for(int i = 0; i < 16; i++){ //b.getSize(); i++){
+  for(int i = 0; i < b.getSize(); i++){
     for(int j = 0; j < 3; j++)
       printf("p(%d, %d) = %s;\n", i + 1, j + 1, basis[i][j].toString().c_str());
     //printf("p(%d) = %s", i, basis[i].toString().c_str());
@@ -276,7 +484,7 @@ int main(void){
   
   printf("d = %f;\nx = [0:d:1];\ny = x;\nz = x;\n\nlx = length(x);\nly = length(y);\nlz = length(z);\n\n", d);
   
-  for(int i = 0; i < 16; i++) //b.getSize(); i++)
+  for(int i = 0; i < b.getSize(); i++)
     for(int j = 0; j < 3; j++)
       printf("p%d%c = zeros(lx, ly, lz);\n", i + 1, x[j]);
 
@@ -286,7 +494,7 @@ int main(void){
   printf("for k = 1:lz\n");
   printf("\n");
 
-  for(int i = 0; i < 16; i++) //b.getSize(); i++)
+  for(int i = 0; i < b.getSize(); i++)
     printf("[p%dX(j, i, k), p%dY(j, i, k), p%dZ(j, i, k)] = p(%d, x(i), y(j), z(k));\n", i + 1, i + 1, i + 1, i + 1);
 
   printf("\n");
@@ -299,7 +507,7 @@ int main(void){
   printf("\n");
   
   printf("\n");
-  for(int i = 0; i < 16; i++) //b.getSize(); i++)
+  for(int i = b.getSize() - 1; i >= 0; i--)
     printf("figure;\nquiver3(x, y, z, p%dX, p%dY, p%dZ);\n", i + 1, i + 1, i + 1);
   
   printf("\n");

FunctionSpace/HexNodeBasis.cpp

@@ -14,6 +14,10 @@ HexNodeBasis::HexNodeBasis(const int order){
   Polynomial* legendre = new Polynomial[order];
   Polynomial* lifting  = new Polynomial[8];
 
+  Polynomial* xi       = new Polynomial[6];
+  Polynomial* eta      = new Polynomial[6];
+  Polynomial* lambda   = new Polynomial[6];
+
   // Legendre Polynomial //
   Legendre::integrated(legendre, order);
   
@@ -124,26 +128,22 @@ HexNodeBasis::HexNodeBasis(const int order){
   }
 
   
-  // Face Based //
+  // Face Based (Preliminary) //
   // Points definig Faces
   int face1[6] = {0, 3, 2, 1, 5, 4};
   int face2[6] = {1, 7, 6, 0, 6, 7};
   int face3[6] = {2, 6, 5, 4, 7, 3};
   int face4[6] = {3, 2, 1, 5, 4, 0};
 
-  // 'Xi' Function
-  Polynomial*  xi    = new Polynomial[6];
+  // 'Xi' Functions
   for(int f = 0; f < 6; f++)
     xi[f]  = lifting[face1[f]] - lifting[face2[f]];
 
-  // 'Eta' Function
-  Polynomial* eta    = new Polynomial[6];
+  // 'Eta' Functions
   for(int f = 0; f < 6; f++)
     eta[f] = lifting[face1[f]] - lifting[face4[f]];
 
-
-  // 'Lambda' Function
-  Polynomial* lambda = new Polynomial[6];
+  // 'Lambda' Functions
   for(int f = 0; f < 6; f++)
     lambda[f] = 
       (*basis)[face1[f]] +
@@ -152,6 +152,7 @@ HexNodeBasis::HexNodeBasis(const int order){
       (*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++){