Merge branch 'splitting' into 'master'

Splitting

See merge request ngsolve/netgen!583

libsrc/include/nginterface_v2.hpp

@@ -389,6 +389,8 @@ namespace netgen
     // also added from nginterface.h, still 1-based, need redesign
     void HPRefinement (int levels, double parameter = 0.125,
                        bool setorders = true,bool ref_level = false);
+    void SplitAlefeld ();
+    
     size_t GetNP() const;
     int GetSurfaceElementSurfaceNumber (size_t ei) const;
     int GetSurfaceElementFDNumber (size_t ei) const;

libsrc/interface/nginterface.cpp

@@ -1122,7 +1122,7 @@ void Ng_HPRefinement (int levels, double parameter, bool setorders,
 {
   NgLock meshlock (mesh->MajorMutex(), true);
   Refinement & ref = const_cast<Refinement&> (mesh->GetGeometry()->GetRefinement());
-  HPRefinement (*mesh, &ref, levels, parameter, setorders, ref_level);
+  HPRefinement (*mesh, &ref, SPLIT_HP, levels, parameter, setorders, ref_level);
   /*
     Refinement * ref;
 

libsrc/interface/nginterface_v2.cpp

@@ -1226,8 +1226,16 @@ namespace netgen
   {
     NgLock meshlock (mesh->MajorMutex(), true);
     Refinement & ref = const_cast<Refinement&> (mesh->GetGeometry()->GetRefinement());
-    ::netgen::HPRefinement (*mesh, &ref, levels, parameter, setorders, ref_level);
+    ::netgen::HPRefinement (*mesh, &ref, SPLIT_HP, levels, parameter, setorders, ref_level);
   }
+
+  void Ngx_Mesh::SplitAlefeld ()
+  {
+    NgLock meshlock (mesh->MajorMutex(), true);
+    Refinement & ref = const_cast<Refinement&> (mesh->GetGeometry()->GetRefinement());
+    ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALEFELD, 1, 0.5, true, true);
+  }
+
   
 int Ngx_Mesh::GetElementOrder (int enr) const
 {

libsrc/meshing/hpref_trig.hpp

@@ -774,3 +774,47 @@ HPRef_Struct reftrig_3singedges =
   reftrig_3singedges_newelstypes, 
   reftrig_3singedges_newels
 };
+
+
+
+
+
+
+
+
+
+
+
+
+// HP_TRIG_3SINGEDGES
+int reftrig_Alefeld_splitedges[][3] =
+{
+  { 0, 0, 0 }
+};
+int reftrig_Alefeld_splitfaces[][4] =
+{
+  { 1, 2, 3, 4 },
+  { 0, 0, 0, 0 }
+};
+
+HPREF_ELEMENT_TYPE reftrig_Alefeld_newelstypes[] =
+{
+  HP_TRIG, HP_TRIG, HP_TRIG,
+  HP_NONE, 
+};
+int reftrig_Alefeld_newels[][8] =
+{
+  { 1, 2, 4 }, 
+  { 2, 3, 4 }, 
+  { 3, 1, 4 },
+};
+HPRef_Struct reftrig_Alefeld =
+{
+  HP_TRIG,
+  reftrig_Alefeld_splitedges, 
+  reftrig_Alefeld_splitfaces,
+  0,
+  reftrig_Alefeld_newelstypes,
+  reftrig_Alefeld_newels
+};
+

libsrc/meshing/hprefinement.cpp

@@ -174,7 +174,10 @@ namespace netgen
       case HP_TRIG_3SINGEDGES:
 	hps = &reftrig_3singedges; break;
  
- 
+      case HP_TRIG_ALEFELD:
+        hps = &reftrig_Alefeld; break;
+
+        
       case HP_QUAD:
 	hps = &refquad; break;
       case HP_DUMMY_QUAD_SINGCORNER:
@@ -554,7 +557,7 @@ namespace netgen
 		       NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
 			INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref); 
 
-  bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels);
+  bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, SplittingType split, int & act_ref, int & levels);
   
   
   void  InitHPElements(Mesh & mesh, NgArray<HPRefElement> & elements) 
@@ -1306,7 +1309,8 @@ namespace netgen
 
   /* ***************************** HPRefinement ********************************** */
 
-  void HPRefinement (Mesh & mesh, Refinement * ref, int levels, double fac1, bool setorders, bool reflevels)
+  void HPRefinement (Mesh & mesh, Refinement * ref, SplittingType split,
+                     int levels, double fac1, bool setorders, bool reflevels)
   {
     PrintMessage (1, "HP Refinement called, levels = ", levels);
 
@@ -1333,7 +1337,7 @@ namespace netgen
     nplevel.Append (mesh.GetNP());
     
     int act_ref=1;
-    bool sing = ClassifyHPElements (mesh,hpelements, act_ref, levels); 
+    bool sing = ClassifyHPElements (mesh,hpelements, split, act_ref, levels); 
 
     sing = true; // iterate at least once
     while(sing) 
@@ -1442,7 +1446,7 @@ namespace netgen
 
 	act_ref++; 
 	
-	sing = ClassifyHPElements(mesh,hpelements, act_ref, levels); 
+	sing = ClassifyHPElements(mesh,hpelements, split, act_ref, levels); 
       }
 
     PrintMessage(3, " HP-Refinement done with ", --act_ref, " refinement steps.");
@@ -1821,7 +1825,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 
 
 
-  bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels)
+  bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, SplittingType split, int & act_ref, int & levels)
   {
     INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1);
     NgBitArray edgepoint(mesh.GetNP());
@@ -1841,7 +1845,9 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
     bool sing = CheckSingularities(mesh, edges, edgepoint_dom, 
 			      cornerpoint, edgepoint, faces, face_edges, 
 			      surf_edges, facepoint, levels, act_ref); 
-  
+
+    if (act_ref == 1 && split == SPLIT_ALEFELD)
+      sing = true;   
     if(sing==0) return(sing); 
 
     int cnt_undef = 0, cnt_nonimplement = 0;
@@ -1859,12 +1865,11 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	HPRefElement old_el = elements[i]; 
 	int dd=3; 
 
-
 	if(act_ref !=1 && (hpel.type == HP_HEX || hpel.type == HP_PRISM || hpel.type == HP_TET 
 			   || hpel.type == HP_PYRAMID || hpel.type == HP_QUAD || hpel.type == HP_TRIG || hpel.type == HP_SEGM)) 
 	  continue; 
 	
-	sing = 1; 
+	sing = 1;
 	switch (hprs->geom)
 	  {
 	  case HP_TET:
@@ -1887,23 +1892,25 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	      break; 
 	    } 
 	  case HP_TRIG: 
-	    { 
+	    {
 	      int dim = mesh.GetDimension(); 
-	      const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());	
+	      const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());
-	      
+
-	      hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, 
+              if (split == SPLIT_HP)
-				       faces, face_edges, surf_edges, facepoint, dim, fd);    
+                hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, 
-	     
+                                         faces, face_edges, surf_edges, facepoint, dim, fd);    
-	      dd = 2; 
+              else
+                hpel.type = HP_TRIG_ALEFELD;
+
+	      dd = 2;
 	      break; 
 	    } 
 	  case HP_QUAD: 
 	    { 
 	      int dim = mesh.GetDimension(); 
-	      const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());	
+	      const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());
 	      hpel.type = ClassifyQuad(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, 
 				  faces, face_edges, surf_edges, facepoint, dim, fd);    
-
 	      dd = 2; 
 	      break; 
 	    }
@@ -1928,15 +1935,13 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	      throw NgException ("hprefinement.cpp: don't know how to set parameters");
 	    }
 	  }
-	    
+
 	if(hpel.type == HP_NONE) 
 	  cnt_undef++; 
 
 	//else 
 	//cout << "elem " << i << " classified type " << hpel.type << endl; 
 
-	
-	
 	if (!Get_HPRef_Struct (hpel.type)) 
 	  {
 	    (*testout) << "hp-element-type " << hpel.type << " not implemented   " << endl;
@@ -1959,8 +1964,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
 	  } 
 
       }
-    
+
-    
     PrintMessage(3, "undefined elements update classification: ", cnt_undef);
     PrintMessage(3, "non-implemented in update classification: ", cnt_nonimplement);
 

libsrc/meshing/hprefinement.hpp

@@ -41,6 +41,8 @@ enum HPREF_ELEMENT_TYPE {
   HP_TRIG_SINGEDGES23,
   HP_TRIG_3SINGEDGES = 40,
 
+  HP_TRIG_ALEFELD,
+  
   HP_QUAD = 50,
   HP_QUAD_SINGCORNER,
   HP_DUMMY_QUAD_SINGCORNER,
@@ -311,10 +313,16 @@ public:
 };
 
 
+enum SplittingType { SPLIT_HP, SPLIT_ALEFELD };
 
-DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, int levels,
+DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, SplittingType split, int levels,
 			  double fac1=0.125, bool setorders=true, bool ref_level = false);
 
+inline void HPRefinement (Mesh & mesh, Refinement * ref, int levels,
+			  double fac1=0.125, bool setorders=true, bool ref_level = false)
+{
+  HPRefinement (mesh, ref, SPLIT_HP, levels, fac1, setorders, ref_level);
+}
 
 #endif