From 82e88f3afb5570c2d1a72642cad9e563ea10d0d5 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 12 Jul 2023 10:50:03 -0700 Subject: [PATCH 1/3] Alefeld split hack --- libsrc/meshing/hpref_trig.hpp | 44 +++++++++++++++++++++++++++++++++ libsrc/meshing/hprefinement.cpp | 36 ++++++++++++++++----------- libsrc/meshing/hprefinement.hpp | 3 ++- 3 files changed, 67 insertions(+), 16 deletions(-) diff --git a/libsrc/meshing/hpref_trig.hpp b/libsrc/meshing/hpref_trig.hpp index 3676ad0f..890d79ff 100644 --- a/libsrc/meshing/hpref_trig.hpp +++ b/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 +}; + diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 8c673efa..a3b38f1e 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/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: @@ -1338,7 +1341,7 @@ namespace netgen sing = true; // iterate at least once while(sing) { - PrintMessage(3, " Start new hp-refinement: step ", act_ref); + PrintMessage(0, " Start new hp-refinement: step ", act_ref); DoRefinement (mesh, hpelements, ref, fac1); DoRefineDummies (mesh, hpelements, ref); @@ -1823,6 +1826,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS bool ClassifyHPElements (Mesh & mesh, NgArray & elements, int & act_ref, int & levels) { + cout << "in classify" << endl; INDEX_2_HASHTABLE edges(mesh.GetNSeg()+1); NgBitArray edgepoint(mesh.GetNP()); INDEX_2_HASHTABLE edgepoint_dom(mesh.GetNSeg()+1); @@ -1841,7 +1845,9 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & 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) + sing = true; // Joachim if(sing==0) return(sing); int cnt_undef = 0, cnt_nonimplement = 0; @@ -1859,12 +1865,11 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & 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 & edges, INDEX_2_HAS break; } case HP_TRIG: - { + { int dim = mesh.GetDimension(); const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); - + hpel.type = HP_TRIG_ALEFELD; + + /* hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint, dim, fd); - - dd = 2; + */ + dd = 2; + continue; 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,6 +1935,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS throw NgException ("hprefinement.cpp: don't know how to set parameters"); } } + + continue; if(hpel.type == HP_NONE) cnt_undef++; @@ -1935,8 +1944,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS //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 +1966,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS } } - - + PrintMessage(3, "undefined elements update classification: ", cnt_undef); PrintMessage(3, "non-implemented in update classification: ", cnt_nonimplement); diff --git a/libsrc/meshing/hprefinement.hpp b/libsrc/meshing/hprefinement.hpp index 0affd82b..c03a2ba9 100644 --- a/libsrc/meshing/hprefinement.hpp +++ b/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,7 +313,6 @@ public: }; - DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, int levels, double fac1=0.125, bool setorders=true, bool ref_level = false); From 5b19ea6451777e75abc24ddfc2307660b2ec75ae Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 12 Jul 2023 17:26:32 -0700 Subject: [PATCH 2/3] enum for macro-based element splitting --- libsrc/include/nginterface_v2.hpp | 2 ++ libsrc/interface/nginterface.cpp | 2 +- libsrc/interface/nginterface_v2.cpp | 10 ++++++++- libsrc/meshing/hprefinement.cpp | 33 +++++++++++++++-------------- libsrc/meshing/hprefinement.hpp | 9 +++++++- 5 files changed, 37 insertions(+), 19 deletions(-) diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index 6c70a9ef..d06c7de9 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/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; diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index 4da0eed9..c76a94d4 100644 --- a/libsrc/interface/nginterface.cpp +++ b/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 (mesh->GetGeometry()->GetRefinement()); - HPRefinement (*mesh, &ref, levels, parameter, setorders, ref_level); + HPRefinement (*mesh, &ref, SPLIT_HP, levels, parameter, setorders, ref_level); /* Refinement * ref; diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index a953dce0..c2f295bb 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1226,8 +1226,16 @@ namespace netgen { NgLock meshlock (mesh->MajorMutex(), true); Refinement & ref = const_cast (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 (mesh->GetGeometry()->GetRefinement()); + ::netgen::HPRefinement (*mesh, &ref, SPLIT_ALEFELD, 1, 0.5, true, true); + } + int Ngx_Mesh::GetElementOrder (int enr) const { diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index a3b38f1e..583ed5cb 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -557,7 +557,7 @@ namespace netgen NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE & faces, INDEX_2_HASHTABLE & face_edges, INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int & levels, int & act_ref); - bool ClassifyHPElements (Mesh & mesh, NgArray & elements, int & act_ref, int & levels); + bool ClassifyHPElements (Mesh & mesh, NgArray & elements, SplittingType split, int & act_ref, int & levels); void InitHPElements(Mesh & mesh, NgArray & elements) @@ -1309,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); @@ -1336,12 +1337,12 @@ 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) { - PrintMessage(0, " Start new hp-refinement: step ", act_ref); + PrintMessage(3, " Start new hp-refinement: step ", act_ref); DoRefinement (mesh, hpelements, ref, fac1); DoRefineDummies (mesh, hpelements, ref); @@ -1445,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."); @@ -1824,9 +1825,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS - bool ClassifyHPElements (Mesh & mesh, NgArray & elements, int & act_ref, int & levels) + bool ClassifyHPElements (Mesh & mesh, NgArray & elements, SplittingType split, int & act_ref, int & levels) { - cout << "in classify" << endl; INDEX_2_HASHTABLE edges(mesh.GetNSeg()+1); NgBitArray edgepoint(mesh.GetNP()); INDEX_2_HASHTABLE edgepoint_dom(mesh.GetNSeg()+1); @@ -1846,8 +1846,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint, levels, act_ref); - if (act_ref == 1) - sing = true; // Joachim + if (act_ref == 1 && split == SPLIT_ALEFELD) + sing = true; if(sing==0) return(sing); int cnt_undef = 0, cnt_nonimplement = 0; @@ -1894,13 +1894,14 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS case HP_TRIG: { int dim = mesh.GetDimension(); - const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); - hpel.type = HP_TRIG_ALEFELD; - - /* - hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, - faces, face_edges, surf_edges, facepoint, dim, fd); - */ + const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); + + if (split == SPLIT_HP) + hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, + faces, face_edges, surf_edges, facepoint, dim, fd); + else + hpel.type = HP_TRIG_ALEFELD; + dd = 2; continue; break; diff --git a/libsrc/meshing/hprefinement.hpp b/libsrc/meshing/hprefinement.hpp index c03a2ba9..3095a38b 100644 --- a/libsrc/meshing/hprefinement.hpp +++ b/libsrc/meshing/hprefinement.hpp @@ -313,9 +313,16 @@ public: }; -DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, int levels, +enum SplittingType { SPLIT_HP, SPLIT_ALEFELD }; + +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 From 06070d49f3002dd2566a247bb194a703693999cf Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Wed, 12 Jul 2023 17:31:14 -0700 Subject: [PATCH 3/3] little cleanup --- libsrc/meshing/hprefinement.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 583ed5cb..22db74d8 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -1903,7 +1903,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS hpel.type = HP_TRIG_ALEFELD; dd = 2; - continue; break; } case HP_QUAD: @@ -1937,8 +1936,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS } } - continue; - if(hpel.type == HP_NONE) cnt_undef++;