dg : implicit 3d benchmark + order 6 to 10 tet polynomial basis

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

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);
     }
   }

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);

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));