diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 6b0e1e201aedf3828db09a3de913c98e230ff953..4859b4ef3c003245d39a2dfd57862bab3bdf776c 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -126,8 +126,13 @@
 #define MSH_TRI_B 68
 #define MSH_POLYG_B 69
 #define MSH_LIN_C 70
+#define MSH_TET_84 71
+#define MSH_TET_120 72
+#define MSH_TET_165 73
+#define MSH_TET_220 74
+#define MSH_TET_286 75
 
-#define MSH_NUM_TYPE 70
+#define MSH_NUM_TYPE 71
 
 // Geometric entities
 #define ENT_NONE     0
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index ea523860ec56ac2e3d5ac7871c55becdcbdb04fc..ca1cc08c550005fe23176be1e96f516697133403 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -126,6 +126,11 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
     case 3: return &polynomialBases::find(MSH_TET_20);
     case 4: return &polynomialBases::find(MSH_TET_35);
     case 5: return &polynomialBases::find(MSH_TET_56);
+    case 6: return &polynomialBases::find(MSH_TET_84);
+    case 7: return &polynomialBases::find(MSH_TET_120);
+    case 8: return &polynomialBases::find(MSH_TET_165);
+    case 9: return &polynomialBases::find(MSH_TET_220);
+    case 10: return &polynomialBases::find(MSH_TET_286);
     default: Msg::Error("Order %d tetrahedron function space not implemented", order);
     }
   }
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 51c1cdee03e286df4028035e0b38d6f3053b7d23..39b93f36d2f2c4fed1ba7d02eb7c2ece3dfeb67f 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1114,6 +1114,36 @@ const polynomialBasis &polynomialBases::find(int tag)
     F.points =    gmshGeneratePointsTetrahedron(5, false);
     generate3dFaceClosure(F.closures, 5);
     break;
+  case MSH_TET_84 :
+    F.numFaces = 4;
+    F.monomials = generatePascalTetrahedron(6);
+    F.points =    gmshGeneratePointsTetrahedron(6, false);
+    generate3dFaceClosure(F.closures, 6);
+    break;
+  case MSH_TET_120 :
+    F.numFaces = 4;
+    F.monomials = generatePascalTetrahedron(7);
+    F.points =    gmshGeneratePointsTetrahedron(7, false);
+    generate3dFaceClosure(F.closures, 7);
+    break;
+  case MSH_TET_165 :
+    F.numFaces = 4;
+    F.monomials = generatePascalTetrahedron(8);
+    F.points =    gmshGeneratePointsTetrahedron(8, false);
+    generate3dFaceClosure(F.closures, 8);
+    break;
+  case MSH_TET_220 :
+    F.numFaces = 4;
+    F.monomials = generatePascalTetrahedron(9);
+    F.points =    gmshGeneratePointsTetrahedron(9, false);
+    generate3dFaceClosure(F.closures, 9);
+    break;
+  case MSH_TET_286 :
+    F.numFaces = 4;
+    F.monomials = generatePascalTetrahedron(10);
+    F.points =    gmshGeneratePointsTetrahedron(10, false);
+    generate3dFaceClosure(F.closures, 10);
+    break;
   case MSH_QUA_4 :
     F.numFaces = 4;
     F.monomials = generatePascalQuad(1);
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 82805bf9a0e14b8f3920c6891cc8fe66a5d6301c..a2fdcdb777fea26972b1646a1594dd297c27e5e8 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -83,7 +83,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows)
   // override the default options with the ones from the option
   // database (if any)
   _try(MatSetFromOptions(_a));
-  _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 4, NULL)); //todo preAllocate off-diagonal part
+  _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 5, NULL)); //todo preAllocate off-diagonal part
   //_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
   _try(VecCreate(PETSC_COMM_WORLD, &_x));
   _try(VecSetSizes(_x, PETSC_DECIDE, nbRows * _blockSize));