quad8 and hex20 geometric elements

libsrc/meshing/curvedelems.cpp

@@ -2887,6 +2887,46 @@ namespace netgen
 	  break;
         }
         
+      case HEX20:
+	{
+	  shapes = 0.0;
+	  T x = xi(0);
+	  T y = xi(1);
+	  T z = xi(2);
+	  
+	  shapes[0] = (1-x)*(1-y)*(1-z);
+	  shapes[1] =    x *(1-y)*(1-z);
+	  shapes[2] =    x *   y *(1-z);
+	  shapes[3] = (1-x)*   y *(1-z);
+	  shapes[4] = (1-x)*(1-y)*(z);
+	  shapes[5] =    x *(1-y)*(z);
+	  shapes[6] =    x *   y *(z);
+	  shapes[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 = shapes[e[i][0]]+shapes[e[i][1]];
+              T xi = sigma[e[i][1]]-sigma[e[i][0]];
+              shapes[8+i] = (1-xi*xi)*lame;
+            }
+          for (int i = 0; i < 12; i++)
+            {
+              shapes[e[i][0]] -= 0.5 * shapes[8+i];
+              shapes[e[i][1]] -= 0.5 * shapes[8+i];
+            }
+          break;
+	}
+
       default:
 	throw NgException("CurvedElements::CalcShape 3d, element type not handled");
 
@@ -3489,12 +3529,49 @@ namespace netgen
 	   *testout << "quad, num dshape = " << endl << dshapes << endl;
 	   */
 	  break;
+	}
+      case HEX20:
+        {
+          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[20];
           
+          ad[0] = (1-x)*(1-y)*(1-z);
+	  ad[1] =    x *(1-y)*(1-z);
+	  ad[2] =    x *   y *(1-z);
+	  ad[3] = (1-x)*   y *(1-z);
+	  ad[4] = (1-x)*(1-y)*(z);
+	  ad[5] =    x *(1-y)*(z);
+	  ad[6] =    x *   y *(z);
+	  ad[7] = (1-x)*   y *(z);
           
+          AutoDiff<3,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++)
+            {
+              auto lame = ad[e[i][0]]+ad[e[i][1]];
+              auto xi = sigma[e[i][1]]-sigma[e[i][0]];
+              ad[8+i] = (1-xi*xi)*lame;
+            }
+          for (int i = 0; i < 12; i++)
+            {
+              ad[e[i][0]] -= 0.5 * ad[8+i];
+              ad[e[i][1]] -= 0.5 * ad[8+i];
+            }
+          for (int i = 0; i < 20; i++)
+            for (int j = 0; j < 3; j++)
+              dshapes(i,j) = ad[i].DValue(j);
           break;
         }
-
       default:
 	throw NgException("CurvedElements::CalcDShape 3d, element type not handled");
       }

libsrc/meshing/meshtype.cpp

@@ -1028,6 +1028,7 @@ namespace netgen
       case 6: typ = PRISM; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
+      case 20: typ = HEX20; break;
       default: cerr << "Element::Element: unknown element with " << np << " points" << endl;
       }
     orderx = ordery = orderz = 1;
@@ -1108,6 +1109,7 @@ namespace netgen
       case 6: typ = PRISM; break;
       case 8: typ = HEX; break;
       case 10: typ = TET10; break;
+      case 20: typ = HEX20; break;
         // 
       default: break;
         cerr << "Element::SetNP unknown element with " << np << " points" << endl;
@@ -1127,6 +1129,7 @@ namespace netgen
       case HEX: np = 8; break;
       case TET10: np = 10; break;
       case PRISM12: np = 12; break;
+      case HEX20: np = 20; break;
 
       default: break;
         cerr << "Element::SetType unknown type  " << int(typ) << endl;

libsrc/meshing/meshtype.hpp

@@ -647,7 +647,7 @@ namespace netgen
     ///
     ELEMENT_TYPE typ;
     /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism)
-    int np:5;
+    int np:6;
     ///
     class flagstruct { 
     public:
@@ -708,7 +708,7 @@ namespace netgen
     ///
     uint8_t GetNV() const
     {
-      __assume(typ >= TET && typ <= HEX);        
+      __assume(typ >= TET && typ <= HEX20);        
       switch (typ)
 	{
         case TET: 
@@ -720,6 +720,7 @@ namespace netgen
 	case PYRAMID:
 	  return 5;
 	case HEX:
+	case HEX20:
 	  return 8;
         default: // not a 3D element
 #ifdef DEBUG
@@ -797,7 +798,8 @@ namespace netgen
 	case PYRAMID: return 5;
 	case PRISM: 
 	case PRISM12: return 5;
-        case HEX: return 6;
+        case HEX: case HEX20:
+          return 6;
 	default:
 #ifdef DEBUG
 	  PrintSysError ("element3d::GetNFaces not implemented for typ", typ)

libsrc/meshing/secondorder.cpp

@@ -223,6 +223,22 @@ namespace netgen
 	    { 3, 4, 10 },
 	    { 4, 5, 11 },
 	  };
+	static int betw_hex[12][3] =
+	  {
+            { 0, 1, 8 },
+            { 2, 3, 9 },
+            { 3, 0, 10 },
+            { 1, 2, 11 },
+            { 4, 5, 12 },
+            { 6, 7, 13 },
+            { 7, 4, 14 },
+            { 5, 6, 15 },
+            { 0, 4, 16 },
+            { 1, 5, 17 },
+            { 2, 6, 18 },
+            { 3, 7, 19 },
+          };
+        
 	int (*betw)[3] = NULL;
 
 	switch (el.GetType())
@@ -243,6 +259,14 @@ namespace netgen
 	      onp = 6;
 	      break;
 	    }
+          case HEX:
+          case HEX20:
+	    {
+	      betw = betw_hex;
+	      newel.SetType (HEX20);
+	      onp = 8;
+	      break;
+	    }
 	  default:
 	    PrintSysError ("MakeSecondOrder, illegal vol type ", int(el.GetType()));
 	  }

libsrc/meshing/topology.cpp

@@ -543,7 +543,7 @@ namespace netgen
         cnt = 0;
 
         ParallelForRange
-          (tm, mesh->Points().Size(),
+          (tm, mesh->GetNV(), // Points().Size(),
            [&] (size_t begin, size_t end)
            {
              INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10); 
@@ -581,7 +581,7 @@ namespace netgen
 
         // accumulate number of edges
         int ned = edge2vert.Size();
-        for (auto v : mesh->Points().Range())
+        for (size_t v = 0; v < mesh->GetNV(); v++)
           {
             auto hv = cnt[v];
             cnt[v] = ned;
@@ -595,7 +595,7 @@ namespace netgen
 	// for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++)
 
         ParallelForRange
-          (tm, mesh->Points().Size(),
+          (tm, mesh->GetNV(), // Points().Size(),
            [&] (size_t begin, size_t end)
            {
              INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10);
@@ -719,7 +719,7 @@ namespace netgen
         // for (auto v : mesh.Points().Range())
         NgProfiler::StartTimer (timer2b1);
         ParallelForRange
-          (tm, mesh->Points().Size(),
+          (tm, mesh->GetNV(), // Points().Size(),
             [&] (size_t begin, size_t end)
             {
               INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10); 
@@ -754,7 +754,8 @@ namespace netgen
         
         // accumulate number of faces
         int nfa = oldnfa;
-        for (auto v : mesh->Points().Range())
+        // for (auto v : Range(mesh->GetNV())) // Points().Range())
+        for (size_t v = 0; v < mesh->GetNV(); v++)          
           {
             auto hv = cnt[v];
             cnt[v] = nfa;
@@ -765,7 +766,7 @@ namespace netgen
         // for (auto v : mesh.Points().Range())
 
         ParallelForRange
-          (tm, mesh->Points().Size(),
+          (tm, mesh->GetNV(), // Points().Size(),
             [&] (size_t begin, size_t end)
             {
               INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10);