From 0daeeb20aa63af1dc41b3543a4afd26d7845d054 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Thu, 3 Oct 2024 18:06:34 +0200 Subject: [PATCH] BLayer - Option to build only volelements without new surface elements/edges/domains --- libsrc/meshing/boundarylayer.cpp | 387 +++++-------------- libsrc/meshing/boundarylayer.hpp | 5 +- libsrc/meshing/boundarylayer2d.cpp | 2 +- libsrc/meshing/boundarylayer_interpolate.cpp | 10 +- libsrc/meshing/boundarylayer_limiter.hpp | 158 ++++---- libsrc/meshing/debugging.cpp | 2 +- libsrc/meshing/improve3.cpp | 1 + libsrc/meshing/meshfunc.cpp | 3 +- libsrc/meshing/meshing.hpp | 2 +- libsrc/meshing/meshtype.cpp | 22 +- libsrc/meshing/meshtype.hpp | 18 +- libsrc/meshing/python_mesh.cpp | 12 +- ng/ngpkg.cpp | 1 + 13 files changed, 235 insertions(+), 388 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index c2dc476f..8c72f219 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -145,7 +145,7 @@ Vec<3> BoundaryLayerTool ::getEdgeTangent(PointIndex pi, int edgenr, cout << pts << endl; for (auto *p_seg : segs) cout << *p_seg << endl; - throw Exception("Something went wrong in getEdgeTangent!"); + throw NG_EXCEPTION("Something went wrong in getEdgeTangent!"); } tangent = mesh[pts[1]] - mesh[pts[0]]; return tangent.Normalize(); @@ -253,36 +253,8 @@ BoundaryLayerTool::BoundaryLayerTool(Mesh &mesh_, static Timer timer("BoundaryLayerTool::ctor"); RegionTimer regt(timer); ProcessParameters(); - - // for(auto & seg : mesh.LineSegments()) - // seg.edgenr = seg.epgeominfo[1].edgenr; - - total_height = 0.0; - for (auto h : par_heights) - total_height += h; - - max_edge_nr = -1; - for (const auto &seg : mesh.LineSegments()) - if (seg.edgenr > max_edge_nr) - max_edge_nr = seg.edgenr; - - int ndom = mesh.GetNDomains(); - ndom_old = ndom; - - new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); - new_mat_nrs = -1; - for (auto [bcname, matname] : par_new_mat) { - mesh.SetMaterial(++ndom, matname); - regex pattern(bcname); - for (auto i : Range(1, mesh.GetNFD() + 1)) { - auto &fd = mesh.GetFaceDescriptor(i); - if (regex_match(fd.GetBCName(), pattern)) - new_mat_nrs[i] = ndom; - } - } - - if (!params.outside) - domains.Invert(); + if (domains.NumSet() == 0) + return; topo.SetBuildVertex2Element(true); mesh.UpdateTopology(); @@ -321,15 +293,17 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() { isIn != domains.Test(fd.DomainOut())) { int new_si = mesh.GetNFD() + 1; surfacefacs[i] = isIn ? 1. : -1.; - // -1 surf nr is so that curving does not do anything - FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(), - isIn ? fd.DomainOut() : new_mat_nrs[i], -1); - new_fd.SetBCProperty(new_si); - new_fd.SetSurfColour(fd.SurfColour()); - mesh.AddFaceDescriptor(new_fd); - si_map[i] = new_si; moved_surfaces.SetBit(i); - mesh.SetBCName(new_si - 1, "mapped_" + name); + if (!insert_only_volume_elements) { + // -1 surf nr is so that curving does not do anything + FaceDescriptor new_fd(-1, isIn ? new_mat_nrs[i] : fd.DomainIn(), + isIn ? fd.DomainOut() : new_mat_nrs[i], -1); + new_fd.SetBCProperty(new_si); + new_fd.SetSurfColour(fd.SurfColour()); + mesh.AddFaceDescriptor(new_fd); + si_map[i] = new_si; + mesh.SetBCName(new_si - 1, "mapped_" + name); + } } // curving of surfaces with boundary layers will often // result in pushed through elements, since we do not (yet) @@ -347,6 +321,8 @@ void BoundaryLayerTool ::CreateNewFaceDescriptors() { } void BoundaryLayerTool ::CreateFaceDescriptorsSides() { + if (insert_only_volume_elements) + return; BitArray face_done(mesh.GetNFD() + 1); face_done.Clear(); for (const auto &sel : mesh.SurfaceElements()) { @@ -717,7 +693,8 @@ void BoundaryLayerTool ::InsertNewElements( sel.GeomInfo()[i].v = 0.0; } sel.SetIndex(si_map[segj.si]); - mesh.AddSurfaceElement(sel); + new_sels.Append(sel); + new_sels_on_moved_bnd.Append(sel); // TODO: Too many, would be enough to only add outermost ones Segment s1; @@ -771,8 +748,6 @@ void BoundaryLayerTool ::InsertNewElements( BitArray fixed_points(np + 1); fixed_points.Clear(); - BitArray moveboundarypoint(np + 1); - moveboundarypoint.Clear(); auto p2el = mesh.CreatePoint2ElementTable(); for (SurfaceElementIndex si = 0; si < nse; si++) { // copy because surfaceels array will be resized! @@ -820,44 +795,7 @@ void BoundaryLayerTool ::InsertNewElements( for (auto i : Range(points)) newel[i] = newPoint(sel[i], -1, groups[i]); newel.SetIndex(si_map[sel.GetIndex()]); - mesh.AddSurfaceElement(newel); - - // also move volume element adjacent to this surface element accordingly - ElementIndex ei = -1; - // if(groups[0] || groups[1] || groups[2]) - // for(auto ei_ : p2el[sel.PNums()[0]]) - // { - // const auto & el = mesh[ei_]; - // // if(!domains.Test(el.GetIndex())) continue; - // cout << "check " << ei_ << "\t" << el << "\t" << sel << endl; - // auto pnums = el.PNums(); - // if(pnums.Contains(sel[1]) && pnums.Contains(sel[2])) { - // ei = ei_; - // break; - // } - // } - if (ei != -1) { - auto &el = mesh[ei]; - for (auto i : Range(el.GetNP())) - for (auto j : Range(3)) { - if (groups[j] && el[i] == sel[j]) { - el[i] = newel[j]; - break; - } - } - } - } else { - bool has_moved = false; - for (auto p : sel.PNums()) - has_moved |= hasMoved(p); - if (has_moved) - for (auto p : sel.PNums()) { - if (hasMoved(p)) { - fixed_points.SetBit(p); - if (is_boundary_moved.Test(sel.GetIndex())) - moveboundarypoint.SetBit(p); - } - } + new_sels.Append(newel); } if (is_boundary_moved.Test(sel.GetIndex())) { for (auto &p : mesh[si].PNums()) @@ -866,25 +804,27 @@ void BoundaryLayerTool ::InsertNewElements( } } - for (SegmentIndex sei = 0; sei < nseg; sei++) { - auto &seg = segments[sei]; - if (is_boundary_moved.Test(seg.si)) - for (auto &p : seg.PNums()) - if (hasMoved(p)) - p = newPoint(p); - // else if(hasMoved(seg[0]) || hasMoved(seg[1])) - // { - // auto tangent = mesh[seg[1]] - mesh[seg[0]]; - // if(hasMoved(seg[0]) && growthvectors[seg[0]] * tangent > 0) - // seg[0] = newPoint(seg[0]); - // if(hasMoved(seg[1]) && growthvectors[seg[1]] * tangent < 0) - // seg[1] = newPoint(seg[1]); - // } - } + // for (SegmentIndex sei = 0; sei < nseg; sei++) { + // auto &seg = segments[sei]; + // if (is_boundary_moved.Test(seg.si)) { + // if(insert_only_volume_elements) throw + // Exception("insert_only_volume_elements and is_boundary_moved are + // incompatible"); for (auto &p : seg.PNums()) + // if (hasMoved(p)) + // p = newPoint(p); + // } + // } // fill holes in surface mesh at special boundary points (i.e. points with >=4 // adjacent boundary faces) - auto p2sel = mesh.CreatePoint2SurfaceElementTable(); + auto p2sel = ngcore::CreateSortedTable( + new_sels.Range(), + [&](auto &table, SurfaceElementIndex ei) { + for (PointIndex pi : new_sels[ei].PNums()) + table.Add(pi, ei); + }, + mesh.GetNP()); + for (auto &[special_pi, special_point] : special_boundary_points) { if (special_point.growth_groups.Size() != 2) throw Exception("special_point.growth_groups.Size() != 2"); @@ -921,7 +861,7 @@ void BoundaryLayerTool ::InsertNewElements( for (auto &pi : sel.PNums()) if (pi != pi_common && pi != new_special_pi0) pi = new_special_pi1; - mesh.AddSurfaceElement(sel); + new_sels.Append(sel); } } } @@ -960,155 +900,17 @@ void BoundaryLayerTool ::InsertNewElements( sel[1] = seg[0]; sel[2] = pi_new_other; sel.SetIndex(face); - mesh.AddSurfaceElement(sel); + new_sels.Append(sel); } } } } } - - for (ElementIndex ei = 0; ei < ne; ei++) { - auto el = mesh[ei]; - ArrayMem fixed; - ArrayMem moved; - bool moved_bnd = false; - for (const auto &p : el.PNums()) { - if (fixed_points.Test(p)) - fixed.Append(p); - if (hasMoved(p)) - moved.Append(p); - if (moveboundarypoint.Test(p)) - moved_bnd = true; - } - - bool do_move, do_insert; - if (changed_domains.Test(el.GetIndex())) { - do_move = fixed.Size() && moved_bnd; - do_insert = do_move; - } else { - do_move = !fixed.Size() || moved_bnd; - do_insert = !do_move; - } - - // if (do_move) { - // for (auto& p : mesh[ei].PNums()) - // if (hasMoved(p)) { - // if (special_boundary_points.count(p)) { - // auto& special_point = special_boundary_points[p]; - // auto& group = special_point.growth_groups[0]; - // p = group.new_points.Last(); - // } else - // p = newPoint(p); - // } - // } - if (do_insert) { - if (el.GetType() == TET) { - if (moved.Size() == 3) // inner corner - { - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - PointIndex p3 = moved[2]; - auto v1 = mesh[p1]; - auto n = Cross(mesh[p2] - v1, mesh[p3] - v1); - auto d = mesh[newPoint(p1, 0)] - v1; - if (n * d > 0) - Swap(p2, p3); - PointIndex p4 = p1; - PointIndex p5 = p2; - PointIndex p6 = p3; - for (auto i : Range(par_heights)) { - Element nel(PRISM); - nel[0] = p4; - nel[1] = p5; - nel[2] = p6; - p4 = newPoint(p1, i); - p5 = newPoint(p2, i); - p6 = newPoint(p3, i); - nel[3] = p4; - nel[4] = p5; - nel[5] = p6; - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - } - } - if (moved.Size() == 2) { - if (fixed.Size() == 1) { - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - for (auto i : Range(par_heights)) { - PointIndex p3 = newPoint(moved[1], i); - PointIndex p4 = newPoint(moved[0], i); - Element nel(PYRAMID); - nel[0] = p1; - nel[1] = p2; - nel[2] = p3; - nel[3] = p4; - nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - - moved[1]; - if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) * - (mesh[nel[4]] - mesh[nel[1]]) > - 0) - Swap(nel[1], nel[3]); - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - p1 = p4; - p2 = p3; - } - } - } - if (moved.Size() == 1 && fixed.Size() == 1) { - PointIndex p1 = moved[0]; - for (auto i : Range(par_heights)) { - Element nel = el; - PointIndex p2 = newPoint(moved[0], i); - for (auto &p : nel.PNums()) { - if (p == moved[0]) - p = p1; - else if (p == fixed[0]) - p = p2; - } - p1 = p2; - mesh.AddVolumeElement(nel); - } - } - } else if (el.GetType() == PYRAMID) { - if (moved.Size() == 2) { - if (fixed.Size() != 2) - throw Exception("This case is not implemented yet! Fixed size = " + - ToString(fixed.Size())); - PointIndex p1 = moved[0]; - PointIndex p2 = moved[1]; - for (auto i : Range(par_heights)) { - PointIndex p3 = newPoint(moved[1], i); - PointIndex p4 = newPoint(moved[0], i); - Element nel(PYRAMID); - nel[0] = p1; - nel[1] = p2; - nel[2] = p3; - nel[3] = p4; - nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - - fixed[1] - moved[0] - moved[1]; - if (Cross(mesh[p2] - mesh[p1], mesh[p4] - mesh[p1]) * - (mesh[nel[4]] - mesh[nel[1]]) > - 0) - Swap(nel[1], nel[3]); - nel.SetIndex(el.GetIndex()); - mesh.AddVolumeElement(nel); - p1 = p4; - p2 = p3; - } - } else if (moved.Size() == 1) - throw Exception("This case is not implemented yet!"); - } else if (do_move) { - throw Exception( - "Boundarylayer only implemented for tets and pyramids outside " - "yet!"); - } - } - } } void BoundaryLayerTool ::SetDomInOut() { + if (insert_only_volume_elements) + return; for (auto i : Range(1, nfd_old + 1)) if (moved_surfaces.Test(i)) { if (auto dom = mesh.GetFaceDescriptor(si_map[i]).DomainIn(); @@ -1121,6 +923,8 @@ void BoundaryLayerTool ::SetDomInOut() { } void BoundaryLayerTool ::SetDomInOutSides() { + if (insert_only_volume_elements) + return; BitArray done(mesh.GetNFD() + 1); done.Clear(); for (auto sei : Range(mesh.SurfaceElements())) { @@ -1146,6 +950,8 @@ void BoundaryLayerTool ::SetDomInOutSides() { } void BoundaryLayerTool ::AddSegments() { + if (insert_only_volume_elements) + return; if (have_single_segments) MergeAndAddSegments(mesh, segments, new_segments); else { @@ -1155,45 +961,11 @@ void BoundaryLayerTool ::AddSegments() { } } -void BoundaryLayerTool ::FixVolumeElements() { - static Timer timer("BoundaryLayerTool::FixVolumeElements"); - RegionTimer rt(timer); - BitArray is_inner_point(mesh.GetNP() + 1); - is_inner_point.Clear(); - - auto changed_domains = domains; - if (!params.outside) - changed_domains.Invert(); - - for (ElementIndex ei : Range(ne)) - if (changed_domains.Test(mesh[ei].GetIndex())) - for (auto pi : mesh[ei].PNums()) - if (mesh[pi].Type() == INNERPOINT) - is_inner_point.SetBit(pi); - - Array points; - for (auto pi : mesh.Points().Range()) - if (is_inner_point.Test(pi)) - points.Append(pi); - - auto p2el = mesh.CreatePoint2ElementTable(is_inner_point); - - // smooth growth vectors to shift additional element layers to the inside and - // fix flipped tets - for ([[maybe_unused]] auto step : Range(0)) { - for (auto pi : points) { - Vec<3> average_gw = 0.0; - auto &els = p2el[pi]; - size_t cnt = 0; - for (auto ei : els) - if (ei < ne) - for (auto pi1 : mesh[ei].PNums()) - if (pi1 <= np) { - average_gw += growthvectors[pi1]; - cnt++; - } - growthvectors[pi] = 1.0 / cnt * average_gw; - } +void BoundaryLayerTool ::AddSurfaceElements() { + for (auto &sel : + insert_only_volume_elements ? new_sels_on_moved_bnd : new_sels) { + cout << "add surface element " << sel << endl; + mesh.AddSurfaceElement(sel); } } @@ -1242,10 +1014,14 @@ void BoundaryLayerTool ::ProcessParameters() { for (auto id : surfids) par_surfid.Append(id); } - if (string *mat = get_if(¶ms.new_material); mat) - par_new_mat = {{".*", *mat}}; - else - par_new_mat = *get_if>(¶ms.new_material); + + insert_only_volume_elements = !params.new_material.has_value(); + if (params.new_material) { + if (string *mat = get_if(&*params.new_material); mat) + par_new_mat = {{".*", *mat}}; + else + par_new_mat = *get_if>(&*params.new_material); + } if (params.project_boundaries.has_value()) { auto proj_bnd = *params.project_boundaries; @@ -1282,6 +1058,53 @@ void BoundaryLayerTool ::ProcessParameters() { for (auto i : *get_if>(¶ms.domain)) domains.SetBit(i); } + if (domains.NumSet() == 0) + return; + total_height = 0.0; + for (auto h : par_heights) + total_height += h; + + max_edge_nr = -1; + for (const auto &seg : mesh.LineSegments()) + if (seg.edgenr > max_edge_nr) + max_edge_nr = seg.edgenr; + + int ndom = mesh.GetNDomains(); + ndom_old = ndom; + + new_mat_nrs.SetSize(mesh.FaceDescriptors().Size() + 1); + new_mat_nrs = -1; + if (insert_only_volume_elements) { + for (auto i : Range(1, mesh.GetNFD() + 1)) { + auto &fd = mesh.GetFaceDescriptor(i); + auto domin = fd.DomainIn(); + auto domout = fd.DomainOut(); + for (int dom : {domin, domout}) + if (domains.Test(dom)) { + if (params.outside) { + dom = domin + domout - dom; + if (dom == 0) + throw NG_EXCEPTION("No new material specified for boundarylayer " + "on the outside of domain"); + } + new_mat_nrs[i] = dom; + } + } + } else { + for (auto [bcname, matname] : par_new_mat) { + mesh.SetMaterial(++ndom, matname); + regex pattern(bcname); + for (auto i : Range(1, mesh.GetNFD() + 1)) { + auto &fd = mesh.GetFaceDescriptor(i); + if (regex_match(fd.GetBCName(), pattern)) + new_mat_nrs[i] = ndom; + } + } + } + cout << "new mat numbers " << endl << new_mat_nrs << endl; + + if (!params.outside) + domains.Invert(); } void BoundaryLayerTool ::Perform() { @@ -1296,6 +1119,7 @@ void BoundaryLayerTool ::Perform() { SetDomInOut(); AddSegments(); + AddSurfaceElements(); mesh.CalcSurfacesOfNode(); topo.SetBuildVertex2Element(true); @@ -1311,6 +1135,7 @@ void BoundaryLayerTool ::Perform() { mesh[pi] += height * (*gw); } + mesh.CalcSurfacesOfNode(); mesh.GetTopology().ClearEdges(); mesh.SetNextMajorTimeStamp(); mesh.UpdateTopology(); diff --git a/libsrc/meshing/boundarylayer.hpp b/libsrc/meshing/boundarylayer.hpp index f83a402c..2430182d 100644 --- a/libsrc/meshing/boundarylayer.hpp +++ b/libsrc/meshing/boundarylayer.hpp @@ -49,6 +49,7 @@ class BoundaryLayerTool MeshTopology & topo; BoundaryLayerParameters params; Array, PointIndex> growthvectors; + std::map> non_bl_growth_vectors; Table p2sel; BitArray domains, is_edge_moved, is_boundary_projected, is_boundary_moved; @@ -62,11 +63,13 @@ class BoundaryLayerTool // These parameters are derived from given BoundaryLayerParameters and the Mesh Array par_heights; Array par_surfid; + bool insert_only_volume_elements; map par_new_mat; Array par_project_boundaries; bool have_single_segments; Array segments, new_segments; + Array new_sels, new_sels_on_moved_bnd; Array, PointIndex> mapto; Array mapfrom; @@ -91,7 +94,7 @@ class BoundaryLayerTool void SetDomInOut(); void SetDomInOutSides(); void AddSegments(); - void FixVolumeElements(); + void AddSurfaceElements(); Vec<3> getNormal(const Element2d & el) { diff --git a/libsrc/meshing/boundarylayer2d.cpp b/libsrc/meshing/boundarylayer2d.cpp index 93f0cf36..77f0083d 100644 --- a/libsrc/meshing/boundarylayer2d.cpp +++ b/libsrc/meshing/boundarylayer2d.cpp @@ -1,5 +1,5 @@ #include -#include "meshing.hpp" +#include "boundarylayer.hpp" #include "meshing2.hpp" #include "../geom2d/csg2d.hpp" diff --git a/libsrc/meshing/boundarylayer_interpolate.cpp b/libsrc/meshing/boundarylayer_interpolate.cpp index 03ebde05..f7b0655b 100644 --- a/libsrc/meshing/boundarylayer_interpolate.cpp +++ b/libsrc/meshing/boundarylayer_interpolate.cpp @@ -2,10 +2,6 @@ namespace netgen { -// TODO: Hack, move this to the header or restructure the whole growth_vectors -// storage -static std::map> non_bl_growth_vectors; - void BoundaryLayerTool ::InterpolateGrowthVectors() { int new_max_edge_nr = max_edge_nr; for (const auto &seg : segments) @@ -113,10 +109,10 @@ void BoundaryLayerTool ::InterpolateGrowthVectors() { if (plast != seg[0] && plast != seg[1]) continue; auto pnew = plast == seg[0] ? seg[1] : seg[0]; - if(pnew == points[0] && points.Size()>1) { - + if (pnew == points[0] && points.Size() > 1) { } - if (points_set.count(pnew) > 0 && (pnew != points[0] || points.Size()==2)) + if (points_set.count(pnew) > 0 && + (pnew != points[0] || points.Size() == 2)) continue; edge_len += (mesh[points.Last()] - mesh[pnew]).Length(); points.Append(pnew); diff --git a/libsrc/meshing/boundarylayer_limiter.hpp b/libsrc/meshing/boundarylayer_limiter.hpp index 3c32be8c..19f0f6c6 100644 --- a/libsrc/meshing/boundarylayer_limiter.hpp +++ b/libsrc/meshing/boundarylayer_limiter.hpp @@ -25,17 +25,31 @@ struct GrowthVectorLimiter { GrowthVectorLimiter(BoundaryLayerTool &tool_) : tool(tool_), params(tool_.params), mesh(tool_.mesh), height(tool_.total_height), growthvectors(tool_.growthvectors), - map_from(mesh.Points().Size()), - p2sel(mesh.CreatePoint2SurfaceElementTable()) { + map_from(mesh.Points().Size()) { changed_domains = tool.domains; if (!params.outside) changed_domains.Invert(); map_from = tool.mapfrom; + p2sel = ngcore::CreateSortedTable( + tool.new_sels.Range(), + [&](auto &table, SurfaceElementIndex ei) { + for (PointIndex pi : tool.new_sels[ei].PNums()) + table.Add(pi, ei); + }, + mesh.GetNP()); + } + + auto SurfaceElementsRange() { return Range(tool.nse + tool.new_sels.Size()); } + + const auto &Get(SurfaceElementIndex sei) { + if (sei < tool.nse) + return mesh[sei]; + return tool.new_sels[sei - tool.nse]; } std::pair GetMinMaxLimit(SurfaceElementIndex sei) { - const auto &sel = mesh[sei]; + const auto &sel = Get(sei); double min_limit = GetLimit(sel[0]); double max_limit = min_limit; for (auto i : IntRange(1, sel.GetNP())) { @@ -75,7 +89,7 @@ struct GrowthVectorLimiter { Point<3> GetPoint(PointIndex pi_to, double shift = 1., bool apply_limit = false) { - if (tool.growth_vector_map.count(pi_to) == 0) + if (pi_to <= tool.np || tool.growth_vector_map.count(pi_to) == 0) return mesh[pi_to]; return mesh[pi_to] + GetVector(pi_to, shift, apply_limit); @@ -97,7 +111,7 @@ struct GrowthVectorLimiter { auto GetTrig(SurfaceElementIndex sei, double shift = 0.0, bool apply_limit = false) { - auto sel = mesh[sei]; + auto sel = Get(sei); std::array, 3> trig; for (auto i : Range(3)) trig[i] = GetPoint(sel[i], shift, apply_limit); @@ -105,7 +119,7 @@ struct GrowthVectorLimiter { } auto GetMappedTrig(SurfaceElementIndex sei, double shift = 0.0) { - auto sel = mesh[sei]; + auto sel = Get(sei); std::array, 3> trig; for (auto i : Range(3)) trig[i] = GetMappedPoint(sel[i], shift); @@ -115,7 +129,7 @@ struct GrowthVectorLimiter { auto GetSideTrig(SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first_vertex = true) { auto trig = GetMappedTrig(sei, 0.0); - auto sel = mesh[sei]; + auto sel = Get(sei); auto index1 = (index + 1) % 3; if (!grow_first_vertex) index1 = (index + 2) % 3; @@ -133,7 +147,7 @@ struct GrowthVectorLimiter { auto seg = GetSeg(pi_to, seg_shift, true); - for (auto pi : mesh[sei].PNums()) { + for (auto pi : Get(sei).PNums()) { if (pi == pi_from) return false; if (map_from[pi] == pi_from) @@ -172,7 +186,7 @@ struct GrowthVectorLimiter { double max_limit = max(GetLimit(pi_to), trig_max_limit); bool result = false; result = SetLimit(pi_to, s * max_limit); - for (auto pi : mesh[sei].PNums()) + for (auto pi : Get(sei).PNums()) result = SetLimit(pi, s * max_limit); return result; } else { @@ -198,8 +212,8 @@ struct GrowthVectorLimiter { return; for (PointIndex pi : IntRange(tool.np, mesh.GetNP())) { std::set pis; - for (auto sel : p2sel[pi]) - for (auto pi_ : mesh[sel].PNums()) + for (auto sei : p2sel[pi]) + for (auto pi_ : Get(sei).PNums()) pis.insert(pi_); ArrayMem limits; for (auto pi : pis) { @@ -230,7 +244,7 @@ struct GrowthVectorLimiter { // from original surface elements if (sei >= tool.nse) return false; - const auto sel = mesh[sei]; + const auto sel = Get(sei); auto np = sel.GetNP(); for (auto i : Range(np)) { if (sel[i] > tool.np) @@ -254,8 +268,6 @@ struct GrowthVectorLimiter { for (SurfaceElementIndex sei : mesh.SurfaceElements().Range()) { auto sel = mesh[sei]; const auto &fd = mesh.GetFaceDescriptor(sel.GetIndex()); - if (sei >= tool.nse) - continue; if (sel.GetNP() == 4) continue; @@ -338,9 +350,9 @@ struct GrowthVectorLimiter { tree = make_unique>(bbox); - for (auto sei : mesh.SurfaceElements().Range()) { - const auto &sel = mesh[sei]; - auto sel_index = mesh[sei].GetIndex(); + for (auto sei : SurfaceElementsRange()) { + const auto &sel = Get(sei); + auto sel_index = sel.GetIndex(); Box<3> box(Box<3>::EMPTY_BOX); for (auto pi : sel.PNums()) { @@ -369,7 +381,7 @@ struct GrowthVectorLimiter { box.Add(GetPoint(pi_to, GetLimit(pi_from))); tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&](SurfaceElementIndex sei) { - const auto &sel = mesh[sei]; + const auto &sel = Get(sei); if (sel.PNums().Contains(pi_from)) return false; if (sel.PNums().Contains(pi_to)) @@ -392,8 +404,8 @@ struct GrowthVectorLimiter { mesh.GetBox(pmin, pmax); BoxTree<3, SurfaceElementIndex> setree(pmin, pmax); - for (auto sei : mesh.SurfaceElements().Range()) { - const Element2d &tri = mesh[sei]; + for (auto sei : SurfaceElementsRange()) { + const Element2d &tri = Get(sei); Box<3> box(Box<3>::EMPTY_BOX); for (PointIndex pi : tri.PNums()) @@ -403,61 +415,71 @@ struct GrowthVectorLimiter { setree.Insert(box, sei); } - for (auto sei : mesh.SurfaceElements().Range()) { - const Element2d &tri = mesh[sei]; + for (auto sei : SurfaceElementsRange()) { + const Element2d &tri = Get(sei); Box<3> box(Box<3>::EMPTY_BOX); for (PointIndex pi : tri.PNums()) box.Add(GetPoint(pi, 1.0, true)); - setree.GetFirstIntersecting( - box.PMin(), box.PMax(), [&](SurfaceElementIndex sej) { - const Element2d &tri2 = mesh[sej]; + setree.GetFirstIntersecting(box.PMin(), box.PMax(), [&](size_t sej) { + const Element2d &tri2 = Get(sej); - if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer()) - return false; + if (mesh[tri[0]].GetLayer() != mesh[tri2[0]].GetLayer()) + return false; - netgen::Point<3> tri1_points[3], tri2_points[3]; - const netgen::Point<3> *trip1[3], *trip2[3]; - for (int k = 0; k < 3; k++) { - trip1[k] = &tri1_points[k]; - trip2[k] = &tri2_points[k]; + netgen::Point<3> tri1_points[3], tri2_points[3]; + const netgen::Point<3> *trip1[3], *trip2[3]; + for (int k = 0; k < 3; k++) { + trip1[k] = &tri1_points[k]; + trip2[k] = &tri2_points[k]; + } + auto set_points = [&]() { + for (int k = 0; k < 3; k++) { + tri1_points[k] = GetPoint(tri[k], 1.0, true); + tri2_points[k] = GetPoint(tri2[k], 1.0, true); + } + }; + + set_points(); + + int counter = 0; + while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) { + changed = true; + PointIndex pi_max_limit = PointIndex::INVALID; + for (PointIndex pi : + {tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]}) + if (pi > tool.np && (!pi_max_limit.IsValid() || + GetLimit(pi) > GetLimit(pi_max_limit))) + pi_max_limit = map_from[pi]; + + if (!pi_max_limit.IsValid()) + break; + + ScaleLimit(pi_max_limit, 0.9); + set_points(); + counter++; + if (counter > 5) { + break; + cerr << "Limit intersecting surface elements: too many " + "limitation steps, sels: " + << Get(sei) << '\t' << Get(sej) << endl; + for (auto si : {sei, sej}) { + auto sel = Get(si); + cerr << "Limits: "; + for (auto pi : sel.PNums()) + cerr << GetLimit(pi) << ",\t"; + cerr << endl; + for (auto pi : sel.PNums()) + cerr << GetPoint(pi, 1.0, true) << "\t"; + cerr << endl; } - auto set_points = [&]() { - for (int k = 0; k < 3; k++) { - tri1_points[k] = GetPoint(tri[k], 1.0, true); - tri2_points[k] = GetPoint(tri2[k], 1.0, true); - } - }; - - set_points(); - - int counter = 0; - while (IntersectTriangleTriangle(&trip1[0], &trip2[0])) { - changed = true; - PointIndex pi_max_limit = PointIndex::INVALID; - for (PointIndex pi : - {tri[0], tri[1], tri[2], tri2[0], tri2[1], tri2[2]}) - if (pi > tool.np && - (!pi_max_limit.IsValid() || - limits[tool.mapfrom[pi]] > limits[pi_max_limit])) - pi_max_limit = tool.mapfrom[pi]; - - if (!pi_max_limit.IsValid()) - break; - - limits[pi_max_limit] *= 0.9; - set_points(); - counter++; - if (counter > 30) { - cerr << "Limit intersecting surface elements: too many " - "limitation steps, sels: " - << mesh[sei] << '\t' << mesh[sej] << endl; - break; - } - } - return false; - }); + cerr << "pi_max_limit " << pi_max_limit << endl; + break; + } + } + return false; + }); } } } @@ -491,7 +513,7 @@ struct GrowthVectorLimiter { [&](PointIndex pi_to, SurfaceElementIndex sei) { if (LimitGrowthVector(pi_to, sei, trig_shift, seg_shift)) limit_counter++; - auto sel = mesh[sei]; + auto sel = Get(sei); bool is_mapped = true; for (auto pi : sel.PNums()) { if (pi >= tool.np) diff --git a/libsrc/meshing/debugging.cpp b/libsrc/meshing/debugging.cpp index 949bb45c..c5dbd906 100644 --- a/libsrc/meshing/debugging.cpp +++ b/libsrc/meshing/debugging.cpp @@ -36,13 +36,13 @@ namespace netgen mesh->SetMaterial(2, "2_points"); mesh->SetMaterial(3, "3_points"); mesh->SetMaterial(4, "4_points"); - mesh->Compress(); mesh->ClearSurfaceElements(); for (auto & el : openelements) mesh->AddSurfaceElement( el ); + mesh->Compress(); return mesh; } diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index ac6fbb4d..955431ed 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -2637,6 +2637,7 @@ double MeshOptimize3d :: SplitImprove2Element ( Element & elem = mesh[ei0]; if (elem.IsDeleted()) return false; if (ei0 == ei) continue; + if (elem.GetType() != TET) return false; if (elem[0] == pi1 || elem[1] == pi1 || elem[2] == pi1 || elem[3] == pi1 || (elem.GetNP()==5 && elem[4]==pi1) ) if(!has_both_points0.Contains(ei0)) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index fc9ca918..4b5116bd 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -3,6 +3,7 @@ #include #include "meshing.hpp" #include "debugging.hpp" +#include "boundarylayer.hpp" namespace netgen { @@ -421,7 +422,7 @@ namespace netgen mesh.FindOpenElements(domain); PrintMessage (5, mesh.GetNOpenElements(), " open faces"); - // GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol"); + // GetOpenElements( mesh, domain )->Save("open_"+ToString(domain)+"_"+ToString(cntsteps)+".vol"); cntsteps++; diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp index 779453b8..e9363ce1 100644 --- a/libsrc/meshing/meshing.hpp +++ b/libsrc/meshing/meshing.hpp @@ -53,7 +53,7 @@ namespace netgen #include "bisect.hpp" #include "hprefinement.hpp" -#include "boundarylayer.hpp" +// #include "boundarylayer.hpp" #include "specials.hpp" #include "validate.hpp" #include "basegeom.hpp" diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 7f49d334..728ac06a 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2928,14 +2928,17 @@ namespace netgen case 1: print_vec(std::get<1>(mp.thickness)); break; } ost <<"\n new_material: "; - switch(mp.new_material.index()) - { - case 0: ost << std::get<0>(mp.new_material); break; - case 1: - for (const auto & [key, value] : std::get<1>(mp.new_material)) - ost << key << " -> " << value << ", "; - break; - } + if(!mp.new_material) ost << "nullopt"; + else { + switch(mp.new_material->index()) + { + case 0: ost << std::get<0>(*mp.new_material); break; + case 1: + for (const auto & [key, value] : std::get<1>(*mp.new_material)) + ost << key << " -> " << value << ", "; + break; + } + } ost << "\n domain: "; switch(mp.domain.index()) { @@ -2957,9 +2960,8 @@ namespace netgen ost << "nullopt"; ost << "\n grow_edges: " << mp.grow_edges; ost << "\n limit_growth_vectors: " << mp.limit_growth_vectors; - ost << "\n sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex; + ost << "\n sides_keep_surfaceindex: " << (mp.sides_keep_surfaceindex ? ToString(*mp.sides_keep_surfaceindex) : "nullopt"); ost << "\n keep_surfaceindex: " << mp.keep_surfaceindex; - ost << "\n limit_safety: " << mp.limit_safety; ost << endl; return ost; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 178dd471..926c1411 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1276,18 +1276,16 @@ namespace netgen struct BoundaryLayerParameters { - std::variant> boundary; std::variant> thickness; - std::variant> new_material; std::variant> domain; - bool outside; - std::optional>> project_boundaries; - bool grow_edges; - bool limit_growth_vectors; - bool sides_keep_surfaceindex; - bool keep_surfaceindex; - - double limit_safety = 0.3; // alloow only 30% of the growth vector length + std::variant> boundary = ".*"; + std::optional>> new_material = nullopt; + std::optional>> project_boundaries = nullopt; + bool outside = false; + bool grow_edges = true; + bool limit_growth_vectors = true; + std::optional sides_keep_surfaceindex = nullopt; // !outside by default + bool keep_surfaceindex = false; }; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 0500ac5c..b36d03a6 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -8,6 +8,7 @@ #include #include "meshing.hpp" +#include "boundarylayer.hpp" // #include // #include #include <../interface/rw_medit.hpp> @@ -1730,9 +1731,8 @@ py::arg("point_tolerance") = -1.) std::optional>> project_boundaries, bool grow_edges, bool limit_growth_vectors, - bool sides_keep_surfaceindex, - bool keep_surfaceindex, - double limit_safety) + std::optional sides_keep_surfaceindex, + bool keep_surfaceindex) { BoundaryLayerParameters blp; blp.boundary = boundary; @@ -1745,14 +1745,13 @@ py::arg("point_tolerance") = -1.) blp.limit_growth_vectors = limit_growth_vectors; blp.sides_keep_surfaceindex = sides_keep_surfaceindex; blp.keep_surfaceindex = keep_surfaceindex; - blp.limit_safety = limit_safety; return blp; }), py::arg("boundary"), py::arg("thickness"), py::arg("new_material"), py::arg("domain") = ".*", py::arg("outside") = false, py::arg("project_boundaries")=nullopt, py::arg("grow_edges")=true, - py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=false, - py::arg("keep_surfaceindex")=false, py::arg("limit_safety")=0.3, + py::arg("limit_growth_vectors") = true, py::arg("sides_keep_surfaceindex")=nullopt, + py::arg("keep_surfaceindex")=false, R"delimiter( Add boundary layer to mesh. @@ -1805,7 +1804,6 @@ project_boundaries : Optional[str] = None .def_readwrite("limit_growth_vectors", &BoundaryLayerParameters::limit_growth_vectors) .def_readwrite("sides_keep_surfaceindex", &BoundaryLayerParameters::sides_keep_surfaceindex) .def_readwrite("keep_surfaceindex", &BoundaryLayerParameters::keep_surfaceindex) - .def_readwrite("limit_safety", &BoundaryLayerParameters::limit_safety) ; py::implicitly_convertible(); diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index 8b6a245b..7d39ba37 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -9,6 +9,7 @@ The interface between the GUI and the netgen library #include #include +#include "../libsrc/meshing/boundarylayer.hpp" #include