diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 0cf3bcf1..7be9252b 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -2194,7 +2194,7 @@ void MeshOptimize3d :: SwapImproveSurface ( double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face, Table & elementsonnode, - DynamicTable & belementsonnode, bool check_only ) + DynamicTable & belementsonnode, bool conform_segments, bool check_only ) { PointIndex pi1, pi2, pi3, pi4, pi5; Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); @@ -2228,28 +2228,31 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face, } - bool bface = 0; - for (int k = 0; k < belementsonnode[pi1].Size(); k++) - { - const Element2d & bel = - mesh[belementsonnode[pi1][k]]; + if(!conform_segments) + { + bool bface = 0; + for (int k = 0; k < belementsonnode[pi1].Size(); k++) + { + const Element2d & bel = + mesh[belementsonnode[pi1][k]]; - bool bface1 = 1; - for (int l = 0; l < 3; l++) - if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3) + bool bface1 = 1; + for (int l = 0; l < 3; l++) + if (bel[l] != pi1 && bel[l] != pi2 && bel[l] != pi3) + { + bface1 = 0; + break; + } + + if (bface1) { - bface1 = 0; + bface = 1; break; } - - if (bface1) - { - bface = 1; - break; } - } - if (bface) return 0.0; + if (bface) return 0.0; + } FlatArray row = elementsonnode[pi1]; @@ -2367,11 +2370,11 @@ double MeshOptimize3d :: SwapImprove2 ( ElementIndex eli1, int face, 2 -> 3 conversion */ -void MeshOptimize3d :: SwapImprove2 () +void MeshOptimize3d :: SwapImprove2 (bool conform_segments) { static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t); - if (goal == OPT_CONFORM) return; + if (!conform_segments && goal == OPT_CONFORM) return; mesh.BuildBoundaryEdges(false); @@ -2433,7 +2436,7 @@ void MeshOptimize3d :: SwapImprove2 () for (int j = 0; j < 4; j++) { - double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, true); + double d_badness = SwapImprove2( eli1, j, elementsonnode, belementsonnode, conform_segments, true); if(d_badness<0.0) my_faces_with_improvement.Append( std::make_tuple(d_badness, eli1, j) ); } @@ -2449,7 +2452,7 @@ void MeshOptimize3d :: SwapImprove2 () { if(mesh[eli].IsDeleted()) continue; - if(SwapImprove2( eli, j, elementsonnode, belementsonnode, false) < 0.0) + if(SwapImprove2( eli, j, elementsonnode, belementsonnode, conform_segments, false) < 0.0) cnt++; } diff --git a/libsrc/meshing/improve3.hpp b/libsrc/meshing/improve3.hpp index 5bd0f1aa..92987dd6 100644 --- a/libsrc/meshing/improve3.hpp +++ b/libsrc/meshing/improve3.hpp @@ -45,8 +45,8 @@ public: void SwapImprove (const TBitArray * working_elements = NULL); void SwapImproveSurface (const TBitArray * working_elements = NULL, const NgArray< idmap_type* > * idmaps = NULL); - void SwapImprove2 (); - double SwapImprove2 (ElementIndex eli1, int face, Table & elementsonnode, DynamicTable & belementsonnode, bool check_only=false ); + void SwapImprove2 (bool conform_segments = false); + double SwapImprove2 (ElementIndex eli1, int face, Table & elementsonnode, DynamicTable & belementsonnode, bool conform_segments, bool check_only=false ); void ImproveMesh() { mesh.ImproveMesh(mp, goal); } diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index ce6a4699..b0f6d6db 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -783,7 +783,7 @@ namespace netgen free_segs.Append(segi); auto get_nonconforming = [&] (const auto & p2el) { - Array nonconforming; + Array nonconforming; for (auto segi : free_segs) { auto seg = mesh[segi]; @@ -805,16 +805,6 @@ namespace netgen auto split_segment = [&] (SegmentIndex segi, const auto & p2el) { auto seg = mesh[segi]; auto p_new = Center(mesh[seg[0]], mesh[seg[1]]); - auto pi_new = mesh.AddPoint(p_new); - auto seg_new0 = seg; - auto seg_new1 = seg; - seg_new0[1] = pi_new; - seg_new1[0] = pi_new; - - mesh[segi][0] = PointIndex::INVALID; - mesh.AddSegment(seg_new0); - mesh.AddSegment(seg_new1); - double lam[3]; ElementIndex ei_start = mesh.GetElementOfPoint(p_new, lam, false, domain); @@ -824,6 +814,9 @@ namespace netgen } ei_start -= 1; + if(mesh[ei_start].IsDeleted()) + return; + double max_inside = -1.; ElementIndex ei_max_inside = -1; @@ -832,6 +825,9 @@ namespace netgen for(auto pi : mesh[ei_start].PNums()) { for(auto ei1 : p2el[pi]) { double lam[3]; + + if(mesh[ei1].IsDeleted()) + return; if(!mesh.PointContainedIn3DElement(p_new, lam, ei1+1)) continue; @@ -843,18 +839,34 @@ namespace netgen } } + if(max_inside < 1e-4) { + PrintMessage(3, "Could not find volume element with new point inside"); + return; + } + // split tet into 4 new tests, with new point inside auto el = mesh[ei_max_inside]; if(el.GetNP() != 4) { - PrintMessage(1, "Only tet elements are supported to split around free segments"); + PrintMessage(3, "Only tet elements are supported to split around free segments"); return; } if(el.IsDeleted()) { - PrintMessage(1,"Element to split is already deleted"); + PrintMessage(3,"Element to split is already deleted"); return; } + auto pi_new = mesh.AddPoint(p_new); + auto seg_new0 = seg; + auto seg_new1 = seg; + seg_new0[1] = pi_new; + seg_new1[0] = pi_new; + + mesh[segi][0] = PointIndex::INVALID; + mesh.AddSegment(seg_new0); + mesh.AddSegment(seg_new1); + + int pmap[4][4] = { {0,1,2,4}, {1,3,2,4}, @@ -897,7 +909,7 @@ namespace netgen for ([[maybe_unused]] auto i : Range(3)) { optmesh.ImproveMesh(); - optmesh.SwapImprove2 (); + optmesh.SwapImprove2(true); optmesh.ImproveMesh(); optmesh.SwapImprove(); optmesh.ImproveMesh(); @@ -910,7 +922,7 @@ namespace netgen auto bad_segs = get_nonconforming(p2el); if(bad_segs.Size() > 0) { - auto bad_seg = mesh[free_segs[bad_segs[0]]]; + auto bad_seg = mesh[bad_segs[0]]; if(debugparam.write_mesh_on_error) mesh.Save("free_segment_not_conformed_dom_"+ToString(domain)+"_seg_"+ToString(bad_seg[0])+"_"+ToString(bad_seg[1])+".vol.gz"); throw Exception("Segment not resolved in volume mesh in domain " + ToString(domain)+ ", seg: " + ToString(bad_seg));