** ReferenceSpace:

     -- Minimal number of reference space (at least for tet, tri and quad)
     -- Automatic generation of these reference space
     -- New interface
     -- Modification of TriNodeBasis accordingly

** Need to compute Dof permutations automaticaly

** All Basis are commented (exept for TriNodeBasis)

** WARNING: Memory leak somewhere

FunctionSpace/Basis.h

@@ -81,6 +81,10 @@ class Basis{
   virtual unsigned int getNOrientation(void) const = 0;
   virtual unsigned int getOrientation(const MElement& element) const = 0;
 
+  virtual void mapFromXYZtoABC(const MElement& element,
+                               const fullVector<double>& xyz,
+                               double abc[3]) const = 0;
+
   // Functions Ordering //
   virtual std::vector<size_t>
     getFunctionOrdering(const MElement& element) const = 0;

FunctionSpace/BasisHierarchical0From.cpp

@@ -59,6 +59,12 @@ getOrientation(const MElement& element) const{
   return refSpace->getReferenceSpace(element);
 }
 
+void BasisHierarchical0From::mapFromXYZtoABC(const MElement& element,
+                                             const fullVector<double>& xyz,
+                                             double abc[3]) const{
+  return refSpace->mapFromXYZtoABC(element, xyz, abc);
+}
+
 vector<size_t> BasisHierarchical0From::
 getFunctionOrdering(const MElement& element) const{
   vector<size_t> ordering(nFunction);
@@ -66,6 +72,33 @@ getFunctionOrdering(const MElement& element) const{
   for(size_t i = 0; i < nFunction; i++)
     ordering[i] = i;
 
+  /*
+  if(element.getNum() == 5){
+    ordering[0] = 0;
+    ordering[1] = 1;
+    ordering[2] = 2;
+    ordering[3] = 3;
+    ordering[4] = 4;
+    ordering[5] = 5;
+    ordering[6] = 6;
+    ordering[7] = 7;
+    ordering[8] = 8;
+    ordering[9] = 9;
+  }
+
+  else{
+    ordering[0] = 2;
+    ordering[1] = 0;
+    ordering[2] = 1;
+    ordering[3] = 7;
+    ordering[4] = 8;
+    ordering[5] = 3;
+    ordering[6] = 4;
+    ordering[7] = 5;
+    ordering[8] = 6;
+    ordering[9] = 9;
+  }
+  */
   return ordering;
 }
 

FunctionSpace/BasisHierarchical0From.h

@@ -39,6 +39,10 @@ class BasisHierarchical0From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
+  virtual void mapFromXYZtoABC(const MElement& element,
+                               const fullVector<double>& xyz,
+                               double abc[3]) const;
+
   virtual std::vector<size_t>
     getFunctionOrdering(const MElement& element) const;
 

FunctionSpace/BasisHierarchical1From.cpp

@@ -59,6 +59,12 @@ getOrientation(const MElement& element) const{
   return refSpace->getReferenceSpace(element);
 }
 
+void BasisHierarchical1From::mapFromXYZtoABC(const MElement& element,
+                                             const fullVector<double>& xyz,
+                                             double abc[3]) const{
+  return refSpace->mapFromXYZtoABC(element, xyz, abc);
+}
+
 vector<size_t> BasisHierarchical1From::
 getFunctionOrdering(const MElement& element) const{
   vector<size_t> ordering(nFunction);

FunctionSpace/BasisHierarchical1From.h

@@ -39,6 +39,10 @@ class BasisHierarchical1From: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
+  virtual void mapFromXYZtoABC(const MElement& element,
+                               const fullVector<double>& xyz,
+                               double abc[3]) const;
+
   virtual std::vector<size_t>
     getFunctionOrdering(const MElement& element) const;
 

FunctionSpace/BasisLagrange.cpp

@@ -78,6 +78,11 @@ static size_t matchClosure(vector<int>& reduced,
   return i;
 }
 
+void BasisLagrange::mapFromXYZtoABC(const MElement& element,
+                                    const fullVector<double>& xyz,
+                                    double abc[3]) const{
+}
+
 vector<size_t> BasisLagrange::
 getFunctionOrdering(const MElement& element) const{
   /*

FunctionSpace/BasisLagrange.h

@@ -43,6 +43,10 @@ class BasisLagrange: public BasisLocal{
   virtual unsigned int getNOrientation(void) const;
   virtual unsigned int getOrientation(const MElement& element) const;
 
+  virtual void mapFromXYZtoABC(const MElement& element,
+                               const fullVector<double>& xyz,
+                               double abc[3]) const;
+
   virtual std::vector<size_t>
     getFunctionOrdering(const MElement& element) const;
 

FunctionSpace/FunctionSpaceScalar.cpp

@@ -14,7 +14,7 @@ double FunctionSpaceScalar::
 interpolate(const MElement& element,
             const std::vector<double>& coef,
             const fullVector<double>& xyz) const{
-
+  /*
   // Const Cast For MElement //
   MElement& eelement =
     const_cast<MElement&>(element);
@@ -34,11 +34,11 @@ interpolate(const MElement& element,
   element.getNode(0, abcPnt[0][0], abcPnt[0][1], abcPnt[0][2]);
   element.getNode(1, abcPnt[1][0], abcPnt[1][1], abcPnt[1][2]);
   element.getNode(2, abcPnt[2][0], abcPnt[2][1], abcPnt[2][2]);
-  /*
+
   std::cout << order[0] << " | "
             << order[1] << " | "
             << order[2] << std::endl;
-  */
+
   abcMat[0][0] = abcPnt[0][0];
   abcMat[0][1] = abcPnt[1][0];
   abcMat[0][2] = abcPnt[2][0];
@@ -50,14 +50,14 @@ interpolate(const MElement& element,
   abcMat[2][0] = abcPnt[0][2];
   abcMat[2][1] = abcPnt[1][2];
   abcMat[2][2] = abcPnt[2][2];
-  /*
+
   for(size_t i = 0; i < 3; i++){
     for(size_t j = 0; j < 3; j++)
       std::cout << abcMat[i][j] << "\t";
     std::cout << std::endl;
   }
   std::cout << std::endl;
-  */
+
   double phiUVW[3];
   element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW);
 
@@ -90,12 +90,17 @@ interpolate(const MElement& element,
   uvw[0] = abc[0];
   uvw[1] = abc[1];
   uvw[2] = abc[2];
+  */
+
+  // Get ABC Space coordinate //
+  double abc[3];
+  (*basis)[0]->mapFromXYZtoABC(element, xyz, abc);
 
   // Get Basis Functions //
   const unsigned int nFun = (*basis)[0]->getNFunction();
   fullMatrix<double>  fun(nFun, 1);
 
-  (*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
+  (*basis)[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
 
   // Interpolate (in Reference Place) //
   double val = 0;

FunctionSpace/LineEdgeBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 LineEdgeBasis::LineEdgeBasis(unsigned int order){
+  /*
   // Reference Space //
   refSpace  = new LineReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -71,9 +72,11 @@ LineEdgeBasis::LineEdgeBasis(unsigned int order){
 
   // Free Temporary Space //
   delete[] intLegendre;
+  */
 }
 
 LineEdgeBasis::~LineEdgeBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -86,4 +89,5 @@ LineEdgeBasis::~LineEdgeBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/LineNedelecBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 LineNedelecBasis::LineNedelecBasis(void){
+  /*
   // Reference Space //
   refSpace  = new LineReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -41,9 +42,11 @@ LineNedelecBasis::LineNedelecBasis(void){
   // Nedelec //
   basis[0][0] = new vector<Polynomial>(first);
   basis[1][0] = new vector<Polynomial>(second);
+  */
 }
 
 LineNedelecBasis::~LineNedelecBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -56,4 +59,5 @@ LineNedelecBasis::~LineNedelecBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/LineNodeBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 LineNodeBasis::LineNodeBasis(unsigned int order){
+  /*
   // Reference Space //
   refSpace  = new LineReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -66,9 +67,11 @@ LineNodeBasis::LineNodeBasis(unsigned int order){
 
   // Free Temporary Sapce //
   delete[] intLegendre;
+  */
 }
 
 LineNodeBasis::~LineNodeBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -81,4 +84,5 @@ LineNodeBasis::~LineNodeBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/LineReferenceSpace.cpp

@@ -4,6 +4,7 @@
 using namespace std;
 
 LineReferenceSpace::LineReferenceSpace(void){
+  /*
   // Vertex Definition //
   nVertex = 2;
 
@@ -21,20 +22,23 @@ LineReferenceSpace::LineReferenceSpace(void){
 
   // Init All //
   init();
+  */
 }
 
 LineReferenceSpace::~LineReferenceSpace(void){
+  /*
   // Delete Ref Edge //
   for(size_t i = 0; i < nEdge; i++)
     delete[] refEdge[i];
 
   delete[] refEdge;
+  */
 }
 
 string LineReferenceSpace::toLatex(void) const{
-  const size_t nPerm = pTree->getNPermutation();
+  //const size_t nPerm = pTree->getNPermutation();
   stringstream stream;
-
+  /*
   stream << "\\documentclass{article}" << endl << endl
 
          << "\\usepackage{longtable}"  << endl
@@ -64,6 +68,6 @@ string LineReferenceSpace::toLatex(void) const{
 
   stream << "\\end{longtable}" << endl
          << "\\end{document}"  << endl;
-
+  */
   return stream.str();
 }

FunctionSpace/QuadEdgeBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
+  /*
   // Reference Space //
   refSpace  = new QuadReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -238,9 +239,11 @@ QuadEdgeBasis::QuadEdgeBasis(unsigned int order){
   // Free Temporary Sapce //
   delete[] legendre;
   delete[] intLegendre;
+  */
 }
 
 QuadEdgeBasis::~QuadEdgeBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -253,4 +256,5 @@ QuadEdgeBasis::~QuadEdgeBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/QuadNedelecBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 QuadNedelecBasis::QuadNedelecBasis(void){
+  /*
   // Reference Space //
   refSpace  = new QuadReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -105,9 +106,11 @@ QuadNedelecBasis::QuadNedelecBasis(void){
       delete old;
     }
   }
+  */
 }
 
 QuadNedelecBasis::~QuadNedelecBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -120,4 +123,5 @@ QuadNedelecBasis::~QuadNedelecBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/QuadNodeBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 QuadNodeBasis::QuadNodeBasis(unsigned int order){
+  /*
   // Reference Space //
   refSpace  = new QuadReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -142,9 +143,11 @@ QuadNodeBasis::QuadNodeBasis(unsigned int order){
 
   // Free Temporary Sapce //
   delete[] legendre;
+  */
 }
 
 QuadNodeBasis::~QuadNodeBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -157,4 +160,5 @@ QuadNodeBasis::~QuadNodeBasis(void){
   }
 
   delete[] basis;
+  */
 }

FunctionSpace/QuadReferenceSpace.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 QuadReferenceSpace::QuadReferenceSpace(void){
+  /*
   // Vertex Definition //
   nVertex = 4;
 
@@ -37,9 +38,11 @@ QuadReferenceSpace::QuadReferenceSpace(void){
 
   // Init All //
   init();
+  */
 }
 
 QuadReferenceSpace::~QuadReferenceSpace(void){
+  /*
   // Delete Ref Edge //
   for(size_t i = 0; i < nEdge; i++)
     delete[] refEdge[i];
@@ -52,13 +55,14 @@ QuadReferenceSpace::~QuadReferenceSpace(void){
 
   delete[] refFace;
   delete[] nNodeInFace;
+  */
 }
 
 string QuadReferenceSpace::toLatex(void) const{
-  const size_t   nPerm = pTree->getNPermutation();
+  //const size_t   nPerm = pTree->getNPermutation();
   stringstream   stream;
-  vector<size_t> perm(nVertex);
-
+  //vector<size_t> perm(nVertex);
+  /*
   stream << "\\documentclass{article}" << endl << endl
 
          << "\\usepackage{longtable}"  << endl
@@ -98,6 +102,6 @@ string QuadReferenceSpace::toLatex(void) const{
 
   stream << "\\end{longtable}" << endl
          << "\\end{document}"  << endl;
-
+  */
   return stream.str();
 }

FunctionSpace/ReferenceSpace.cpp

@@ -7,55 +7,12 @@
 using namespace std;
 
 ReferenceSpace::ReferenceSpace(void){
-  // Init to NULL                  //
-  nEdge           = 0;
-  refEdge         = NULL;
-  //permutedRefEdge = NULL;
-  edge            = NULL;
-
-  nFace           = 0;
-  nNodeInFace     = NULL;
-  refFace         = NULL;
-  //permutedRefFace = NULL;
-  face            = NULL;
-
   // Defining Ref Edge and Face in //
   // Derived Class                 //
 }
 
 ReferenceSpace::~ReferenceSpace(void){
-  const size_t nPerm = pTree->getNPermutation();
-
-  // Destroy Tree //
   delete pTree;
-
-  // Delete Permutated Edge //
-  for(size_t p = 0; p < nPerm; p++){
-    for(size_t i = 0; i < nEdge; i++){
-      //delete[] permutedRefEdge[p][i];
-      delete   (*(*edge)[p])[i];
-    }
-
-    //delete[] permutedRefEdge[p];
-    delete   (*edge)[p];
-  }
-
-  //delete[] permutedRefEdge;
-  delete   edge;
-
-  // Delete Permutated Face //
-  for(size_t p = 0; p < nPerm; p++){
-    for(size_t i = 0; i < nFace; i++){
-      //delete[] permutedRefFace[p][i];
-      delete (*(*face)[p])[i];
-    }
-
-    //delete[] permutedRefFace[p];
-    delete (*face)[p];
-  }
-
-  //delete[] permutedRefFace;
-  delete face;
 }
 
 void ReferenceSpace::init(void){
@@ -67,22 +24,51 @@ void ReferenceSpace::init(void){
 
   pTree = new PermutationTree(vertexSeq);
 
-  // Edges & Faces //
-  getEdge();
-  getFace();
-
   // Cyclic Permutation //
-  findCyclicPermutation();
+  list<size_t>          listOfTrueReferenceSpace;
+  list<vector<size_t> > listOfRefNodeIndexPermutation;
+  list<vector<size_t> > listOfReverseNodeIndexPermutation;
+
+  findCyclicPermutation(listOfTrueReferenceSpace,
+                        listOfRefNodeIndexPermutation,
+                        listOfReverseNodeIndexPermutation);
+
+  // Iterators //
+  list<size_t>::iterator refSpaceIt  = listOfTrueReferenceSpace.begin();
+
+  // Reference Spaces Node Id //
+  const size_t nRefSpace = listOfTrueReferenceSpace.size();
+  refSpaceNodeId.resize(nRefSpace);
+
+  for(size_t i = 0; i < nRefSpace; i++, refSpaceIt++){
+    refSpaceNodeId[i].resize(nVertex);
+
+    pTree->fillWithPermutation(*refSpaceIt, refSpaceNodeId[i]);
+  }
+
+  // ABC space to UVW space shape function index mapping //
+  const size_t nPermutation = pTree->getNPermutation();
+  ABCtoUVWIndex.resize(nPermutation);
+  ABCtoUVWIndex.assign(listOfReverseNodeIndexPermutation.begin(),
+                       listOfReverseNodeIndexPermutation.end());
+
+  // UVW space to ABC space shape function index mapping //
+  UVWtoABCIndex.resize(nPermutation);
+  UVWtoABCIndex.assign(listOfRefNodeIndexPermutation.begin(),
+                       listOfRefNodeIndexPermutation.end());
+
+  // Ordered Edges & Faces Node Index//
+  getOrderedEdge();
+  getOrderedFace();
 }
 
-void ReferenceSpace::findCyclicPermutation(void){
+void ReferenceSpace::
+findCyclicPermutation(list<size_t>&          listOfTrueReferenceSpace,
+                      list<vector<size_t> >& listOfRefNodeIndexPermutation,
+                      list<vector<size_t> >& listOfReverseNodeIndexPermutation){
   // Alloc Some Data //
   const size_t nPerm = pTree->getNPermutation();
 
-  list<size_t>          listOfTrueReferenceSpace;
-  list<vector<size_t> > listOfRefNodeIndexPermutation;
-  list<vector<size_t> > listOfReverseNodeIndexPermutation;
-
   vector<size_t> pTest(nVertex);
   vector<size_t> pRef(nVertex);
 
@@ -90,12 +76,13 @@ void ReferenceSpace::findCyclicPermutation(void){
   list<size_t>::iterator end;
   triple match;
 
+  // Values For unPermutedIndex //
   vector<size_t> unPermutedIndex(nVertex);
 
   for(size_t i = 0; i < nVertex; i++)
     unPermutedIndex[i] = i;
 
-  // Find Cyclic Permutation
+  // Find Cyclic Permutation //
   for(size_t i = 0; i < nPerm; i++){
     // No match
     match.first = false;
@@ -193,19 +180,22 @@ static bool haveSameNode(vector<size_t>& face0,
 size_t ReferenceSpace::findCorrespondingFace(vector<size_t>& face,
                                              vector<size_t>& node){
   // Init Stuff //
+  const size_t nFace    = refFaceNodeIdx.size();
   const size_t faceSize = face.size();
   bool match = false;
   size_t f = 0;
 
   vector<size_t> testFace(faceSize);
+  size_t         nNodeInFace;
 
   // Test All Face until match
   while(!match && f < nFace){
+    nNodeInFace = refFaceNodeIdx[f].size();
 
-    if(nNodeInFace[f] == faceSize){
+    if(nNodeInFace == faceSize){
       // Get face f nodes
       for(size_t i = 0; i < faceSize; i++)
-        testFace[i] = node[refFace[f][i]];
+        testFace[i] = node[refFaceNodeIdx[f][i]];
 
       // Look if match
       match = haveSameNode(testFace, face);
@@ -259,19 +249,21 @@ isCyclicPermutation(vector<size_t>& pTest,
                     vector<size_t>& pRef){
 
   // Node IDs of Reference Space first Face
-  vector<size_t> refNode(nNodeInFace[0]);
+  const size_t   nNodeInFaceZero = refFaceNodeIdx[0].size();
+  vector<size_t> refNode(nNodeInFaceZero);
 
-  for(size_t i = 0; i < nNodeInFace[0]; i++)
-    refNode[i] = pRef[refFace[0][i]];
+  for(size_t i = 0; i < nNodeInFaceZero; i++)
+    refNode[i] = pRef[refFaceNodeIdx[0][i]];
 
   // Corresponding Face in Test Permutation
   size_t testFaceId = findCorrespondingFace(refNode, pTest);
 
   // Node IDs of Test Permutation correspnding Face
-  vector<size_t> testNode(nNodeInFace[testFaceId]);
+  const size_t   nNodeInFaceTest = refFaceNodeIdx[testFaceId].size();
+  vector<size_t> testNode(nNodeInFaceTest);
 
-  for(size_t i = 0; i < nNodeInFace[testFaceId]; i++)
-    testNode[i] = pTest[refFace[testFaceId][i]];
+  for(size_t i = 0; i < nNodeInFaceTest; i++)
+    testNode[i] = pTest[refFaceNodeIdx[testFaceId][i]];
 
   // Return Triple //
   triple tri = {
@@ -283,229 +275,125 @@ isCyclicPermutation(vector<size_t>& pTest,
   return tri;
 }
 
-void ReferenceSpace::getEdge(void){
-  const size_t nPerm = pTree->getNPermutation();
+void ReferenceSpace::getOrderedEdge(void){
+  // Fill orderedEdgeNodeIdx[s][e]                         //
+  // (for all reference space 's' and edge 'e') such that: //
+  //                                                       //
+  // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][0]] <      //
+  // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][1]]        //
 
   // Alloc
-  vector<const vector<unsigned int>*>* tmp;
-  edge = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
-  /*
-  // All Edges have two nodes
-  unsigned int* nNodeInEdge = new unsigned int[nEdge];
-
-  for(unsigned int e = 0; e < nEdge; e++)
-    nNodeInEdge[e] = 2;
-
-  // Get Edge nodes depending on the orientation
-  //    (edge 1 = [1 2] for orientation 1)
-  //    (       = [4 1] for orientation 2)
-  getPermutedRefEntity(&permutedRefEdge,
-                       refEdge,
-                       nNodeInEdge,
-                       nEdge);
-
-  delete[] nNodeInEdge;
-  */
+  const size_t nEdge     = refEdgeNodeIdx.size();
+  const size_t nRefSpace = refSpaceNodeId.size();
+
+  orderedEdgeNodeIdx.resize(nRefSpace);
+
   // Populate Edge
-  for(size_t p = 0; p < nPerm; p++){
-    tmp = new vector<const vector<unsigned int>*>(nEdge);
+  for(size_t s = 0; s < nRefSpace; s++){
+    orderedEdgeNodeIdx[s].resize(nEdge);
 
-    for(size_t e = 0; e < nEdge; e++)
-      (*tmp)[e] = getOrderedEdge(p, e);
+    for(size_t e = 0; e < nEdge; e++){
+      orderedEdgeNodeIdx[s][e].resize(2);
 
-    (*edge)[p] = tmp;
+      orderRefEntityForGivenRefSpace(refEdgeNodeIdx[e],
+                                     refSpaceNodeId[s],
+                                     orderedEdgeNodeIdx[s][e]);
+    }
   }
 }
 
-void ReferenceSpace::getFace(void){
-  const size_t nPerm = pTree->getNPermutation();
+void ReferenceSpace::getOrderedFace(void){
+  // Fill orderedEdgeNodeIdx[s][f]                                //
+  // (for all reference space 's' and face 'f') such that:        //
+  //                                                              //
+  // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][0]] <             //
+  // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][1]] <             //
+  // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][2]] <             //
+  // (refSpaceNodeId[s][orderedFaceNodeIdx[s][f][3]]) (Quad Face) //
+  //                                                              //
+  // If we have a Quad face, a correction is needed such that:    //
+  //                                                              //
+  // orderedFaceNodeIdx[s][f][2]                                  //
+  // is *opposite* to                                             //
+  // orderedFaceNodeIdx[s][f][0]                                  //
 
   // Alloc
-  vector<const vector<unsigned int>*>* tmp;
-  face = new vector<const vector<const vector<unsigned int>*>*>(nPerm);
-  /*
-  // Get Face nodes depending on the orientation
-  //    (face 1 = [1 2 3 4] for orientation 1)
-  //    (       = [4 1 2 3] for orientation 2)
-  getPermutedRefEntity(&permutedRefFace,
-                       refFace,
-                       nNodeInFace,
-                       nFace);
-  */
-  // Populate Face
-  for(size_t p = 0; p < nPerm; p++){
-    tmp = new vector<const vector<unsigned int>*>(nFace);
+  const size_t nFace     = refFaceNodeIdx.size();
+  const size_t nRefSpace = refSpaceNodeId.size();
 
-    for(size_t f = 0; f < nFace; f++){
-      // Dending on the number of node per face
-      // The ordering strategy is different
+  orderedFaceNodeIdx.resize(nRefSpace);
 
-      switch(nNodeInFace[f]){
-      case 3: (*tmp)[f] = getOrderedTriFace(p, f);  break;
-      case 4: (*tmp)[f] = getOrderedQuadFace(p, f); break;
+  // Populate Face
+  for(size_t s = 0; s < nRefSpace; s++){
+    orderedFaceNodeIdx[s].resize(nFace);
 
-      default:
-        throw Exception("I can handle face with %d nodes",
-                        nNodeInFace[f]);
-      }
-    }
+    for(size_t f = 0; f < nFace; f++){
+      size_t nNodeInFace = refFaceNodeIdx[f].size();
 
-    (*face)[p] = tmp;
-  }
-}
-/*
-void ReferenceSpace::getPermutedRefEntity(unsigned int**** permutedRefEntity,
-                                          unsigned int**   refEntity,
-                                          unsigned int*    nNodeInEntity,
-                                          unsigned int     nEntity){
-  // Alloc PermutedRefEntity
-  (*permutedRefEntity) = new unsigned int**[nPerm];
-
-  for(unsigned int p = 0; p < nPerm; p++){
-    (*permutedRefEntity)[p] = new unsigned int*[nEntity];
-
-    for(unsigned int e = 0; e < nEntity; e++)
-      (*permutedRefEntity)[p][e] = new unsigned int[nNodeInEntity[e]];
-  }
+      orderedFaceNodeIdx[s][f].resize(nNodeInFace);
 
-  // Generate Permuted Ref Entity
-  for(unsigned int p = 0; p < nPerm; p++){
-    for(unsigned int e = 0; e < nEntity; e++){
-      for(unsigned int n = 0; n < nNodeInEntity[e]; n++)
-        (*permutedRefEntity)[p][e][n] = perm[p][refEntity[e][n]];
+      orderRefEntityForGivenRefSpace(refFaceNodeIdx[f],
+                                     refSpaceNodeId[s],
+                                     orderedFaceNodeIdx[s][f]);
+      if(nNodeInFace == 4)
+        correctQuadFaceNodeIdx(orderedFaceNodeIdx[s][f]);
     }
   }
 }
-*/
-const vector<unsigned int>* ReferenceSpace::
-getOrderedEdge(unsigned int permutation,
-               unsigned int edge){
-
-  /**************************************************************************
-   * I want to know how I need to take the values of refEdge[edge]          *
-   * (i.e the node *index* of a given edge) so that the corresponding       *
-   * values in permutation (i.e perm[permutation][refEdge[edge]]) are going *
-   * from the smallest value to the biggest.                                *
-   **************************************************************************/
-
-  // Alloc
-  vector<unsigned int>* ordered = new vector<unsigned int>(2);
-
-  // Permutation
-  vector<size_t> perm(nVertex);
-  pTree->fillWithPermutation(permutation, perm);
-
-  // Order refEdge
-  orderRefEntityForGivenPermutation(refEdge[edge],
-                                    perm,
-                                    *ordered);
-
-  // Return ordered
-  return ordered;
-}
 
-const std::vector<unsigned int>* ReferenceSpace::
-getOrderedTriFace(unsigned int permutation,
-                  unsigned int face){
-
-  /***********************************************
-   * Same as getOrderedEdge, but I with refFace. *
-   ***********************************************/
+void ReferenceSpace::
+orderRefEntityForGivenRefSpace(vector<size_t>& refEntityNodeIdx,
+                               vector<size_t>& refSpaceNodeId,
+                               vector<size_t>& orderedEntityNodeIdx){
+  // Get Size
+  const size_t size = orderedEntityNodeIdx.size();
 
-  // Alloc
-  vector<unsigned int>* ordered = new vector<unsigned int>(3);
+  // Ref Entity Node *ID*
+  vector<pair<size_t, size_t> > refEntityNodeId(size);
 
-  // Permutation
-  vector<size_t> perm(nVertex);
-  pTree->fillWithPermutation(permutation, perm);
+  for(size_t i = 0; i < size; i++){
+    refEntityNodeId[i].first  = i;                             // Node Id INDEX
+    refEntityNodeId[i].second = refSpaceNodeId[refEntityNodeIdx[i]]; // Node Id
+  }
 
-  // Order refFace
-  orderRefEntityForGivenPermutation(refFace[face],
-                                    perm,
-                                    *ordered);
+  // Sort it with repsect to second entry (Node ID)
+  std::sort(refEntityNodeId.begin(), refEntityNodeId.end(), sortPredicate);
 
-  // Return ordered
-  return ordered;
+  // Populate orderedEntityNodeIdx (Usign sorted Node ID old Index)
+  for(size_t i = 0; i < size; i++)
+    orderedEntityNodeIdx[i] = refEntityNodeIdx[refEntityNodeId[i].first];
 }
 
-const std::vector<unsigned int>* ReferenceSpace::
-getOrderedQuadFace(unsigned int permutation,
-                   unsigned int face){
-
-  /*********************************************************************
-   * Same as getOrderedFace, but refFace has 4 nodes.                  *
-   *                                                                   *
-   * In that case,                                                     *
-   * I need node at ordered[2] to be *opposite* to node at ordered[0]. *
-   *  --> If not make a permutation such                               *
-   *      that node opposite at ordered[2]                             *
-   *      is at ordered[0]                                             *
-   *********************************************************************/
-
-  // Alloc
-  vector<unsigned int>* ordered = new vector<unsigned int>(4);
-
-  // Permutation
-  vector<size_t> perm(nVertex);
-  pTree->fillWithPermutation(permutation, perm);
+void ReferenceSpace::
+correctQuadFaceNodeIdx(vector<size_t>& correctedQuadFaceNodeIdx){
 
-  // Order refFace
-  orderRefEntityForGivenPermutation(refFace[face],
-                                    perm,
-                                    *ordered);
+  // Get :                       //
+  // correctedQuadFaceNodeIdx[2] //
+  // is *opposite* to            //
+  // correctedQuadFaceNodeIdx[0] //
 
-  // Get ordered[2] opposite to ordered[0]
-  unsigned int opposite = ((*ordered)[0] + 2) % 4;
+  size_t opposite = (correctedQuadFaceNodeIdx[0] + 2) % 4;
 
-  if((*ordered)[2] != opposite){
+  if(correctedQuadFaceNodeIdx[2] != opposite){
     // Find opposite node index
-    unsigned int tmp;
-    unsigned int idx = 1;
-    while((*ordered)[idx] != opposite)
-      idx++;
-
-    // Swap ordered[2] and ordered[idx]
-    tmp = (*ordered)[2];
-    (*ordered)[2]   = opposite;
-    (*ordered)[idx] = tmp;
-  }
-
-  // Return
-  return ordered;
-}
-
-void ReferenceSpace::
-orderRefEntityForGivenPermutation(size_t* refEntity,
-                                  vector<size_t>& permutation,
-                                  vector<unsigned int>& orderedEntity){
-  // Get Size
-  const unsigned int size = orderedEntity.size();
+    size_t tmp;
+    size_t idx = 1;
 
-  // Copy srcSort in sorted
-  vector<pair<unsigned int, unsigned int> > sorted(size);
+    while(correctedQuadFaceNodeIdx[idx] != opposite)
+      idx++;
 
-  for(unsigned int i = 0; i < size; i++){
-    sorted[i].first  = i;
-    sorted[i].second = permutation[refEntity[i]];
+    // Swap correctedQuadFaceNodeIdx[2] and correctedQuadFaceNodeIdx[idx]
+    tmp = correctedQuadFaceNodeIdx[2];
+    correctedQuadFaceNodeIdx[2]   = opposite;
+    correctedQuadFaceNodeIdx[idx] = tmp;
   }
-
-  // Sort it with repsect to second entry
-  std::sort(sorted.begin(), sorted.end(), sortPredicate);
-
-  // Swap with respect to srcSwap
-  for(unsigned int i = 0; i < size; i++)
-    orderedEntity[i] = refEntity[sorted[i].first];
 }
 
-size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
-  // Const_Cast //
-  MElement& element = const_cast<MElement&>(elem);
-
+size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
   // Get Primary Vertices Global ID //
-  const int nVertex = element.getNumPrimaryVertices();
-  vector<pair<unsigned int, unsigned int> > vertexGlobalId(nVertex);
+  vector<pair<size_t, size_t> > vertexGlobalId(nVertex);
 
-  for(int i = 0; i < nVertex; i++){
+  for(size_t i = 0; i < nVertex; i++){
     vertexGlobalId[i].first  = i;
     vertexGlobalId[i].second = element.getVertex(i)->getNum();
   }
@@ -517,7 +405,7 @@ size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
   // Reduce Vertex Global ID //
   vector<size_t> vertexReducedId(nVertex);
 
-  for(int i = 0; i < nVertex; i++)
+  for(size_t i = 0; i < nVertex; i++)
     vertexReducedId[vertexGlobalId[i].first] = i;
 
   // Tree Lookup //
@@ -531,37 +419,114 @@ size_t ReferenceSpace::getReferenceSpace(const MElement& elem) const{
   }
 }
 
+void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
+                                     const fullVector<double>& xyz,
+                                     double abc[3]){
+  // Get UVW coordinate //
+  double phys[3] = {xyz(0), xyz(1), xyz(2)};
+  double uvw[3];
+
+  element.xyz2uvw(phys, uvw);
+
+  // Compute coordinate in ABC Space //
+  // ABC node coordinate
+  double** abcMat = new double*[nVertex];
+
+  for(size_t i = 0; i < nVertex; i++)
+    abcMat[i] = new double[3];
+
+  for(size_t i = 0; i < nVertex; i++)
+    element.getNode(i, abcMat[i][0], abcMat[i][1], abcMat[i][2]);
+
+  // UVW (order 1) shape functions
+  double* phiUVW = new double[nVertex];
+  element.getShapeFunctions(uvw[0], uvw[1], uvw[2], phiUVW, 1);
+
+  // Element Permutation Index //
+  const size_t permutationIdx = getPermutationIdx(element);
+
+  // Map From UVW to ABC //
+  abc[0] = 0;
+  for(size_t i = 0; i < nVertex; i++)
+    abc[0] += abcMat[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+  abc[1] = 0;
+  for(size_t i = 0; i < nVertex; i++)
+    abc[1] += abcMat[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+  abc[2] = 0;
+  for(size_t i = 0; i < nVertex; i++)
+    abc[2] += abcMat[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
+
+  // Free //
+  delete[] phiUVW;
+
+  for(size_t i = 0; i < nVertex; i++)
+    delete[] abcMat[i];
+
+  delete[] abcMat;
+}
+
+void ReferenceSpace::getJacobian(const MElement& element,
+                                 const fullVector<double>& xyz,
+                                 fullMatrix<double>& jac) const{
+}
+
 string ReferenceSpace::toString(void) const{
-  const size_t nPerm = pTree->getNPermutation();
+  const size_t nRefSpace = refSpaceNodeId.size();
+  const size_t nEdge     = refEdgeNodeIdx.size();
+  const size_t nFace     = refFaceNodeIdx.size();
   stringstream stream;
 
-  // ReferenceSpaces //
-  stream << "Reference Spaces:" << endl
+  // Permutation Tree //
+  stream << "Permutation Tree:" << endl
          << pTree->toString()   << endl;
 
-  stream << "Edges Permutations:" << endl;
+  // Reference Space //
+  stream << "Reference Space Node IDs:" << endl;
 
-  for(size_t i = 0; i < nPerm; i++){
-    stream << "  * RefSpace #" << i + 1 << ":" << endl;
+  for(size_t s = 0; s < nRefSpace; s++){
+    stream << "  ** Reference Space #" << s + 1 << ": [";
 
-    for(size_t j = 0; j < nEdge; j++)
-      stream << "      -- ["
-             << edge->at(i)->at(j)->at(0) << ", "
-             << edge->at(i)->at(j)->at(1) << "]" << endl;
+    for(size_t v = 0; v < nVertex - 1; v++)
+      stream << refSpaceNodeId[s][v] << ", ";
+
+    stream << refSpaceNodeId[s][nVertex - 1] << "]" << endl;
   }
 
-  stream << "Faces Permutations:" << endl;
+  stream << endl << "Ordered Edge Node Index:" << endl;
 
-  for(size_t i = 0; i < nPerm; i++){
-    stream << "  * RefSpace #" << i + 1 << ":" << endl;
+  for(size_t s = 0; s < nRefSpace; s++){
+    stream << "  ** Reference Space #" << s + 1 << ":"
+           << endl;
 
-    for(size_t j = 0; j < nFace; j++)
+    for(size_t e = 0; e < nEdge; e++)
       stream << "      -- ["
-             << face->at(i)->at(j)->at(0) << ", "
-             << face->at(i)->at(j)->at(1) << ", "
-             << face->at(i)->at(j)->at(2) << "]" << endl;
+             << orderedEdgeNodeIdx[s][e][0] << ", "
+             << orderedEdgeNodeIdx[s][e][1] << "]"
+             << endl;
+  }
+
+  stream << endl << "Ordered Face Node Index:" << endl;
+
+  for(size_t s = 0; s < nRefSpace; s++){
+    stream << "  ** Reference Space #" << s + 1 << ":"
+           << endl;
+
+    for(size_t f = 0; f < nFace; f++){
+      const size_t nNodeInFace = orderedFaceNodeIdx[s][f].size();
+      stream << "      -- [";
+
+      for(size_t n = 0; n < nNodeInFace - 1; n++)
+        stream << orderedFaceNodeIdx[s][f][n] << ", ";
+
+      stream << orderedFaceNodeIdx[s][f][nNodeInFace - 1] << "]"
+             << endl;
+    }
   }
 
+  stream << endl;
+
   return stream.str();
 }
 

FunctionSpace/ReferenceSpace.h

@@ -3,10 +3,10 @@
 
 #include <vector>
 #include <list>
-#include <stack>
 #include <string>
 
 #include "MElement.h"
+#include "fullMatrix.h"
 #include "PermutationTree.h"
 
 /**
@@ -28,70 +28,51 @@ class ReferenceSpace{
   } triple;
 
  protected:
+  // Element Definition //
+  size_t nVertex;
+  std::vector<std::vector<size_t> > refEdgeNodeIdx;
+  std::vector<std::vector<size_t> > refFaceNodeIdx;
+
   // Permutation Tree //
   PermutationTree* pTree;
 
-  //////////////////////////////////////////////////////////////////////////
-  // !!! CHANGE VARIABLES NAME !!!                                        //
-  //                                                                      //
-  // perm              = nodeReducedGlobalId                              //
-  // refEntity         = entityNodeIndex                                  //
-  // entity            = orientedEntityNodeIndex                          //
-  //                                                                      //
-  // !!! I'm Working on *BOTH* indexes and global IDs !!!                 //
-  //                                                                      //
-  // 16 +                                                                 //
-  //    |\        ---> NodeIndex       = [v0 v1 v2]                       //
-  //    | \       ---> GlobalID        = [17 42 16]                       //
-  //    |  \      ---> ReducedGlobalId = [1  2  0 ] ---> Permutation = 3  //
-  // 17 +---+ 42                                                          //
-  //                                                                      //
-  // *  refEdge[e][i] = ith INDEX of edge[e]                              //
-  //    For example refEdge[1] = [1 2] because edge[1]                    //
-  //                              - -                                     //
-  //                                           is made of nodes v1 and v2 //
-  //                                                             -      - //
-  //                                                                      //
-  // *  edge[p][e][i] = ith INDEX of edge[e] in orientation[p]            //
-  //          such that GlobalID[edge[p][e][0]] > GlobalId[edge[p][e][0]] //
-  //    For example edge[3][1] = [2 1] because edge[1] goes from v2 to v1 //
-  //                              - -                             -     - //
-  //////////////////////////////////////////////////////////////////////////
-
-  // Number of Vertices //
-  size_t nVertex;
+  // Reference Spaces Node Ids //
+  std::vector<std::vector<size_t> > refSpaceNodeId;
 
-  // Edge Permutation //
-  size_t    nEdge;
-  size_t**  refEdge;
-  //unsigned int*** permutedRefEdge;
-  std::vector<const std::vector<const std::vector<unsigned int>*>*>* edge;
+  // ABC space to UVW space shape function index mapping //
+  std::vector<std::vector<size_t> > ABCtoUVWIndex;
 
-  // Face Permutation //
-  size_t    nFace;
-  size_t*   nNodeInFace;
-  size_t**  refFace;
-  //unsigned int*** permutedRefFace;
-  std::vector<const std::vector<const std::vector<unsigned int>*>*>* face;
+  // UVW space to ABC space shape function index mapping //
+  std::vector<std::vector<size_t> > UVWtoABCIndex;
+
+  // Ordered Edge & Face Node Index //
+  std::vector<std::vector<std::vector<size_t> > > orderedEdgeNodeIdx;
+  std::vector<std::vector<std::vector<size_t> > > orderedFaceNodeIdx;
 
  public:
   virtual ~ReferenceSpace(void);
 
   size_t getNReferenceSpace(void) const;
-
   size_t getReferenceSpace(const MElement& element) const;
-  /*
-  std::vector<std::vector<size_t> > getEdgeIndex(size_t referenceSpaceId) const;
-  std::vector<std::vector<size_t> > getFaceIndex(size_t referenceSpaceId) const;
-  */
-  const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-    getAllEdge(void) const;
 
-  const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-    getAllFace(void) const;
+  const std::vector<std::vector<size_t> >&
+    getNodeId(void) const;
 
-  virtual std::string toString(void) const;
+  const std::vector<std::vector<std::vector<size_t> > >&
+    getEdgeNodeIndex(void) const;
+
+  const std::vector<std::vector<std::vector<size_t> > >&
+    getFaceNodeIndex(void) const;
+
+  void mapFromXYZtoABC(const MElement& element,
+                       const fullVector<double>& xyz,
+                       double abc[3]);
+
+  void getJacobian(const MElement& element,
+                   const fullVector<double>& xyz,
+                   fullMatrix<double>& jac) const;
 
+  virtual std::string toString(void) const;
   virtual std::string toLatex(void) const;
 
  protected:
@@ -99,38 +80,32 @@ class ReferenceSpace{
 
   void init(void);
 
-  void getEdge(void);
-  void getFace(void);
+  void getOrderedEdge(void);
+  void getOrderedFace(void);
+
+  void findCyclicPermutation
+    (std::list<size_t>&               listOfTrueReferenceSpace,
+     std::list<std::vector<size_t> >& listOfRefNodeIndexPermutation,
+     std::list<std::vector<size_t> >& listOfReverseNodeIndexPermutation);
 
-  void   findCyclicPermutation(void);
   triple isCyclicPermutation(std::vector<size_t>& pTest,
                              std::vector<size_t>& pRef);
 
   size_t findCorrespondingFace(std::vector<size_t>& face,
                                std::vector<size_t>& node);
 
-  /*
-  void getPermutedRefEntity(unsigned int**** permutedRefEntity,
-                            unsigned int**   refEntity,
-                            unsigned int*    nNodeInEntity,
-                            unsigned int     nEntity);
-  */
-  const std::vector<unsigned int>* getOrderedEdge(unsigned int permutation,
-                                                  unsigned int edge);
-
-  const std::vector<unsigned int>* getOrderedTriFace(unsigned int permutation,
-                                                     unsigned int face);
-
-  const std::vector<unsigned int>* getOrderedQuadFace(unsigned int permutation,
-                                                      unsigned int face);
+  static void
+    orderRefEntityForGivenRefSpace(std::vector<size_t>& refEntityNodeIdx,
+                                   std::vector<size_t>& refSpaceNodeId,
+                                   std::vector<size_t>& orderedEntityNodeIdx);
 
   static void
-    orderRefEntityForGivenPermutation(size_t* refEntity,
-                                      std::vector<size_t>& permutation,
-                                      std::vector<unsigned int>& orderedEntity);
+    correctQuadFaceNodeIdx(std::vector<size_t>& correctedQuadFaceNodeIdx);
+
+  size_t getPermutationIdx(const MElement& element) const;
 
-  static bool sortPredicate(const std::pair<unsigned int, unsigned int>& a,
-                            const std::pair<unsigned int, unsigned int>& b);
+  static bool sortPredicate(const std::pair<size_t, size_t>& a,
+                            const std::pair<size_t, size_t>& b);
 };
 
 
@@ -198,26 +173,30 @@ class ReferenceSpace{
 //////////////////////
 
 inline size_t ReferenceSpace::getNReferenceSpace(void) const{
-  return pTree->getNPermutation();
+  return refSpaceNodeId.size();
 }
 
 inline
-const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-ReferenceSpace::getAllEdge(void) const{
-  return *edge;
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpace::getEdgeNodeIndex(void) const{
+  return orderedEdgeNodeIdx;
 }
 
 inline
-const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
-ReferenceSpace::getAllFace(void) const{
-  return *face;
+const std::vector<std::vector<std::vector<size_t> > >&
+ReferenceSpace::getFaceNodeIndex(void) const{
+  return orderedFaceNodeIdx;
 }
 
+inline
+size_t ReferenceSpace::getReferenceSpace(const MElement& element) const{
+  return pTree->getTagFromPermutation(getPermutationIdx(element));
+}
 
 inline
 bool
-ReferenceSpace::sortPredicate(const std::pair<unsigned int, unsigned int>& a,
-                              const std::pair<unsigned int, unsigned int>& b){
+ReferenceSpace::sortPredicate(const std::pair<size_t, size_t>& a,
+                              const std::pair<size_t, size_t>& b){
   return a.second < b.second;
 }
 

FunctionSpace/ReferenceSpaceLagrange.cpp

@@ -21,6 +21,7 @@ ReferenceSpaceLagrange::~ReferenceSpaceLagrange(void){
 }
 
 void ReferenceSpaceLagrange::getLagrangeNode(void){
+  /*
   const size_t nPerm = pTree->getNPermutation();
 
   // Alloc //
@@ -58,6 +59,7 @@ void ReferenceSpaceLagrange::getLagrangeNode(void){
     // Insert in node
     (*node)[p] = tmp;
   }
+  */
 }
 
 void ReferenceSpaceLagrange::edgeSeq(vector<unsigned int>& vec,

FunctionSpace/TetEdgeBasis.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 
 TetEdgeBasis::TetEdgeBasis(unsigned int order){
+  /*
   // Reference Space //
   refSpace  = new TetReferenceSpace;
   nRefSpace = refSpace->getNReferenceSpace();
@@ -321,9 +322,11 @@ TetEdgeBasis::TetEdgeBasis(unsigned int order){
   delete[] legendre;
   delete[] sclLegendre;
   delete[] intLegendre;
+  */
 }
 
 TetEdgeBasis::~TetEdgeBasis(void){
+  /*
   // ReferenceSpace //
   delete refSpace;
 
@@ -336,4 +339,5 @@ TetEdgeBasis::~TetEdgeBasis(void){
   }
 
   delete[] basis;
+  */
 }