Minor Fix: TriNodeBasis face functions -- Add QuadFaces to ReferenceSpace --...

Minor Fix: TriNodeBasis face functions -- Add QuadFaces to ReferenceSpace -- QuadNodeBasis: Edges OK -- QuadFace: Orientation is *WRONG*

FunctionSpace/BasisGenerator.cpp

@@ -7,14 +7,14 @@
 #include "LineNedelecBasis.h"
 #include "LineLagrangeBasis.h"
 
-#include "QuadNodeBasis.h"
-#include "QuadEdgeBasis.h"
-
 #include "TriNodeBasis.h"
 #include "TriEdgeBasis.h"
 #include "TriNedelecBasis.h"
 #include "TriLagrangeBasis.h"
 
+#include "QuadNodeBasis.h"
+#include "QuadEdgeBasis.h"
+
 #include "TetNodeBasis.h"
 #include "TetEdgeBasis.h"
 #include "TetNedelecBasis.h"
@@ -113,8 +113,11 @@ BasisLocal* BasisGenerator::triHierarchicalGen(unsigned int basisType,
 BasisLocal* BasisGenerator::quaHierarchicalGen(unsigned int basisType,
 					       unsigned int order){
   switch(basisType){
-    //case  0: return new QuadNodeBasis(order);
-    //case  1: return new QuadEdgeBasis(order);
+  case  0: return new QuadNodeBasis(order);
+  case  1:
+    if (order == 0) throw Exception("Nedelec not implemented on Quads");
+    else            throw Exception("1-form not implemented on Quads");
+
   case  2: throw Exception("2-form not implemented on Quads");
   case  3: throw Exception("3-form not implemented on Quads");
 

FunctionSpace/CMakeLists.txt

@@ -26,22 +26,22 @@ set(SRC
   LineNedelecBasis.cpp
   LineLagrangeBasis.cpp
 
-#  QuadNodeBasis.cpp
-#  QuadEdgeBasis.cpp
-
   TriNodeBasis.cpp
   TriEdgeBasis.cpp
   TriNedelecBasis.cpp
   TriLagrangeBasis.cpp
 
-#  HexNodeBasis.cpp
-#  HexEdgeBasis.cpp
+  QuadNodeBasis.cpp
+#  QuadEdgeBasis.cpp
 
   TetNodeBasis.cpp
   TetEdgeBasis.cpp
   TetNedelecBasis.cpp
   TetLagrangeBasis.cpp
 
+#  HexNodeBasis.cpp
+#  HexEdgeBasis.cpp
+
   FunctionSpace.cpp
   FunctionSpaceScalar.cpp
   FunctionSpaceVector.cpp

FunctionSpace/LineReferenceSpace.cpp

@@ -19,8 +19,10 @@ LineReferenceSpace::LineReferenceSpace(void){
   nFace   = 0;
   refFace = NULL;
 
-  // Init All //
-  init();
+  // Init All (Use Tri Face --     //
+  //           unused --           //
+  //           just for interface) //
+  init(0);
 }
 
 LineReferenceSpace::~LineReferenceSpace(void){

FunctionSpace/QuadNodeBasis.cpp

 #include "QuadNodeBasis.h"
+#include "QuadReferenceSpace.h"
 #include "Legendre.h"
 
 using namespace std;
 
-QuadNodeBasis::QuadNodeBasis(int order){
+QuadNodeBasis::QuadNodeBasis(unsigned int order){
+  // Reference Space //
+  refSpace  = new QuadReferenceSpace;
+  nRefSpace = refSpace->getNPermutation();
+
+  const vector<const vector<const vector<unsigned int>*>*>&
+    edgeV = refSpace->getAllEdge();
+
+  const vector<const vector<const vector<unsigned int>*>*>&
+    faceV = refSpace->getAllFace();
+
   // Set Basis Type //
   this->order = order;
 
   type = 0;
   dim  = 2;
 
-  nVertex = 4;
-  nEdge   = 4 * (order - 1);
-  nFace   =     (order - 1) * (order - 1);
-  nCell   = 0;
-
-  nEdgeClosure = 2;
-  nFaceClosure = 0;
-
-  size = nVertex + nEdge + nFace + nCell;
-
-  // Alloc Temporary Space //
-  Polynomial* legendre = new Polynomial[order];
-  Polynomial  lifting[4];
-  Polynomial  liftingSub[2][4];
+  nVertex   = 4;
+  nEdge     = 4 * (order - 1);
+  nFace     =     (order - 1) * (order - 1);
+  nCell     = 0;
+  nFunction = nVertex + nEdge + nFace + nCell;
 
   // Legendre Polynomial //
+  Polynomial* legendre = new Polynomial[order];
   Legendre::integrated(legendre, order);
 
-  // Vertices definig Edges & Permutations //
-  const int edgeV[2][4][2] =
+  // Lagrange & Lifting //
+  const Polynomial lagrange[4] =
     {
-      { {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);
+      Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+                 (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
 
-  (*edge)[0] = new vector<Polynomial*>(nEdge);
-  (*edge)[1] = new vector<Polynomial*>(nEdge);
+      Polynomial((Polynomial(1, 1, 0, 0)) *
+                 (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
 
+      Polynomial((Polynomial(1, 1, 0, 0)) *
+                 (Polynomial(1, 0, 1, 0))),
 
-  // Lifting //
-  lifting[0] =
-    (Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
-    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
-
-  lifting[1] =
-    (Polynomial(1, 1, 0, 0)) +
-    (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0));
-
-  lifting[2] =
-    (Polynomial(1, 1, 0, 0)) +
-    (Polynomial(1, 0, 1, 0));
-
-  lifting[3] =
-    (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]];
+      Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
+                 (Polynomial(1, 0, 1, 0)))
+    };
 
-    liftingSub[1][e] =
-      lifting[edgeV[1][e][0]] -
-      lifting[edgeV[1][e][1]];
-  }
+  const Polynomial lifting[4] =
+    {
+      Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+                 (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
 
+      Polynomial((Polynomial(1, 1, 0, 0)) +
+                 (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0))),
 
-  // Vertex Based (Lagrange) //
-  (*node)[0] =
-    new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
-		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
+      Polynomial((Polynomial(1, 1, 0, 0)) +
+                 (Polynomial(1, 0, 1, 0))),
 
-  (*node)[1] =
-    new Polynomial((Polynomial(1, 1, 0, 0)) *
-		   (Polynomial(1, 0, 0, 0) - Polynomial(1, 0, 1, 0)));
+      Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) +
+                 (Polynomial(1, 0, 1, 0)))
+    };
 
-  (*node)[2] =
-    new Polynomial((Polynomial(1, 1, 0, 0)) *
-		   (Polynomial(1, 0, 1, 0)));
+  // Basis //
+  basis = new Polynomial**[nRefSpace];
 
-  (*node)[3] =
-    new Polynomial((Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0)) *
-		   (Polynomial(1, 0, 1, 0)));
+  for(unsigned int s = 0; s < nRefSpace; s++)
+    basis[s] = new Polynomial*[nFunction];
 
+  // Vertex Based //
+  for(unsigned int s = 0; s < nRefSpace; s++){
+    basis[s][0] = new Polynomial(lagrange[0]);
+    basis[s][1] = new Polynomial(lagrange[1]);
+    basis[s][2] = new Polynomial(lagrange[2]);
+    basis[s][3] = new Polynomial(lagrange[3]);
+  }
 
   // Edge Based //
-  for(int c = 0; c < 2; c++){
-    unsigned int i = 0;
-
-    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]]));
+  for(unsigned int s = 0; s < nRefSpace; s++){
+    unsigned int i = nVertex;
+
+    for(unsigned int e = 0; e < 4; e++){
+      for(unsigned int l = 1; l < order; l++){
+	basis[s][i] =
+	  new Polynomial(legendre[l].compose(lifting[(*(*edgeV[s])[e])[1]] -
+                                             lifting[(*(*edgeV[s])[e])[0]])
+                         *
+			 (lagrange[(*(*edgeV[s])[e])[0]] +
+                          lagrange[(*(*edgeV[s])[e])[1]]));
 
 	i++;
       }
     }
   }
 
+  // Face Based //
 
-  // Cell Based //
-  Polynomial px = Polynomial(2, 1, 0, 0);
-  Polynomial py = Polynomial(2, 0, 1, 0);
+  // NB: We use (*(*faceV[s])[f])[]
+  //     where f = 0, because triangles
+  //     have only ONE face: the face '0'
 
-  px = px - Polynomial(1, 0, 0, 0);
-  py = py - Polynomial(1, 0, 0, 0);
+  for(unsigned int s = 0; s < nRefSpace; s++){
+    unsigned int i = nVertex + nEdge;
 
-  unsigned int i = 0;
+    for(unsigned int l1 = 1; l1 < order; l1++){
+      for(unsigned int l2 = 1; l2 < order; l2++){
+       basis[s][i] =
+         new Polynomial(legendre[l1].compose(lifting[(*(*faceV[s])[0])[0]] -
+                                             lifting[(*(*faceV[s])[0])[1]])
+                         *
 
-  for(int l1 = 1; l1 < order; l1++){
-    for(int l2 = 1; l2 < order; l2++){
-      (*cell)[i] =
-	new Polynomial(legendre[l1].compose(px) * legendre[l2].compose(py));
+                        legendre[l2].compose(lifting[(*(*faceV[s])[0])[0]] -
+                                             lifting[(*(*faceV[s])[0])[3]]));
 
-      i++;
+        i++;
+      }
     }
   }
 
+  cout << (*(*faceV[1])[0])[0] << ", "
+       << (*(*faceV[1])[0])[1] << endl
+
+       << (*(*faceV[1])[0])[0] << ", "
+       << (*(*faceV[1])[0])[3] << endl
+
+       << (*(*faceV[1])[0])[2] << endl;
 
   // Mapping to Gmsh Quad //
   // x = (u + 1) / 2
@@ -140,47 +139,30 @@ QuadNodeBasis::QuadNodeBasis(int order){
   Polynomial  mapY(Polynomial(0.5, 0, 1, 0) +
 		   Polynomial(0.5, 0, 0, 0));
 
-  for(int i = 0; i < nVertex; i++)
-    *(*node)[i] = (*node)[i]->compose(mapX, mapY);
-
-  for(int i = 0; i < nEdgeClosure; i++)
-    for(int j = 0; j < nEdge; j++)
-      *(*(*edge)[i])[j] = (*(*edge)[i])[j]->compose(mapX, mapY);
-
-  for(int i = 0; i < nCell; i++)
-    *(*cell)[i] = (*cell)[i]->compose(mapX, mapY);
-
+  for(unsigned int s = 0; s < nRefSpace; s++){
+    for(unsigned int i = 0; i < nFunction; i++){
+      Polynomial* tmp;
+      tmp = basis[s][i];
+      basis[s][i] = new Polynomial(tmp->compose(mapX, mapY));
+      delete tmp;
+    }
+  }
 
   // Free Temporary Sapce //
   delete[] legendre;
 }
 
 QuadNodeBasis::~QuadNodeBasis(void){
-  // Vertex Based //
-  for(int i = 0; i < nVertex; i++)
-    delete (*node)[i];
+  // ReferenceSpace //
+  delete refSpace;
 
-  delete node;
-
-
-  // Edge Based //
-  for(int c = 0; c < 2; c++){
-    for(int i = 0; i < nEdge; i++)
-      delete (*(*edge)[c])[i];
+  // Basis //
+  for(unsigned int i = 0; i < nRefSpace; i++){
+    for(unsigned int j = 0; j < nFunction; j++)
+      delete basis[i][j];
 
-    delete (*edge)[c];
+    delete[] basis[i];
   }
 
-  delete edge;
-
-
-  // Face Based //
-  delete face;
-
-
-  // Cell Based //
-  for(int i = 0; i < nCell; i++)
-    delete (*cell)[i];
-
-  delete cell;
+  delete[] basis;
 }

FunctionSpace/QuadNodeBasis.h

@@ -24,7 +24,7 @@ class QuadNodeBasis: public BasisHierarchical0From{
   //! @param order The order of the Basis
   //!
   //! Returns a new Node-Basis for Quads of the given order
-  QuadNodeBasis(int order);
+  QuadNodeBasis(unsigned int order);
 
   //! @return Deletes this Basis
   //!

FunctionSpace/QuadReferenceSpace.cpp

@@ -28,8 +28,8 @@ QuadReferenceSpace::QuadReferenceSpace(void){
   refFace[0][2] = 2;
   refFace[0][3] = 3;
 
-  // Init All //
-  init();
+  // Init All (Quad Face) //
+  init(1);
 }
 
 QuadReferenceSpace::~QuadReferenceSpace(void){

FunctionSpace/ReferenceSpace.cpp

@@ -36,7 +36,12 @@ ReferenceSpace::ReferenceSpace(void){
 
   // Defining Ref Edge and Face in //
   // Dervived Class                //
-  // And CALL INIT()               //
+  // And CALL init(faceType)       //
+  //                               //
+  // faceType is:                  //
+  //  * 0 if triangular faces      //
+  //  * 1 if quadragular faces     //
+  //  * 3 if NOT IMPLEMENTED       //
 }
 
 ReferenceSpace::~ReferenceSpace(void){
@@ -65,7 +70,13 @@ ReferenceSpace::~ReferenceSpace(void){
   delete face;
 }
 
-void ReferenceSpace::init(void){
+void ReferenceSpace::init(unsigned int faceType){
+  if(faceType > 3)
+    throw Exception("ReferenceSpace: unknown face type");
+
+  if(faceType == 3)
+    throw Exception("ReferenceSpace: prism not implemented");
+
   // Init Root //
   nPerm      = 0;
   nextLeafId = 0;
@@ -96,7 +107,13 @@ void ReferenceSpace::init(void){
 
   // Get Edges & Faces //
   getEdge();
-  getFace();
+
+  switch(faceType){
+  case 0: getTriFace(); break;
+  case 1: getQuaFace(); break;
+
+  default: throw Exception("ReferenceSpace: impossible error");
+  }
 }
 
 void ReferenceSpace::populate(node* pTreeRoot){
@@ -188,7 +205,7 @@ void ReferenceSpace::getEdge(void){
   }
 }
 
-void ReferenceSpace::getFace(void){
+void ReferenceSpace::getTriFace(void){
   vector<const vector<unsigned int>*>* tmp;
   face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
 
@@ -204,6 +221,23 @@ void ReferenceSpace::getFace(void){
   }
 }
 
+void ReferenceSpace::getQuaFace(void){
+  vector<const vector<unsigned int>*>* tmp;
+  face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
+
+  for(unsigned int p = 0; p < nPerm; p++){
+    tmp = new vector<const vector<unsigned int>*>(nFace);
+
+    for(unsigned int i = 0; i < nFace; i++)
+      (*tmp)[i] = inOrder(p,
+			  refFace[i][0],
+			  refFace[i][1],
+			  refFace[i][2],
+                          refFace[i][3]);
+    (*face)[p] = tmp;
+  }
+}
+
 const vector<unsigned int>* ReferenceSpace::
 inOrder(unsigned int permutation,
 	unsigned int a,
@@ -249,6 +283,30 @@ inOrder(unsigned int permutation,
   return inorder;
 }
 
+const std::vector<unsigned int>* ReferenceSpace::
+inOrder(unsigned int permutation,
+	unsigned int a,
+	unsigned int b,
+	unsigned int c,
+        unsigned int d){
+
+  unsigned int v;
+  unsigned int k = 0;
+  vector<unsigned int>* inorder =
+    new vector<unsigned int>(4);
+
+  for(unsigned int i = 0; i < nVertex; i++){
+    v = perm[permutation][i];
+
+    if(v == a || v == b || v == c || v == d){
+      (*inorder)[k] = v;
+      k++;
+    }
+  }
+
+  return inorder;
+}
+
 unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
   // Const_Cast //
   MElement& element = const_cast<MElement&>(elem);

FunctionSpace/ReferenceSpace.h

@@ -80,7 +80,7 @@ class ReferenceSpace{
  protected:
   ReferenceSpace(void);
 
-  void init(void);
+  void init(unsigned int faceType);
   void populate(node* pTreeRoot);
   void destroy(node* node);
 
@@ -89,7 +89,8 @@ class ReferenceSpace{
   void unconnect(void);                   // Unconnects leafs marked before
 
   void getEdge(void);
-  void getFace(void);
+  void getTriFace(void);
+  void getQuaFace(void);
 
   const std::vector<unsigned int>* inOrder(unsigned int permutation,
 					   unsigned int a,
@@ -100,6 +101,12 @@ class ReferenceSpace{
 					   unsigned int b,
 					   unsigned int c);
 
+  const std::vector<unsigned int>* inOrder(unsigned int permutation,
+					   unsigned int a,
+					   unsigned int b,
+					   unsigned int c,
+                                           unsigned int d);
+
   static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a,
 			    const std::pair<unsigned int, MVertex*>& b);
 

FunctionSpace/TetReferenceSpace.cpp

@@ -29,8 +29,8 @@ TetReferenceSpace::TetReferenceSpace(void){
     refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
   }
 
-  // Init All //
-  init();
+  // Init All (Tri Face) //
+  init(0);
 }
 
 TetReferenceSpace::~TetReferenceSpace(void){

FunctionSpace/TriNodeBasis.cpp

@@ -27,13 +27,12 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
   nCell     = 0;
   nFunction = nVertex + nEdge + nFace + nCell;
 
-  // Alloc Some Space //
+  // Legendre Polynomial //
   const int orderMinus = order - 1;
 
   Polynomial* legendre    = new Polynomial[order];
   Polynomial* intLegendre = new Polynomial[order];
 
-  // Legendre Polynomial //
   Legendre::legendre(legendre, orderMinus);
   Legendre::intScaled(intLegendre, order);
 
@@ -66,7 +65,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
   for(unsigned int s = 0; s < nRefSpace; s++){
     unsigned int i = nVertex;
 
-    for(int e = 0; e < 3; e++){
+    for(unsigned int e = 0; e < 3; e++){
       for(unsigned int l = 1; l < order; l++){
 	basis[s][i] =
 	  new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
@@ -84,10 +83,6 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
   // NB: We use (*(*faceV[s])[f])[]
   //     where f = 0, because triangles
   //     have only ONE face: the face '0'
-
-  // TO CHECK: Are Triangles face matching tets ?
-
-  const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0);
   const int orderMinusTwo = order - 2;
 
   for(unsigned int s = 0; s < nRefSpace; s++){
@@ -103,8 +98,11 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
                                                  lagrange[(*(*faceV[s])[0])[1]])
 			 *
 
-			 legendre[l2].compose(lagrange[(*(*faceV[s])[0])[2]])
+			 legendre[l2].compose((lagrange[(*(*faceV[s])[0])[2]] * 2)
+                                              -
+                                              Polynomial(1, 0, 0, 0))
                          *
+
                          lagrange[(*(*faceV[s])[0])[2]]);
 	i++;
       }

FunctionSpace/TriReferenceSpace.cpp

@@ -27,8 +27,8 @@ TriReferenceSpace::TriReferenceSpace(void){
   refFace[0][1] = 1;
   refFace[0][2] = 2;
 
-  // Init All //
-  init();
+  // Init All (Tri Face) //
+  init(0);
 }
 
 TriReferenceSpace::~TriReferenceSpace(void){