hex-filling of thin domains (first steps ...)

libsrc/csg/python_csg.cpp

@@ -502,8 +502,15 @@ DLL_HEADER void ExportCSG()
              SetGlobalMesh (dummy);
              dummy->SetGeometry(geo);
 	     ng_geometry = geo;
-             geo->FindIdenticSurfaces(1e-8 * geo->MaxSize()); 
-             geo->GenerateMesh (dummy, param, 0, 6);
+             geo->FindIdenticSurfaces(1e-8 * geo->MaxSize());
+             try
+               {
+                 geo->GenerateMesh (dummy, param, 0, 6);
+               }
+             catch (NgException ex)
+               {
+                 cout << "Caught NgException: " << ex.What() << endl;
+               }
              return dummy;
            }))
     ;

libsrc/meshing/CMakeLists.txt

@@ -20,7 +20,7 @@ add_library(mesh ${NG_LIB_TYPE}
         smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp		    
         topology.cpp triarls.cpp validate.cpp bcfunctions.cpp   
         parallelmesh.cpp  paralleltop.cpp  paralleltop.hpp basegeom.cpp     
-        python_mesh.cpp
+        python_mesh.cpp hexarls.cpp
         ${mesh_object_libs}
         )
 

libsrc/meshing/curvedelems.cpp

@@ -2590,6 +2590,40 @@ namespace netgen
 	  shapes[5] =    x *(1-y)*(z);
 	  shapes[6] =    x *   y *(z);
 	  shapes[7] = (1-x)*   y *(z);
+
+	  if (info.order == 1) return;
+	  
+	  double mu[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),
+          };
+	    
+	  int ii = 8;
+	  const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (HEX);
+	  
+	  for (int i = 0; i < 8; i++)
+	    {
+	      int eorder = edgeorder[info.edgenrs[i]];
+	      if (eorder >= 2)
+		{
+		  int vi1 = edges[i][0]-1, vi2 = edges[i][1]-1;
+		  if (el[vi1] > el[vi2]) swap (vi1, vi2);
+
+		  CalcEdgeShape (eorder, mu[vi1]-mu[vi2], &shapes(ii));
+		  double lame = shapes(vi1)+shapes(vi2);
+		  for (int j = 0; j < order-1; j++)
+		    shapes(ii+j) *= lame;
+		  ii += eorder-1;
+		}
+	    }
+
+          
 	  break;
 	}
 

libsrc/meshing/meshfunc.cpp

@@ -8,6 +8,7 @@ namespace netgen
   extern const char * prismrules2[];
   extern const char * pyramidrules[];
   extern const char * pyramidrules2[];
+  extern const char * hexrules[];
 
 
   // extern double teterrpow; 
@@ -94,8 +95,9 @@ namespace netgen
 	   }
 	 domain_bbox.Increase (0.01 * domain_bbox.Diam());
 	 
-	
-        for (int qstep = 1; qstep <= 3; qstep++)
+
+         for (int qstep = 1; qstep <= 3; qstep++)
+         // for (int qstep = 0; qstep <= 3; qstep++)  // for hex-filling
 	  {
 	    // cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
 	    if (mesh3d.HasOpenQuads())
@@ -105,6 +107,10 @@ namespace netgen
 		const char ** rulep = NULL;
 		switch (qstep)
 		  {
+		  case 0:
+		    // rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls";
+		    rulep = hexrules;
+		    break;
 		  case 1:
 		    rulefile += "/rules/prisms2.rls";
 		    rulep = prismrules2;

libsrc/meshing/meshtype.cpp

@@ -1165,7 +1165,7 @@ namespace netgen
         { 3, 2, 3, 5, 0 },
         { 3, 3, 4, 5, 0 },
         { 3, 4, 1, 5, 0 } };
-
+    
     static const int prismfaces[][5] =
       {
         { 3, 1, 3, 2, 0 },
@@ -1174,6 +1174,17 @@ namespace netgen
         { 4, 2, 3, 6, 5 },
         { 4, 3, 1, 4, 6 }
       };
+    
+    static const int hexfaces[][5] =
+      {
+        { 4, 4, 3, 2, 1 },
+        { 4, 3, 7, 6, 2 },
+        { 4, 7, 8, 5, 6 },
+        { 4, 8, 4, 1, 5 },
+        { 4, 1, 2, 6, 5 },
+        { 4, 3, 4, 8, 7 }
+      };
+    
 
     switch (np)
       {
@@ -1209,6 +1220,13 @@ namespace netgen
             face.PNum(j) = PNum(prismfaces[i-1][j]);
           break;
         }
+      case 8:
+        {
+          face.SetType(QUAD);
+          for (int j = 1; j <= 4; j++)
+            face.PNum(j) = PNum(hexfaces[i-1][j]);
+          break;
+        }
       }
   }
 

libsrc/meshing/meshtype.hpp

@@ -679,6 +679,7 @@ namespace netgen
 	case PYRAMID: return 5;
 	case PRISM: 
 	case PRISM12: return 5;
+        case HEX: return 6;
 	default:
 #ifdef DEBUG
 	  PrintSysError ("element3d::GetNFaces not implemented for typ", typ)

libsrc/meshing/parser3.cpp

@@ -530,6 +530,21 @@ void vnetrule :: LoadRule (istream & ist)
 		  ist >> elements.Last().PNum(6);
 		  ist >> ch;    // ','
 		}
+              
+	      if (ch == COMMASIGN)
+		{
+		  //		  elements.Last().SetNP(6);
+		  elements.Last().SetType(HEX);
+		  ist >> elements.Last().PNum(7);
+		  ist >> ch;    // ','
+		}
+	      if (ch == COMMASIGN)
+		{
+		  //		  elements.Last().SetNP(6);
+		  elements.Last().SetType(HEX);
+		  ist >> elements.Last().PNum(8);
+		  ist >> ch;    // ','
+		}
 
 	      /*
 	      orientations.Append (fourint());