gmsh + dg : mini tet

Common/GmshDefines.h

@@ -212,8 +212,9 @@
 #define MSH_TET_SUB 136
 #define MSH_TET_16  137
 #define MSH_TRI_MINI 138
+#define MSH_TET_MINI 139
 
-#define MSH_NUM_TYPE 138
+#define MSH_NUM_TYPE 139
 
 // Geometric entities
 #define ENT_NONE     0

Numeric/BasisFactory.cpp

@@ -30,9 +30,10 @@ const nodalBasis* BasisFactory::getNodalBasis(int tag)
   // Get the parent type to see which kind of basis
   // we want to create
   nodalBasis* F = NULL;
-  if (tag == MSH_TRI_MINI) {
-    F = new miniBasis();
-  }
+  if (tag == MSH_TRI_MINI)
+    F = new miniBasisTri();
+  else if (tag == MSH_TET_MINI)
+    F = new miniBasisTet();
   else {
     int parentType = ElementType::ParentTypeFromTag(tag);
     switch(parentType) {

Numeric/ElementType.cpp

@@ -91,6 +91,7 @@ int ElementType::ParentTypeFromTag(int tag)
     case(MSH_TRI_SUB):  case(MSH_TET_SUB):
       return TYPE_XFEM;
     case(MSH_TRI_MINI):
+    case(MSH_TET_MINI):
       return TYPE_MINI;
     default:
       Msg::Error("Unknown element tag %i for parent type, returning -1.", tag);
@@ -132,6 +133,7 @@ int ElementType::OrderFromTag(int tag)
   case MSH_TRI_24  : return 8;
   case MSH_TRI_27  : return 9;
   case MSH_TRI_30  : return 10;
+  case MSH_TRI_MINI  : return 1;
   case MSH_POLYG_  : return 1;
   case MSH_POLYG_B : return 1;
   case MSH_TET_1   : return 0;
@@ -153,6 +155,7 @@ int ElementType::OrderFromTag(int tag)
   case MSH_TET_46  : return 8;
   case MSH_TET_52  : return 9;
   case MSH_TET_58  : return 10;
+  case MSH_TET_MINI  : return 1;
   case MSH_QUA_1   : return 0;
   case MSH_QUA_4   : return 1;
   case MSH_QUA_9   : return 2;
@@ -259,7 +262,7 @@ int ElementType::DimensionFromTag(int tag)
     case(MSH_TRI_21I):  case(MSH_TRI_24):
     case(MSH_TRI_27):   case(MSH_TRI_30):
     case(MSH_TRI_B):    case(MSH_TRI_1):
-    case(MSH_TRI_SUB):
+    case(MSH_TRI_SUB):  case(MSH_TRI_MINI):
 
     case(MSH_QUA_4):    case(MSH_QUA_9):
     case(MSH_QUA_8):    case(MSH_QUA_16):
@@ -285,6 +288,7 @@ int ElementType::DimensionFromTag(int tag)
     case(MSH_TET_46):   case(MSH_TET_52):
     case(MSH_TET_58):   case(MSH_TET_1):
     case(MSH_TET_16):   case(MSH_TET_SUB):
+    case(MSH_TET_MINI):
 
     case(MSH_PYR_5):    case(MSH_PYR_14):
     case(MSH_PYR_13):   case(MSH_PYR_30):

Numeric/miniBasis.cpp

 #include "miniBasis.h"
 #include "BasisFactory.h"
-miniBasis::miniBasis()
+miniBasisTri::miniBasisTri()
 {
   type = MSH_TRI_MINI;
   parentType = TYPE_TRI;
@@ -21,6 +21,20 @@ miniBasis::miniBasis()
   points(1, 0) = 1.; points(1, 1) = 0.;
   points(2, 0) = 0.; points(2, 1) = 1.;
   points(3, 0) = 1./3; points(3, 1) = 1./3;
+  /*
+  monomials :
+    0 : 1
+    1 : x
+    2 : y
+    3 : xy
+    4 : x²y
+    5 : xy²
+  fs :
+    0 : 1 -1 -1  0  0  0 : 1 - x - y
+    1 : 1  0  0  0  0  0 : x
+    2 : 0  1  0  0  0  0 : y
+    3 : 0  0  0  1 -1 -1 : xy - x²y - xy² = (1 - x - y) x y
+  */
   monomials.resize(6, 2);
   monomials(0, 0) = 0.; monomials(0, 1) = 0.;
   monomials(1, 0) = 1.; monomials(1, 1) = 0.;
@@ -35,3 +49,59 @@ miniBasis::miniBasis()
   coefficients(2, 2) = 1.;
   coefficients(3, 3) = 1.; coefficients(3, 4) = -1.; coefficients(3, 5) = -1.; 
 }
+
+miniBasisTet::miniBasisTet()
+{
+  type = MSH_TET_MINI;
+  parentType = TYPE_TET;
+  order = 3;
+  dimension = 3;
+  numFaces = 4;
+  serendip = false;
+  const nodalBasis &p1 = *BasisFactory::getNodalBasis(MSH_TET_4);
+  closures = p1.closures;
+  fullClosures = p1.fullClosures;
+  for(size_t i = 0; i < fullClosures.size(); ++i) {
+    fullClosures[i].push_back(4);
+  }
+  closureRef = p1.closureRef;
+  points.resize(5, 3);
+  for (int i = 0; i < 4; ++i)
+    for (int j = 0; j < 3; ++j)
+      points(i, j) = p1.points(i, j);
+  points(4, 0) = points(4, 1) = points(4, 2) = 1./4;
+
+  /*
+  monomials :
+    0 : 1
+    1 : x
+    2 : y
+    3 : z
+    4 : xyz
+    5 : x²yz
+    6 : xy²z
+    7 : xyz²
+  sf :
+    0 : 1 -1 -1 -1  0  0  0  0 : 1 - x - y - z
+    1 : 0  1  0  0  0  0  0  0 : x
+    2 : 0  0  1  0  0  0  0  0 : y
+    3 : 0  0  0  1  0  0  0  0 : z
+    4 : 0  0  0  0  1 -1 -1 -1 : xyz - x²yz - xy²z - xyz² = (1 - x - y - z) x y z
+  */
+  monomials.resize(8, 3);
+  monomials(0, 0) = 0.; monomials(0, 1) = 0.; monomials(0, 2) = 0.;
+  monomials(1, 0) = 1.; monomials(1, 1) = 0.; monomials(1, 2) = 0.;
+  monomials(2, 0) = 0.; monomials(2, 1) = 1.; monomials(2, 2) = 0.;
+  monomials(3, 0) = 0.; monomials(3, 1) = 0.; monomials(3, 2) = 1.;
+  monomials(4, 0) = 1.; monomials(4, 1) = 1.; monomials(4, 2) = 1.;
+  monomials(5, 0) = 2.; monomials(5, 1) = 1.; monomials(5, 2) = 1.;
+  monomials(6, 0) = 1.; monomials(6, 1) = 2.; monomials(6, 2) = 1.;
+  monomials(7, 0) = 1.; monomials(7, 1) = 1.; monomials(7, 2) = 2.;
+  coefficients.resize(5, 8);
+  coefficients.setAll(0.);
+  coefficients(0, 0) = 1.; coefficients(0, 1) = -1.; coefficients(0, 2) = -1.; coefficients(0, 3) = -1.;
+  coefficients(1, 1) = 1.;
+  coefficients(2, 2) = 1.;
+  coefficients(3, 3) = 1.;
+  coefficients(4, 4) = 1.; coefficients(4, 5) = -1.; coefficients(4, 6) = -1.; coefficients(4, 7) = -1.;
+}

Numeric/miniBasis.h

 #ifndef _MINI_BASIS_H_
 #define _MINI_BASIS_H_
 #include "polynomialBasis.h"
-class miniBasis : public polynomialBasis { // mini is NOT a real nodal basis but in GMSH, only the nodal basis have closure and mini have closure so...
+ // mini is NOT a real nodal basis but in GMSH, only the nodal basis have closure and mini have closure so...
+class miniBasisTri : public polynomialBasis {
   public:
-    miniBasis();
+    miniBasisTri();
+};
+class miniBasisTet : public polynomialBasis {
+  public:
+    miniBasisTet();
 };
 #endif