Alfeld split in 3D

Author: JSchoeberl

libsrc/interface/nginterface_v2.cpp

@@ -1178,7 +1178,7 @@ namespace netgen
   {
     NgLock meshlock (mesh->MajorMutex(), true);
     Refinement & ref = const_cast<Refinement&> (mesh->GetGeometry()->GetRefinement());
-    ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALFELD, 1, 1.0/3.0, true, true);
+    ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALFELD, 1, 1.0 / (mesh->GetDimension()+1), true, true);
   }
 
 

libsrc/meshing/hpref_tet.hpp

@@ -4184,6 +4184,47 @@ HPRef_Struct reftet_3f_0e_0v =
 
 
 
+
+// HP_TET_ALFELD
+int reftet_Alfeld_splitedges[][3] =
+{
+  { 0, 0, 0 }
+};
+int reftet_Alfeld_splitfaces[][4] =
+{
+  { 0, 0, 0, 0 }
+};
+
+int reftet_Alfeld_splitelements[][5] =
+  {
+    { 1, 2, 3, 4, 5},
+    { 0 }
+  };
+
+HPREF_ELEMENT_TYPE reftet_Alfeld_newelstypes[] =
+{
+  HP_TET, HP_TET, HP_TET, HP_TET,
+  HP_NONE, 
+};
+
+int reftet_Alfeld_newels[][8] =
+{
+  { 1, 2, 3, 5 }, 
+  { 1, 4, 2, 5 }, 
+  { 1, 3, 4, 5 },
+  { 3, 2, 4, 5 }
+};
+HPRef_Struct reftet_Alfeld =
+{
+  HP_TET,
+  reftet_Alfeld_splitedges, 
+  reftet_Alfeld_splitfaces,
+  reftet_Alfeld_splitelements,  
+  reftet_Alfeld_newelstypes,
+  reftet_Alfeld_newels
+};
+
+
 /*
 
 */

libsrc/meshing/hprefinement.cpp

@@ -429,6 +429,10 @@ namespace netgen
         hps = &reftet_2f_1e_0vb; break;
       case HP_TET_3F_0E_0V:
         hps = &reftet_3f_0e_0v; break;
+
+      case HP_TET_ALFELD:
+        hps = &reftet_Alfeld; break;
+
 
       case HP_PRISM:
 	hps = &refprism; break;
@@ -1952,6 +1956,9 @@ namespace netgen
               else
                 cout << "type = " << hpel.type << endl;
               */
+              if (split == SPLIT_ALFELD && mesh.GetDimension()==3)
+                hpel.type = HP_TET_ALFELD;
+
 	      break;
 	    }
 	  case HP_PRISM:
@@ -1982,9 +1989,9 @@ namespace netgen
               if (split == SPLIT_HP)
                 hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, 
                                          faces, face_edges, surf_edges, facepoint, dim, fd);    
-              else if (split == SPLIT_ALFELD)
+              else if (split == SPLIT_ALFELD && mesh.GetDimension()==2)
                 hpel.type = HP_TRIG_ALFELD;
-              else if (split == SPLIT_POWELL)
+              else if (split == SPLIT_POWELL && mesh.GetDimension()==2)
                 hpel.type = HP_TRIG_POWELL;
 
 	      dd = 2;

libsrc/meshing/hprefinement.hpp

@@ -189,6 +189,8 @@ enum HPREF_ELEMENT_TYPE {
 
   HP_TET_3F_0E_0V = 700,  // 3 singular faces, no additional points or edges  NEW   done 
 
+  HP_TET_ALFELD,
+  
   HP_PRISM = 1000,
   HP_PRISM_SINGEDGE,
   HP_PRISM_SINGEDGE_V12,