diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index b1e751a81417953a18de7b21ef30ba1c535d98d3..f7d73285d618d0ad7ad3442bd14a1d873cfe9b31 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -5,6 +5,7 @@
 //
 // Contributor(s):
 //   Koen Hillewaert
+//   Jonathan Lambrechts
 //
 
 #include "GmshDefines.h"
@@ -12,6 +13,25 @@
 #include "polynomialBasis.h"
 #include "Gauss.h"
 
+#include "stdlib.h"
+static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int> &closureRef, polynomialBasis::clCont &closures) {
+  for (int i = 0; i <closures.size(); i++) {
+    printf("%3i  (%2i): ", i, closureRef[i]);
+    if(closureRef[i]==-1){
+      printf("\n");
+      continue;
+    }
+    for (int j = 0; j < fullClosure[i].size(); j++) {
+      printf("%i ", fullClosure[i][j]);
+    }
+    printf ("  --  ");
+    for (int j = 0; j < closures[closureRef[i]].size(); j++) {
+      printf("%i-%i ", fullClosure[i][closures[closureRef[i]][j]], closures[i][j]);
+    }
+    printf("\n");
+  }
+}
+
 static int getTriangleType (int order) {
   switch(order) {
     case 0 : return MSH_TRI_1;
@@ -861,6 +881,33 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
   return coefficient;
 }
 
+static void generate1dVertexClosure(polynomialBasis::clCont &closure)
+{
+  closure.clear();
+  closure.resize(2);
+  closure[0].push_back(0);
+  closure[1].push_back(1);
+  closure[0].type = MSH_PNT;
+  closure[1].type = MSH_PNT;
+}
+
+static void generate1dVertexClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order)
+{
+  closure.clear();
+  closure.resize(2);
+  closure[0].push_back(0);
+  closure[0].push_back(1);
+  closure[1].push_back(1);
+  closure[1].push_back(0);
+  for (int i = 0; i < order - 1; i++) {
+    closure[0].push_back(2 + i);
+    closure[1].push_back(2 + order - 2 - i);
+  }
+  closureRef.resize(2);
+  closureRef[0] = 0;
+  closureRef[1] = 0;
+}
+
 static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis::closure &closure,
                            int order)
 {
@@ -927,10 +974,130 @@ static void getFaceClosureTet(int iFace, int iSign, int iRotate, polynomialBasis
     break;
   }
 }
+static void generate2dEdgeClosureFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order, int nNod, bool serendip)
+{
+  closure.clear();
+  closure.resize(2*nNod);
+  closureRef.resize(2*nNod);
+  int shift = 0;
+  for (int corder = order; corder>=0; corder -=3) {
+    if (corder == 0) {
+      for (int r = 0; r < nNod ; r++){ 
+        closure[r].push_back(shift);
+        closure[r+nNod].push_back(shift);
+      }
+      break;
+    }
+    for (int r = 0; r < nNod ; r++){ 
+      for (int j = 0; j < nNod; j++) {
+        closure[r].push_back(shift + (r + j) % nNod);
+        closure[r + nNod].push_back(shift + (r - j + 1 + nNod) % nNod);
+      }
+    }
+    shift += nNod;
+    int n = nNod*(corder-1);
+    for (int r = 0; r < nNod ; r++){ 
+      for (int j = 0; j < n; j++){
+        closure[r].push_back(shift + (j + (corder - 1) * r) % n);
+        closure[r + nNod].push_back(shift + (n - j - 1 + (corder - 1) * (r + 1)) % n);
+      }
+    }
+    shift += n;
+    if (serendip) break;
+  }
+  for (int r = 0; r < nNod*2 ; r++) {
+    closure[r].type = getLineType(order);
+    closureRef[r] = 0;
+  }
+}
 
-static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
+static void generateFaceClosureTetFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order, bool serendip)
 {
+  closure.clear();
+  //input :
+  static const short int faces[4][3] = {{0,1,2}, {0,3,1}, {0,2,3}, {3,2,1}};
+  static const short int edges[6][2] = {{0,1}, {1,2}, {2,0}, {3,0}, {3,2}, {3,1}};
+  static const int faceOrientation[6] = {0,1,2,5,3,4};
+  closure.resize(24);
+  closureRef.resize(24);
+  for (int i = 0; i < 24; i ++)
+    closureRef[i] = 0;
+  if (order == 0) {
+    for (int i = 0; i < 24; i ++) {
+      closure[i].push_back(0);
+    }
+    return;
+  }
+  //Mapping for the p1 nodes
+  for (int iRotate = 0; iRotate < 3; iRotate ++) {
+    for (int iFace = 0; iFace < 4; iFace ++) {
+      for (int iNode = 0; iNode < 3; iNode ++) {
+        closure[iFace + iRotate * 8 + 0].push_back(faces[iFace][(6 + iNode - iRotate) % 3]);
+        closure[iFace + iRotate * 8 + 4].push_back(faces[iFace][(6 - iNode - iRotate) % 3]);
+      }
+    }
+  }
+  short int nodes2Edges [4][4];
+  for (int i = 0; i < 6; i++) {
+    nodes2Edges[edges[i][0]][edges[i][1]] = i * 2;
+    nodes2Edges[edges[i][1]][edges[i][0]] = i * 2 + 1;
+  }
+  int nodes2Faces[4][4][4];
+  for (int i = 0; i < 4; i++) {
+    for (int iRotate = 0; iRotate < 3; iRotate ++) {
+      short int n0 = faces[i][(3 - iRotate) % 3];
+      short int n1 = faces[i][(4 - iRotate) % 3];
+      short int n2 = faces[i][(5 - iRotate) % 3];
+      nodes2Faces[n0][n1][n2] = i * 6 + iRotate;
+      nodes2Faces[n0][n2][n1] = i * 6 + iRotate + 3;
+    }
+  }
+  polynomialBasis::clCont closureTriangles;
+  std::vector<int> closureTrianglesRef;
+  if (order >= 3)
+    generate2dEdgeClosureFull(closureTriangles, closureTrianglesRef, order - 3, 3, false); 
+  for (int iClosure = 0; iClosure < closure.size(); iClosure++) {
+    std::vector<int> &cl = closure[iClosure];
+    //last node
+    cl.push_back(6 - cl[0] - cl[1] - cl[2]);
+    //edges
+    for (int iEdge = 0; iEdge < 6; iEdge ++) {
+      int n0 = cl[edges[iEdge][0]];
+      int n1 = cl[edges[iEdge][1]];
+      short int &oEdge = nodes2Edges[n0][n1];
+      for (int i = 0 ; i < order - 1; i++)
+        cl.push_back(4 + oEdge/2 * (order - 1) + ((oEdge % 2) ?  order - 2 - i : i));
+    }
+    //faces
+    if (order >= 3) {
+      for (int iFace = 0; iFace < 4; iFace ++) {
+        int n0 = cl[faces[iFace][0]];
+        int n1 = cl[faces[iFace][1]];
+        int n2 = cl[faces[iFace][2]];
+        short int id = nodes2Faces[n0][n1][n2];
+        short int iTriClosure = faceOrientation [id % 6];
+        short int idFace = id/6;
+        int nOnFace = closureTriangles[iTriClosure].size();
+        for (int i = 0; i < nOnFace; i++) {
+          cl.push_back(4 + 6 * (order - 1) + idFace * nOnFace + closureTriangles[iTriClosure][i]);
+        }
+      }
+    }
+  }
+  if (order >= 4 && !serendip) {
+    polynomialBasis::clCont insideClosure;
+    std::vector<int> fakeClosureRef;
+    generateFaceClosureTetFull(insideClosure, fakeClosureRef, order - 4, false);
+    for (int i = 0; i < closure.size(); i ++) {
+      int shift = closure[i].size();
+      for (int j = 0; j < insideClosure[i].size(); j++)
+        closure[i].push_back(insideClosure[i][j] + shift);
+    }
+  }
+}
 
+static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
+{
   closure.clear();
   for (int iRotate = 0; iRotate < 3; iRotate++){
     for (int iSign = 1; iSign >= -1; iSign -= 2){
@@ -942,6 +1109,7 @@ static void generateFaceClosureTet(polynomialBasis::clCont &closure, int order)
     }
   }
 }
+
 static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
                                 polynomialBasis::closure &closure, int order)
 {
@@ -977,6 +1145,48 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate,
       closure[nNodes-1] = order2node[iFace][4]; // center
   }
 }
+static void generateFaceClosurePrismFull(polynomialBasis::clCont &closure, std::vector<int> &closureRef, int order) {
+  closure.clear();
+  //input :
+  static const int triRotation[3][6] = {{0,1,2,3,4,5},{1,2,0,4,5,3},{2,0,1,5,3,4}};
+  static const int triSign[2][6] = {{0,1,2,3,4,5},{0,2,1,3,5,4}};
+  static const int triFace[5][6] = {{0,2,1,3,5,4},{3,4,5,0,1,2}};
+  static const int quadFace[3][6] = {{0,1,2,3,4,5},{5,2,4,3,0,1},{1,2,0,4,5,3}};
+  static const int quadRotation[4][6] = {{0,1,2,3,4,5},{3,0,5,4,1,2},{4,3,5,1,0,2},{1,4,5,0,3,2}};
+  static const int quadSign[2][6] = {{0,1,2,3,4,5},{4,1,5,3,0,2}};
+  //static const short int edges[6][2] = {{0,1}, {1,2}, {2,0}, {3,0}, {3,2}, {3,1}};
+  closure.resize(40);
+  closureRef.resize(40);
+  if (order == 0) {
+    for (int i = 0; i < 40; i ++) {
+      closure[i].push_back(0);
+    }
+    return;
+  }
+  for (int i = 0; i < 40; i++)
+    closure[i].resize(6);
+  //Mapping for the p1 nodes
+  //triangular faces
+  for (int iRotate = 0; iRotate < 3; iRotate ++)
+    for (int iFace = 0; iFace < 2; iFace ++)
+      for (int iSign = 0; iSign < 2; iSign ++)
+        for (int iNode = 0; iNode < 6; iNode ++) {
+          int i = iFace + iRotate * 10 + iSign * 5;
+          closureRef[i] = 0;
+          closure[i][triFace[0][iNode]] =
+            triFace[iFace][triRotation[iRotate][triSign[iSign][iNode]]];
+        }
+  //quad Horiz face
+  for (int iRotate = 0; iRotate < 4; iRotate ++)
+    for (int iFace = 0; iFace < 3; iFace ++)
+      for (int iSign = 0; iSign < 2; iSign ++)
+        for (int iNode = 0; iNode < 6; iNode ++) {
+          int i = (2 + iFace) + iRotate * 10 + iSign * 5;
+          closureRef[i] = (iFace + iRotate + iSign) % 2 ? 3 : 2;
+          closure[i][quadFace[closureRef[i]-2][quadRotation[0][iNode]]] =
+            quadFace[iFace][quadRotation[iRotate][quadSign[iSign][iNode]]];
+        }
+}
 
 static void generateFaceClosurePrism(polynomialBasis::clCont &closure, int order)
 {
@@ -993,8 +1203,7 @@ static void generateFaceClosurePrism(polynomialBasis::clCont &closure, int order
   }
 }
 
-static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order,
-                                  int nNod = 3)
+static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
 {
   closure.clear();
   closure.resize(2*nNod);
@@ -1011,16 +1220,6 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order,
   }
 }
 
-static void generate1dVertexClosure(polynomialBasis::clCont &closure)
-{
-  closure.clear();
-  closure.resize(2);
-  closure[0].push_back(0);
-  closure[1].push_back(1);
-  closure[0].type = MSH_PNT;
-  closure[1].type = MSH_PNT;
-}
-
 static void generateClosureOrder0(polynomialBasis::clCont &closure, int nb)
 {
   closure.clear();
@@ -1039,7 +1238,6 @@ const polynomialBasis *polynomialBases::find(int tag)
   if (it != fs.end())     return &it->second;
   polynomialBasis F;
   F.numFaces = -1;
-
   switch (tag){
   case MSH_PNT:
     F.numFaces = 1;
@@ -1052,6 +1250,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.monomials = generate1DMonomials(0);
     F.points    = generate1DPoints(0);
     generateClosureOrder0(F.closures,2);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 0);
     F.parentType = TYPE_LIN;
     break;
   case MSH_LIN_2 :
@@ -1059,6 +1258,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.monomials = generate1DMonomials(1);
     F.points    = generate1DPoints(1);
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 1);
     F.parentType = TYPE_LIN;
     break;
   case MSH_LIN_3 :
@@ -1066,6 +1266,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.monomials = generate1DMonomials(2);
     F.points    = generate1DPoints(2);
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 2);
     F.parentType = TYPE_LIN;
     break;
   case MSH_LIN_4:
@@ -1074,6 +1275,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(3);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 3);
     break;
   case MSH_LIN_5:
     F.numFaces = 2;
@@ -1081,6 +1283,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(4);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 4);
     break;
   case MSH_LIN_6:
     F.numFaces = 2;
@@ -1088,6 +1291,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(5);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 5);
     break;
   case MSH_LIN_7:
     F.numFaces = 2;
@@ -1095,6 +1299,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(6);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 6);
     break;
   case MSH_LIN_8:
     F.numFaces = 2;
@@ -1102,6 +1307,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(7);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 7);
     break;
   case MSH_LIN_9:
     F.numFaces = 2;
@@ -1109,6 +1315,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(8);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 8);
     break;
   case MSH_LIN_10:
     F.numFaces = 2;
@@ -1116,6 +1323,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(9);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 9);
     break;
   case MSH_LIN_11:
     F.numFaces = 2;
@@ -1123,13 +1331,16 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points    = generate1DPoints(10);
     F.parentType = TYPE_LIN;
     generate1dVertexClosure(F.closures);
+    generate1dVertexClosureFull(F.fullClosures, F.closureRef, 10);
     break;
   case MSH_TRI_1 :
     F.numFaces = 3;
     F.monomials = generatePascalTriangle(0);
     F.points =    gmshGeneratePointsTriangle(0, false);
     F.parentType = TYPE_TRI;
-    generateClosureOrder0(F.closures,6);
+    generateClosureOrder0(F.closures, 6);
+    generateClosureOrder0(F.fullClosures, 6);
+    F.closureRef.resize(6, 0);
     break;
   case MSH_TRI_3 :
     F.numFaces = 3;
@@ -1137,6 +1348,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(1, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 1);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 1, 3, false);
     break;
   case MSH_TRI_6 :
     F.numFaces = 3;
@@ -1144,6 +1356,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(2, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 2);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 3, false);
     break;
   case MSH_TRI_9 :
     F.numFaces = 3;
@@ -1151,6 +1364,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(3, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 3);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 3, true);
     break;
   case MSH_TRI_10 :
     F.numFaces = 3;
@@ -1158,6 +1372,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(3, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 3);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 3, false);
     break;
   case MSH_TRI_12 :
     F.numFaces = 3;
@@ -1165,6 +1380,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(4, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 3, true);
     break;
   case MSH_TRI_15 :
     F.numFaces = 3;
@@ -1172,6 +1388,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(4, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 3, false);
     break;
   case MSH_TRI_15I :
     F.numFaces = 3;
@@ -1179,6 +1396,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(5, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 5);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 3, true);
     break;
   case MSH_TRI_21 :
     F.numFaces = 3;
@@ -1186,6 +1404,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(5, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 5);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 3, false);
     break;
   case MSH_TRI_28 :
     F.numFaces = 3;
@@ -1193,6 +1412,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(6, false);    
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 6);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 3, false);
     break;
   case MSH_TRI_36 :
     F.numFaces = 3;
@@ -1200,6 +1420,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(7, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 7);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 3, false);
     break;
   case MSH_TRI_45 :
     F.numFaces = 3;
@@ -1207,6 +1428,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(8, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 8);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 3, false);
     break;
   case MSH_TRI_55 :
     F.numFaces = 3;
@@ -1214,6 +1436,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(9, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 9);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 3, false);
     break;
   case MSH_TRI_66 :
     F.numFaces = 3;
@@ -1221,6 +1444,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(10, false);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 10);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 3, false);
     break;
   case MSH_TRI_18 :
     F.numFaces = 3;
@@ -1228,6 +1452,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(6, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 6);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 3, true);
     break;
   case MSH_TRI_21I :
     F.numFaces = 3;
@@ -1235,6 +1460,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(7, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 7);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 3, true);
     break;
   case MSH_TRI_24 :
     F.numFaces = 3;
@@ -1242,6 +1468,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(8, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 8);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 3, true);
     break;
   case MSH_TRI_27 :
     F.numFaces = 3;
@@ -1249,6 +1476,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(9, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 9);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 3, true);
     break;
   case MSH_TRI_30 :
     F.numFaces = 3;
@@ -1256,6 +1484,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTriangle(10, true);
     F.parentType = TYPE_TRI;
     generate2dEdgeClosure(F.closures, 10);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 3, true);
     break;
   case MSH_TET_1 :
     F.numFaces = 4;
@@ -1263,6 +1492,8 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(0, false);
     F.parentType = TYPE_TET;
     generateClosureOrder0(F.closures,24);
+    generateClosureOrder0(F.fullClosures, 24);
+    F.closureRef.resize(24, 0);
     break;
   case MSH_TET_4 :
     F.numFaces = 4;
@@ -1270,6 +1501,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(1, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 1);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 1, false);
     break;
   case MSH_TET_10 :
     F.numFaces = 4;
@@ -1277,6 +1509,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(2, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 2);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 2, false);
     break;
   case MSH_TET_20 :
     F.numFaces = 4;
@@ -1284,6 +1517,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(3, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 3);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 3, false);
     break;
   case MSH_TET_35 :
     F.numFaces = 4;
@@ -1291,6 +1525,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(4, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 4);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 4, false);
     break;
   case MSH_TET_34 :
     F.numFaces = 4;
@@ -1298,6 +1533,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(4, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 4);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 4, true);
     break;
   case MSH_TET_52 :
     F.numFaces = 4;
@@ -1305,6 +1541,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(5, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 5);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 5, true);
     break;
   case MSH_TET_56 :
     F.numFaces = 4;
@@ -1312,6 +1549,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(5, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 5);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 5, false);
     break;
   case MSH_TET_74 :
     F.numFaces = 4;
@@ -1319,6 +1557,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(6, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 6);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 6, true);
     break;
   case MSH_TET_84 :
     F.numFaces = 4;
@@ -1326,6 +1565,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(6, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 6);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 6, false);
     break;
   case MSH_TET_100 :
     F.numFaces = 4;
@@ -1333,6 +1573,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(7, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 7);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 7, true);
     break;
   case MSH_TET_120 :
     F.numFaces = 4;
@@ -1340,6 +1581,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(7, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 7);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 7, false);
     break;
   case MSH_TET_130 :
     F.numFaces = 4;
@@ -1347,6 +1589,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(8, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 8);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 8, true);
     break;
   case MSH_TET_164 :
     F.numFaces = 4;
@@ -1354,6 +1597,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(9, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 9);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 9, true);
     break;
   case MSH_TET_165 :
     F.numFaces = 4;
@@ -1361,6 +1605,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(8, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 8);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 8, false);
     break;
   case MSH_TET_202 :
     F.numFaces = 4;
@@ -1368,6 +1613,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(10, true);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 10);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 10, true);
     break;
   case MSH_TET_220 :
     F.numFaces = 4;
@@ -1375,6 +1621,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(9, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 9);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 9, false);
     break;
   case MSH_TET_286 :
     F.numFaces = 4;
@@ -1382,6 +1629,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(10, false);
     F.parentType = TYPE_TET;
     generateFaceClosureTet(F.closures, 10);
+    generateFaceClosureTetFull(F.fullClosures, F.closureRef, 10, false);
     break;
   case MSH_QUA_1 :
     F.numFaces = 4;
@@ -1389,6 +1637,8 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(0,false);
     F.parentType = TYPE_QUA;
     generateClosureOrder0(F.closures,8);
+    generateClosureOrder0(F.fullClosures,8);
+    F.closureRef.resize(8, 0);
     break;
   case MSH_QUA_4 :
     F.numFaces = 4;
@@ -1396,6 +1646,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(1,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 1, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 1, 4, false);
     break;
   case MSH_QUA_9 :
     F.numFaces = 4;
@@ -1403,6 +1654,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(2,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 2, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 4, false);
     break;
   case MSH_QUA_16 :
     F.numFaces = 4;
@@ -1410,6 +1662,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(3,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 3, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 4, false);
     break;
   case MSH_QUA_25 :
     F.numFaces = 4;
@@ -1417,6 +1670,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(4,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 4, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 4, false);
     break;
   case MSH_QUA_36 :
     F.numFaces = 4;
@@ -1424,6 +1678,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(5,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 5, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 4, false);
     break;
   case MSH_QUA_49 :
     F.numFaces = 4;
@@ -1431,6 +1686,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(6,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 6, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 4, false);
     break;
   case MSH_QUA_64 :
     F.numFaces = 4;
@@ -1438,6 +1694,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(7,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 7, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 4, false);
     break;
   case MSH_QUA_81 :
     F.numFaces = 4;
@@ -1445,6 +1702,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(8,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 8, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 4, false);
     break;
   case MSH_QUA_100 :
     F.numFaces = 4;
@@ -1452,6 +1710,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(9,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 9, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 4, false);
     break;
   case MSH_QUA_121 :
     F.numFaces = 4;
@@ -1459,6 +1718,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(10,false);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 10, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 4, false);
     break;
   case MSH_QUA_8 :
     F.numFaces = 4;
@@ -1466,6 +1726,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(2,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 2, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 2, 4, true);
     break;
   case MSH_QUA_12 :
     F.numFaces = 4;
@@ -1473,6 +1734,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(3,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 3, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 3, 4, true);
     break;
   case MSH_QUA_16I :
     F.numFaces = 4;
@@ -1480,6 +1742,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(4,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 4, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 4, 4, true);
     break;
   case MSH_QUA_20 :
     F.numFaces = 4;
@@ -1487,6 +1750,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(5,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 5, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 5, 4, true);
     break;
   case MSH_QUA_24 :
     F.numFaces = 4;
@@ -1494,6 +1758,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(6,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 6, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 6, 4, true);
     break;
   case MSH_QUA_28 :
     F.numFaces = 4;
@@ -1501,6 +1766,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(7,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 7, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 7, 4, true);
     break;
   case MSH_QUA_32 :
     F.numFaces = 4;
@@ -1508,6 +1774,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(8,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 8, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 8, 4, true);
     break;
   case MSH_QUA_36I :
     F.numFaces = 4;
@@ -1515,6 +1782,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(9,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 9, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 9, 4, true);
     break;
   case MSH_QUA_40 :
     F.numFaces = 4;
@@ -1522,6 +1790,7 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsQuad(10,true);
     F.parentType = TYPE_QUA;
     generate2dEdgeClosure(F.closures, 10, 4);
+    generate2dEdgeClosureFull(F.fullClosures, F.closureRef, 10, 4, true);
     break;
   case MSH_PRI_1 :
     F.numFaces = 5;
@@ -1529,6 +1798,8 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsPrism(0, false);
     F.parentType = TYPE_PRI;
     generateClosureOrder0(F.closures,48);
+    generateClosureOrder0(F.fullClosures,48);
+    F.closureRef.resize(48, 0);
     break;
   case MSH_PRI_6 :
     F.numFaces = 5;
@@ -1536,14 +1807,15 @@ const polynomialBasis *polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsPrism(1, false);
     F.parentType = TYPE_PRI;
     generateFaceClosurePrism(F.closures, 1);
+    generateFaceClosurePrismFull(F.fullClosures, F.closureRef, 1);
     break;
-  case MSH_PRI_18 :
+  /*case MSH_PRI_18 :
     F.numFaces = 5;
     F.monomials = generatePascalPrism(2);
     F.points =    gmshGeneratePointsPrism(2, false);
     F.parentType = TYPE_PRI;
     generateFaceClosurePrism(F.closures, 2);
-    break;
+    break;*/
   default :
     Msg::Error("Unknown function space %d: reverting to TET_4", tag);
     F.numFaces = 4;
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 566a224ad596858f787a167c35b19d92ffc78261..6a05e2add178dc40a97d0dda57ab082c0d49299c 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -80,7 +80,13 @@ class polynomialBasis
     int type;
   };
   typedef std::vector<closure> clCont;
-  clCont closures;
+  // closures is the list of the nodes of each face, in the order of the polynomialBasis of the face
+  // fullClosures is  mapping of the nodes of the element that rotates the element so that the considered face becomes the first one
+  // in the right orientation
+  // For element, like prisms that have different kind of faces, fullCLosure[i] rotates the element so that the considered face
+  // becomes the closureRef[i]-th face (the first tringle or the first quad face)
+   clCont closures, fullClosures;
+  std::vector<int> closureRef;
   fullMatrix<double> points;
   fullMatrix<double> monomials;
   fullMatrix<double> coefficients;