read tet10, pyramid13, prism15 and hex20 from gmsh

libsrc/general/array.hpp

@@ -243,10 +243,11 @@ namespace netgen
     }
 
     explicit Array(size_t asize)
-      : FlatArray<T, BASE, TIND> (asize, new T[asize])
+      : FlatArray<T, BASE, TIND> (asize, asize ? new T[asize] : nullptr)
     {
-      allocsize = asize; 
+      allocsize = asize;
-      ownmem = 1;
+      if(asize)
+        ownmem = 1;
     }
 
     /// Generate array in user data

libsrc/include/nginterface.h

@@ -32,7 +32,7 @@
 
 
 // max number of nodes per element
-#define NG_ELEMENT_MAXPOINTS 12
+#define NG_ELEMENT_MAXPOINTS 20
 
 // max number of nodes per surface element
 #define NG_SURFACE_ELEMENT_MAXPOINTS 8
@@ -49,7 +49,7 @@ enum NG_ELEMENT_TYPE {
   NG_SEGM = 1, NG_SEGM3 = 2,
   NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14,
   NG_TET = 20, NG_TET10 = 21, 
-  NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24,
+  NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28,
   NG_HEX = 25, NG_HEX20 = 26
 };
 

libsrc/include/nginterface_v2.hpp

@@ -19,7 +19,7 @@ enum NG_ELEMENT_TYPE {
   NG_SEGM = 1, NG_SEGM3 = 2,
   NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14,
   NG_TET = 20, NG_TET10 = 21, 
-  NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24,
+  NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28,
   NG_HEX = 25, NG_HEX20 = 26
 };
 

libsrc/meshing/curvedelems.cpp

@@ -2790,6 +2790,32 @@ namespace netgen
 	  break;
 	}
 
+      case PRISM15:
+        {
+	  shapes = 0.0;
+	  T x = xi(0);
+	  T y = xi(1);
+	  T z = xi(2);
+          T lam = 1-x-y;
+          T lamz = 1-z;
+          shapes[0] = (2*x*x-x) * (2*lamz*lamz-lamz);
+          shapes[1] = (2*y*y-y) * (2*lamz*lamz-lamz);
+          shapes[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz);
+          shapes[3] = (2*x*x-x) * (2*z*z-z);
+          shapes[4] = (2*y*y-y) * (2*z*z-z);
+          shapes[5] = (2*lam*lam-lam) * (2*z*z-z);
+          shapes[6] = 4 * x * y * (2*lamz*lamz-lamz);
+          shapes[7] = 4 * x * lam * (2*lamz*lamz-lamz);
+          shapes[8] = 4 * y * lam * (2*lamz*lamz-lamz);
+          shapes[9] = x * 4 * z * (1-z);
+          shapes[10] = y * 4 * z * (1-z);
+          shapes[11] = lam * 4 * z * (1-z);
+          shapes[12] = 4 * x * y * (2*z*z-z);
+          shapes[13] = 4 * x * lam * (2*z*z-z);
+          shapes[14] = 4 * y * lam * (2*z*z-z);
+          break;
+        }
+
       case PYRAMID:
 	{
 	  shapes = 0.0;
@@ -2839,6 +2865,29 @@ namespace netgen
 	  break;
 	}
 
+      case PYRAMID13:
+        {
+	  shapes = 0.0;
+	  T x = xi(0);
+	  T y = xi(1);
+	  T z = xi(2);
+          z *= 1-1e-12;
+          shapes[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25);
+          shapes[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2));
+          shapes[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75);
+          shapes[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2));
+          shapes[4] = z*(2*z - 1);
+          shapes[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          shapes[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2);
+          shapes[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          shapes[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2);
+          shapes[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1);
+          shapes[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1);
+          shapes[11] = 4*x*y*z/(-z + 1);
+          shapes[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1);
+          break;
+        }
+
       case HEX:
 	{
 	  shapes = 0.0;
@@ -3300,6 +3349,36 @@ namespace netgen
 
 	}
 
+      case PRISM15:
+        {
+          AutoDiff<3,T> x(xi(0), 0);
+          AutoDiff<3,T> y(xi(1), 1);
+          AutoDiff<3,T> z(xi(2), 2);
+          AutoDiff<3,T> ad[15];
+          AutoDiff<3,T> lam = 1-x-y;
+          AutoDiff<3,T> lamz = 1-z;
+
+          ad[0] = (2*x*x-x) * (2*lamz*lamz-lamz);
+          ad[1] = (2*y*y-y) * (2*lamz*lamz-lamz);
+          ad[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz);
+          ad[3] = (2*x*x-x) * (2*z*z-z);
+          ad[4] = (2*y*y-y) * (2*z*z-z);
+          ad[5] = (2*lam*lam-lam) * (2*z*z-z);
+          ad[6] = 4 * x * y * (2*lamz*lamz-lamz);
+          ad[7] = 4 * x * lam * (2*lamz*lamz-lamz);
+          ad[8] = 4 * y * lam * (2*lamz*lamz-lamz);
+          ad[9] = x * 4 * z * (1-z);
+          ad[10] = y * 4 * z * (1-z);
+          ad[11] = lam * 4 * z * (1-z);
+          ad[12] = 4 * x * y * (2*z*z-z);
+          ad[13] = 4 * x * lam * (2*z*z-z);
+          ad[14] = 4 * y * lam * (2*z*z-z);
+
+          for(int i=0; i<15; i++)
+            for(int j=0; j<3; j++)
+              dshapes(i,j) = ad[i].DValue(j);
+          break;
+        }
       case PYRAMID:
 	{
           // if (typeid(T) == typeid(SIMD<double>)) return;
@@ -3397,6 +3476,54 @@ namespace netgen
 	  break;
 	}
 
+      case PYRAMID13:
+        {
+	  T x = xi(0);
+	  T y = xi(1);
+	  T z = xi(2);
+          z *= 1-1e-12;
+          dshapes(0,0) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1) - 4);
+          dshapes(0,1) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1) - 4);
+          dshapes(0,2) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 5 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1));
+          dshapes(1,0) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(-2*y - z + 2) + (0.5*x - 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1) + 4);
+          dshapes(1,1) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*x + z)*(-2*y - z + 2) + (-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x - 0.5*y - 0.25);
+          dshapes(1,2) = (0.5*x - 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1));
+          dshapes(2,0) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75);
+          dshapes(2,1) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75);
+          dshapes(2,2) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (0.5*x + 0.5*y + 0.5*z - 0.75)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 1 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1));
+          dshapes(3,0) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1));
+          dshapes(3,1) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1) + 4);
+          dshapes(3,2) = (-0.5*x + 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1));
+          dshapes(4,0) = 0;
+          dshapes(4,1) = 0;
+          dshapes(4,2) = 4*z - 1;
+          dshapes(5,0) = -4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          dshapes(5,1) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2);
+          dshapes(5,2) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2) - 4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 4*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2));
+          dshapes(6,0) = -8*x*y/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)/(-2*z + 2);
+          dshapes(6,1) = 4*x*(-2*x - 2*z + 2)/(-2*z + 2);
+          dshapes(6,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*x - 2*z + 2)/((-2*z + 2)*(-2*z + 2));
+          dshapes(7,0) = -4*y*(-2*y - 2*z + 2)/(-2*z + 2);
+          dshapes(7,1) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          dshapes(7,2) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) - 4*y*(-2*y - 2*z + 2)/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2));
+          dshapes(8,0) = 4*y*(-2*y - 2*z + 2)/(-2*z + 2);
+          dshapes(8,1) = -8*x*y/(-2*z + 2) + 4*x*(-2*y - 2*z + 2)/(-2*z + 2);
+          dshapes(8,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2));
+          dshapes(9,0) = -2*z*(-2*y - 2*z + 2)/(-z + 1);
+          dshapes(9,1) = -2*z*(-2*x - 2*z + 2)/(-z + 1);
+          dshapes(9,2) = -2*z*(-2*x - 2*z + 2)/(-z + 1) - 2*z*(-2*y - 2*z + 2)/(-z + 1) + z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + (-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1);
+          dshapes(10,0) = 2*z*(-2*y - 2*z + 2)/(-z + 1);
+          dshapes(10,1) = -4*x*z/(-z + 1);
+          dshapes(10,2) = -4*x*z/(-z + 1) + 2*x*z*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*x*(-2*y - 2*z + 2)/(-z + 1);
+          dshapes(11,0) = 4*y*z/(-z + 1);
+          dshapes(11,1) = 4*x*z/(-z + 1);
+          dshapes(11,2) = 4*x*y*z/((-z + 1)*(-z + 1)) + 4*x*y/(-z + 1);
+          dshapes(12,0) = -4*y*z/(-z + 1);
+          dshapes(12,1) = 2*z*(-2*x - 2*z + 2)/(-z + 1);
+          dshapes(12,2) = -4*y*z/(-z + 1) + 2*y*z*(-2*x - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*y*(-2*x - 2*z + 2)/(-z + 1);
+          break;
+        }
+
       case HEX:
 	{
           // if (typeid(T) == typeid(SIMD<double>)) return;

libsrc/meshing/meshclass.cpp

@@ -4712,7 +4712,7 @@ namespace netgen
                   lam(2) > -eps && 
                   lam(0) + lam(1) + lam(2) < 1+eps);
       }
-    else if (el.GetType() == PRISM)
+    else if (el.GetType() == PRISM || el.GetType() == PRISM15)
       {
         retval = (lam(0) > -eps &&
                   lam(1) > -eps &&
@@ -4720,7 +4720,7 @@ namespace netgen
                   lam(2) < 1+eps &&
                   lam(0) + lam(1) < 1+eps);
       }
-    else if (el.GetType() == PYRAMID)
+    else if (el.GetType() == PYRAMID || el.GetType() == PYRAMID13)
       {
         retval = (lam(0) > -eps &&
                   lam(1) > -eps &&
@@ -4728,7 +4728,7 @@ namespace netgen
                   lam(0) + lam(2) < 1+eps &&
                   lam(1) + lam(2) < 1+eps);
       }
-    else if (el.GetType() == HEX)
+    else if (el.GetType() == HEX || el.GetType() == HEX20)
       {
         retval = (lam(0) > -eps && lam(0) < 1+eps &&
                   lam(1) > -eps && lam(1) < 1+eps &&

libsrc/meshing/meshtype.cpp

@@ -1028,6 +1028,8 @@ namespace netgen
       case 6: typ = PRISM; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
+      case 13: typ = PYRAMID13; break;
+      case 15: typ = PRISM15; break;
       case 20: typ = HEX20; break;
       default: cerr << "Element::Element: unknown element with " << np << " points" << endl;
       }
@@ -1109,6 +1111,8 @@ namespace netgen
       case 6: typ = PRISM; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
+      case 13: typ = PYRAMID13; break;
+      case 15: typ = PRISM15; break;
       case 20: typ = HEX20; break;
         // 
       default: break;
@@ -1128,7 +1132,9 @@ namespace netgen
       case PRISM: np = 6; break;
       case HEX: np = 8; break;
       case TET10: np = 10; break;
+      case PYRAMID13: np = 13; break;
       case PRISM12: np = 12; break;
+      case PRISM15: np = 15; break;
       case HEX20: np = 20; break;
 
       default: break;
@@ -1958,7 +1964,27 @@ namespace netgen
           shape(4) = p(2);
           break;
         }
-
+      case PYRAMID13:
+        {
+	  T x = p(0);
+	  T y = p(1);
+	  T z = p(2);
+          z *= 1-1e-12;
+          shape[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25);
+          shape[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2));
+          shape[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75);
+          shape[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2));
+          shape[4] = z*(2*z - 1);
+          shape[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          shape[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2);
+          shape[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2);
+          shape[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2);
+          shape[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1);
+          shape[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1);
+          shape[11] = 4*x*y*z/(-z + 1);
+          shape[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1);
+          break;
+        }
       case PRISM:
         {
           shape(0) = p(0) * (1-p(2));
@@ -1969,6 +1995,30 @@ namespace netgen
           shape(5) = (1-p(0)-p(1)) * p(2);
           break;
         }
+      case PRISM15:
+        {
+	  T x = p(0);
+	  T y = p(1);
+	  T z = p(2);
+          T lam = 1-x-y;
+          T lamz = 1-z;
+          shape[0] = (2*x*x-x) * (2*lamz*lamz-lamz);
+          shape[1] = (2*y*y-y) * (2*lamz*lamz-lamz);
+          shape[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz);
+          shape[3] = (2*x*x-x) * (2*z*z-z);
+          shape[4] = (2*y*y-y) * (2*z*z-z);
+          shape[5] = (2*lam*lam-lam) * (2*z*z-z);
+          shape[6] = 4 * x * y * (2*lamz*lamz-lamz);
+          shape[7] = 4 * x * lam * (2*lamz*lamz-lamz);
+          shape[8] = 4 * y * lam * (2*lamz*lamz-lamz);
+          shape[9] = x * 4 * z * (1-z);
+          shape[10] = y * 4 * z * (1-z);
+          shape[11] = lam * 4 * z * (1-z);
+          shape[12] = 4 * x * y * (2*z*z-z);
+          shape[13] = 4 * x * lam * (2*z*z-z);
+          shape[14] = 4 * y * lam * (2*z*z-z);
+          break;
+        }
       case HEX:
         {
           shape(0) = (1-p(0))*(1-p(1))*(1-p(2));
@@ -1981,6 +2031,43 @@ namespace netgen
           shape(7) = (1-p(0))*(  p(1))*(  p(2));
           break;
         }
+      case HEX20:
+	{
+	  T x = p(0);
+	  T y = p(1);
+	  T z = p(2);
+	  shape[0] = (1-x)*(1-y)*(1-z);
+	  shape[1] =    x *(1-y)*(1-z);
+	  shape[2] =    x *   y *(1-z);
+	  shape[3] = (1-x)*   y *(1-z);
+	  shape[4] = (1-x)*(1-y)*(z);
+	  shape[5] =    x *(1-y)*(z);
+	  shape[6] =    x *   y *(z);
+	  shape[7] = (1-x)*   y *(z);
+
+          T sigma[8]={(1-x)+(1-y)+(1-z),x+(1-y)+(1-z),x+y+(1-z),(1-x)+y+(1-z),
+                      (1-x)+(1-y)+z,x+(1-y)+z,x+y+z,(1-x)+y+z};
+
+          static const int e[12][2] =
+            {
+              { 0, 1 }, { 2, 3 }, { 3, 0 }, { 1, 2 },
+              { 4, 5 }, { 6, 7 }, { 7, 4 }, { 5, 6 },
+              { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 },
+            };
+
+          for (int i = 0; i < 12; i++)
+            {
+              T lame = shape[e[i][0]]+shape[e[i][1]];
+              T xi = sigma[e[i][1]]-sigma[e[i][0]];
+              shape[8+i] = (1-xi*xi)*lame;
+            }
+          for (int i = 0; i < 12; i++)
+            {
+              shape[e[i][0]] -= 0.5 * shape[8+i];
+              shape[e[i][1]] -= 0.5 * shape[8+i];
+            }
+          break;
+	}
       default:
         throw NgException("Element :: GetNewShape not implemented for that element");
       }

libsrc/meshing/meshtype.hpp

@@ -21,7 +21,7 @@ namespace netgen
     SEGMENT = 1, SEGMENT3 = 2,
     TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14,
     TET = 20, TET10 = 21, 
-    PYRAMID = 22, PRISM = 23, PRISM12 = 24,
+    PYRAMID = 22, PRISM = 23, PRISM12 = 24, PRISM15 = 27, PYRAMID13 = 28,
     HEX = 25, HEX20 = 26
   };
 
@@ -708,16 +708,18 @@ namespace netgen
     ///
     uint8_t GetNV() const
     {
-      __assume(typ >= TET && typ <= HEX20);        
+      __assume(typ >= TET && typ <= PYRAMID13);
       switch (typ)
 	{
         case TET: 
         case TET10: 
           return 4;
-        case PRISM12: 
+        case PRISM12:
-        case PRISM: 
+        case PRISM15:
+        case PRISM:
 	  return 6; 
 	case PYRAMID:
+        case PYRAMID13:
 	  return 5;
 	case HEX:
 	case HEX20:
@@ -795,8 +797,9 @@ namespace netgen
 	{
 	case TET: 
 	case TET10: return 4;
-	case PYRAMID: return 5;
+	case PYRAMID: case PYRAMID13: return 5;
-	case PRISM: 
+	case PRISM:
+        case PRISM15:
 	case PRISM12: return 5;
         case HEX: case HEX20:
           return 6;

libsrc/meshing/python_mesh.cpp

@@ -256,44 +256,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
   py::class_<Element>(m, "Element3D")
     .def(py::init([](int index, py::list vertices)
                   {
-                    Element * newel = nullptr;
+                    std::map<int, ELEMENT_TYPE> types = {{4, TET},
-                    if (py::len(vertices) == 4)
+                                                         {5, PYRAMID},
-                      {
+                                                         {6, PRISM},
-                        newel = new Element(TET);
+                                                         {8, HEX},
-                        for (int i = 0; i < 4; i++)
+                                                         {10, TET10},
-                          (*newel)[i] = py::extract<PointIndex>(vertices[i])();
+                                                         {13, PYRAMID13},
-                        newel->SetIndex(index);
+                                                         {15, PRISM15},
-                      }
+                                                         {20, HEX20}};
-                    else if (py::len(vertices) == 5)
+                    int np = py::len(vertices);
-                      {
+                    auto newel = new Element(types[np]);
-                        newel = new Element(PYRAMID);
+                    for(int i=0; i<np; i++)
-                        for (int i = 0; i < 5; i++)
+                      (*newel)[i] = py::cast<PointIndex>(vertices[i]);
-                          (*newel)[i] = py::extract<PointIndex>(vertices[i])();
+                    newel->SetIndex(index);
-                        newel->SetIndex(index);
-                      }
-                    else if (py::len(vertices) == 6)
-                      {
-                        newel = new Element(PRISM);
-                        for (int i = 0; i < 6; i++)
-                          (*newel)[i] = py::extract<PointIndex>(vertices[i])();
-                        newel->SetIndex(index);
-                      }
-                    else if (py::len(vertices) == 8)
-                      {
-                        newel = new Element(HEX);
-                        for (int i = 0; i < 8; i++)
-                          (*newel)[i] = py::extract<PointIndex>(vertices[i])();
-                        newel->SetIndex(index);
-                      }
-                    else if (py::len(vertices) == 20)
-                      {
-                        newel = new Element(HEX20);
-                        for (int i = 0; i < 20; i++)
-                          (*newel)[i] = py::extract<PointIndex>(vertices[i])();
-                        newel->SetIndex(index);
-                      }
-                    else
-                      throw NgException ("cannot create element");
                     return newel;
                   }),
           py::arg("index")=1,py::arg("vertices"),

libsrc/meshing/secondorder.cpp

@@ -223,6 +223,29 @@ namespace netgen
 	    { 3, 4, 10 },
 	    { 4, 5, 11 },
 	  };
+        static int betw_prism15[9][3] =
+          {
+           { 0, 1, 6 },
+           { 0, 2, 7 },
+           { 1, 2, 8 },
+           { 0, 3, 9 },
+           { 1, 4, 10 },
+           { 2, 5, 11 },
+           { 3, 4, 12 },
+           { 3, 5, 13 },
+           { 4, 5, 14 }
+          };
+        static int betw_pyramid[8][3] =
+          {
+           { 0, 1, 5 },
+           { 3, 2, 6 },
+           { 3, 0, 7 },
+           { 1, 2, 8 },
+           { 0, 4, 9 },
+           { 1, 4, 10 },
+           { 2, 4, 11 },
+           { 3, 4, 12 }
+          };
 	static int betw_hex[12][3] =
 	  {
             { 0, 1, 8 },
@@ -259,6 +282,21 @@ namespace netgen
 	      onp = 6;
 	      break;
 	    }
+          case PRISM15:
+            {
+              betw = betw_prism15;
+              newel.SetType(PRISM15);
+              onp = 6;
+              break;
+            }
+          case PYRAMID:
+          case PYRAMID13:
+            {
+              betw = betw_pyramid;
+              newel.SetType(PYRAMID13);
+              onp = 5;
+              break;
+            }
           case HEX:
           case HEX20:
 	    {

libsrc/meshing/topology.hpp

@@ -212,10 +212,12 @@ inline short int MeshTopology :: GetNVertices (ELEMENT_TYPE et)
       return 4;
 
     case PYRAMID:
+    case PYRAMID13:
       return 5;
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return 6;
 
     case HEX:
@@ -257,10 +259,15 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
 
     case PYRAMID:
       return 5;
+    case PYRAMID13:
+      return 13;
 
     case PRISM:
-    case PRISM12:
       return 6;
+    case PRISM12:
+      return 12;
+    case PRISM15:
+      return 15;
 
     case HEX:
       return 8;
@@ -277,7 +284,7 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et)
 
 inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
 {
-  __assume(et >= SEGMENT && et <= HEX);
+  __assume(et >= SEGMENT && et <= PYRAMID13);
   switch (et)
     {
     case SEGMENT:
@@ -298,10 +305,12 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
       return 6;
 
     case PYRAMID:
+    case PYRAMID13:
       return 8;
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return 9;
 
     case HEX:
@@ -318,7 +327,7 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et)
 
 inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
 {
-  __assume(et >= SEGMENT && et <= HEX);  
+  __assume(et >= SEGMENT && et <= PYRAMID13);
   switch (et)
     {
     case SEGMENT:
@@ -339,10 +348,12 @@ inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et)
       return 4;
 
     case PYRAMID:
+    case PYRAMID13:
       return 5;
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return 5;
 
     case HEX:
@@ -443,10 +454,12 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et)
       return tet_edges;
 
     case PYRAMID:
+    case PYRAMID13:
       return pyramid_edges;
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return prism_edges;
 
     case HEX:
@@ -542,10 +555,12 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et)
       return tet_edges;
 
     case PYRAMID:
+    case PYRAMID13:
       return pyramid_edges;
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return prism_edges;
 
     case HEX:
@@ -627,9 +642,11 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et)
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return prism_faces;
 
     case PYRAMID:
+    case PYRAMID13:
       return pyramid_faces;
 
     case SEGMENT:
@@ -710,9 +727,11 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et)
 
     case PRISM:
     case PRISM12:
+    case PRISM15:
       return prism_faces;
 
     case PYRAMID:
+    case PYRAMID13:
       return pyramid_faces;
 
     case SEGMENT:

libsrc/visualization/vsmesh.cpp

@@ -2682,7 +2682,7 @@ namespace netgen
     for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++)
       {
 	const Element & el = (*mesh)[ei];
-	if (el.GetType() == PYRAMID && !el.IsDeleted())
+	if ((el.GetType() == PYRAMID || el.GetType() == PYRAMID13) && !el.IsDeleted())
 	  {
             int i = ei + 1;
 

libsrc/visualization/vssolution.cpp

@@ -2857,6 +2857,7 @@ namespace netgen
               }
             case PRISM:
             case PRISM12:
+            case PRISM15:
               {
                 lami[0] = (1-lam3) * (1-lam1-lam2);
                 lami[1] = (1-lam3) * lam1;
@@ -2907,6 +2908,7 @@ namespace netgen
               }
             case PRISM:
             case PRISM12:
+            case PRISM15:
               {
                 lami[0] = (1-lam3) * (1-lam1-lam2);
                 lami[1] = (1-lam3) * lam1;
@@ -2918,6 +2920,7 @@ namespace netgen
                 break;
               }
             case PYRAMID:
+            case PYRAMID13:
               {
                 if (lam3 > 1-1e-5)
                   {
@@ -4019,7 +4022,7 @@ namespace netgen
         if(vispar.donotclipdomain > 0 && vispar.donotclipdomain == (*mesh)[ei].GetIndex()) continue;
 
         ELEMENT_TYPE type = (*mesh)[ei].GetType();
-        if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID)
+        if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID || type == PYRAMID13 || type == PRISM15 || type == HEX20)
           {
             const Element & el = (*mesh)[ei];
 
@@ -4078,6 +4081,8 @@ namespace netgen
                       switch (type)
                         {
                         case PRISM:
+                        case PRISM12:
+                        case PRISM15:
                           if (ix+iy <= n)
                             {
                               ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
@@ -4088,9 +4093,11 @@ namespace netgen
                             compress[ii] = -1;
                           break;
                         case HEX:
+                        case HEX20:
                           ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n);
                           break;
                         case PYRAMID:
+                        case PYRAMID13:
                           ploc = Point<3> (double(ix) / n * (1-double(iz)/n),
                                            double(iy) / n * (1-double(iz)/n),
                                            double(iz)/n);
@@ -4104,7 +4111,7 @@ namespace netgen
                         locgrid[compress[ii]] = ploc;
                     }
 
-            if (type != TET && type != TET10 && type != PRISM) cnt_valid = n3;
+            if (type != TET && type != TET10 && type != PRISM && type != PRISM12 && type != PRISM15) cnt_valid = n3;
 
 	    locgrid.SetSize(cnt_valid);
 

python/read_gmsh.py

@@ -1,9 +1,10 @@
 from netgen.meshing import *
 
-
 def ReadGmsh(filename):
+    if not filename.endswith(".msh"):
+        filename += ".msh"
     meshdim = 1
-    with open(filename + ".msh", 'r') as f:
+    with open(filename, 'r') as f:
         while f.readline().split()[0] != "$Elements":
             pass
         nelem = int(f.readline())
@@ -16,15 +17,58 @@ def ReadGmsh(filename):
                 meshdim = 3
                 break
 
-    f = open(filename + ".msh", 'r')
+    f = open(filename, 'r')
     mesh = Mesh(dim=meshdim)
 
     pointmap = {}
     facedescriptormap = {}
-    namemap = {}
+    namemap = { 0 : "default" }
     materialmap = {}
     bbcmap = {}
 
+    segm = 1
+    trig = 2
+    quad = 3
+    tet = 4
+    hex = 5
+    prism = 6
+    pyramid = 7
+    segm3 = 8      # 2nd order line
+    trig6 = 9      # 2nd order trig
+    tet10 = 11     # 2nd order tet
+    point = 15
+    quad8 = 16     # 2nd order quad
+    hex20 = 17     # 2nd order hex
+    prism15 = 18   # 2nd order prism
+    pyramid13 = 19 # 2nd order pyramid
+    segms = [segm, segm3]
+    trigs = [trig, trig6]
+    quads = [quad, quad8]
+    tets = [tet, tet10]
+    hexes = [hex, hex20]
+    prisms = [prism, prism15]
+    pyramids = [pyramid, pyramid13]
+    elem0d = [point]
+    elem1d = segms
+    elem2d = trigs + quads
+    elem3d = tets + hexes + prisms + pyramids
+
+    num_nodes_map = { segm : 2,
+                      trig : 3,
+                      quad : 4,
+                      tet : 4,
+                      hex : 8,
+                      prism : 6,
+                      pyramid : 5,
+                      segm3 : 3,
+                      trig6 : 6,
+                      tet10 : 10,
+                      point : 1,
+                      quad8 : 8,
+                      hex20 : 20,
+                      prism15 : 18,
+                      pyramid13 : 19 }
+
     while True:
         line = f.readline()
         if line == "":
@@ -57,29 +101,14 @@ def ReadGmsh(filename):
                 # the first tag is the physical group nr, the second tag is the group nr of the dim
                 tags = [int(line[3 + k]) for k in range(numtags)]
 
-                if elmtype == 1:  # 2-node line
+                if elmtype not in num_nodes_map:
-                    num_nodes = 2
-                elif elmtype == 2:  # 3-node trig
-                    num_nodes = 3
-                elif elmtype == 3:  # 4-node quad
-                    num_nodes = 4
-                elif elmtype == 4:  # 4-node tet
-                    num_nodes = 4
-                elif elmtype == 5:  # 8-node hex
-                    num_nodes = 8
-                elif elmtype == 6:  # 6-node prism
-                    num_nodes = 6
-                elif elmtype == 7:  # 5-node pyramid
-                    num_nodes = 5
-                elif elmtype == 15:  # 1-node point
-                    num_nodes = 1
-                else:
                     raise Exception("element type", elmtype, "not implemented")
+                num_nodes = num_nodes_map[elmtype]
 
                 nodenums = line[3 + numtags:3 + numtags + num_nodes]
                 nodenums2 = [pointmap[int(nn)] for nn in nodenums]
 
-                if elmtype == 1:
+                if elmtype in elem1d:
                     if meshdim == 3:
                         if tags[1] in bbcmap:
                             index = bbcmap[tags[1]]
@@ -117,7 +146,7 @@ def ReadGmsh(filename):
 
                     mesh.Add(Element1D(index=index, vertices=nodenums2))
 
-                if elmtype in [2, 3]:  # 2d elements
+                if elmtype in elem2d:  # 2d elements
                     if meshdim == 3:
                         if tags[1] in facedescriptormap.keys():
                             index = facedescriptormap[tags[1]]
@@ -142,9 +171,17 @@ def ReadGmsh(filename):
                                 mesh.SetMaterial(index, "surf" + str(tags[1]))
                             materialmap[tags[1]] = index
 
-                    mesh.Add(Element2D(index, nodenums2))
+                    if elmtype in trigs:
+                        ordering = [i for i in range(3)]
+                        if elmtype == trig6:
+                            ordering += [4,5,3]
+                    if elmtype in quads:
+                        ordering = [i for i in range(4)]
+                        if elmtype == quad8:
+                            ordering += [4, 6, 7, 5]
+                    mesh.Add(Element2D(index, [nodenums2[i] for i in ordering]))
 
-                if elmtype in [4, 5, 6, 7]:  # volume elements
+                if elmtype in elem3d:  # volume elements
                     if tags[1] in materialmap:
                         index = materialmap[tags[1]]
                     else:
@@ -157,9 +194,22 @@ def ReadGmsh(filename):
 
                     nodenums2 = [pointmap[int(nn)] for nn in nodenums]
 
-                    if elmtype == 4:
+                    if elmtype in tets:
-                        mesh.Add(Element3D(index, [nodenums2[0], nodenums2[1], nodenums2[3], nodenums2[2]]))
+                        ordering = [0,1,2,3]
-                    elif elmtype in [5, 6, 7]:
+                        if elmtype == tet10:
-                        mesh.Add(Element3D(index, nodenums2))
+                            ordering += [4,6,7,5,9,8]
+                    elif elmtype in hexes:
+                        ordering = [0,1,5,4,3,2,6,7]
+                        if elmtype == hex20:
+                            ordering += [8,16,10,12,13,19,15,14,9,11,18,17]
+                    elif elmtype in prisms:
+                        ordering = [0,2,1,3,5,4]
+                        if elmtype == prism15:
+                            ordering += [7,6,9,8,11,10,13,12,14]
+                    elif elmtype in pyramids:
+                        ordering = [3,2,1,0,4]
+                        if elmtype == pyramid13:
+                            ordering += [10,5,6,8,12,11,9,7]
+                    mesh.Add(Element3D(index, [nodenums2[i] for i in ordering]))
 
     return mesh