diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index e0b3ec76..5908aa7e 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -1476,14 +1476,23 @@ struct GrowthVectorLimiter { mapto.SetSize(0); mapto.SetSize(np); + auto & identifications = mesh.GetIdentifications(); + const int identnr = identifications.GetNr("boundarylayer"); + identifications.SetType(identnr, Identifications::CLOSESURFACES); + auto add_points = [&](PointIndex pi, Vec<3> & growth_vector, Array & new_points) { Point<3> p = mesh[pi]; + PointIndex pi_last = pi; for(auto i : Range(params.heights)) { // p += params.heights[i] * growth_vector; - new_points.Append(mesh.AddPoint(p)); - growth_vector_map[new_points.Last()] = { &growth_vector, params.heights[i] }; + auto pi_new = mesh.AddPoint(p); + new_points.Append(pi_new); + growth_vector_map[pi_new] = { &growth_vector, params.heights[i] }; + if(special_boundary_points.count(pi) == 0) + mesh.GetIdentifications().Add(pi_last, pi_new, identnr); + pi_last = pi_new; } }; @@ -1660,14 +1669,34 @@ struct GrowthVectorLimiter { } } + auto getClosestGroup = [&] (PointIndex pi, SurfaceElementIndex sei) + { + auto n = numGroups(pi); + if(n == 1) return 0; + const auto & sel = mesh[sei]; + auto igroup = 0; + double distance = 1e99; + for(auto j : Range(n)) { + auto g = getGroups(pi, sel.GetIndex()); + auto vcenter = Center(mesh[sel[0]], mesh[sel[1]], mesh[sel[2]]); + auto dist = (vcenter-(mesh[pi]+special_boundary_points[pi].growth_groups[j].growth_vector)).Length2(); + if(dist < distance) { + distance = dist; + igroup = j; + } + } + return igroup; + }; + 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! - auto sel = mesh[si]; + const auto sel = mesh[si]; if(moved_surfaces.Test(sel.GetIndex())) { Array points(sel.PNums()); @@ -1696,7 +1725,13 @@ struct GrowthVectorLimiter { cout << getGroups(sel[i], sel.GetIndex()) << endl; } groups[i] = getGroups(sel[i], sel.GetIndex())[igroup]; + // groups[i] = getClosestGroup(sel[i], sel.GetIndex()); + } + bool add_volume_element = true; + for(auto pi : sel.PNums()) + if(numGroups(pi)>1) + add_volume_element = false; for(auto j : Range(params.heights)) { auto eltype = points.Size() == 3 ? PRISM : HEX; @@ -1709,13 +1744,42 @@ struct GrowthVectorLimiter { for(auto i : Range(points)) el[sel.PNums().Size() + i] = points[i]; el.SetIndex(new_mat_nrs[sel.GetIndex()]); - mesh.AddVolumeElement(el); + if(add_volume_element) + mesh.AddVolumeElement(el); } Element2d newel = sel; 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) { + cout << "move " << ei << mesh[ei] << endl; + 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; + } + } + cout << "after " << ei << mesh[ei] << endl; + } } else { @@ -1822,8 +1886,15 @@ struct GrowthVectorLimiter { if(do_move) { for(auto& p : mesh[ei].PNums()) - if(hasMoved(p)) - p = newPoint(p); + 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) { @@ -2052,6 +2123,7 @@ struct GrowthVectorLimiter { FixVolumeElements(); + InsertNewElements(segmap, in_surface_direction); SetDomInOut(); @@ -2060,9 +2132,9 @@ struct GrowthVectorLimiter { mesh.CalcSurfacesOfNode(); topo.SetBuildVertex2Element(true); mesh.UpdateTopology(); + InterpolateGrowthVectors(); // cout << "growthvectors before " << endl<< growthvectors << endl; - InterpolateGrowthVectors(); // cout << "growthvectors after " << endl << growthvectors << endl; if(params.limit_growth_vectors) diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 8022c7d7..ff2af770 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -275,6 +275,7 @@ namespace netgen auto & mesh = *md.mesh; auto domain = md.domain; MeshingParameters & mp = md.mp; + cout << "try to close open quads " << md.domain << endl; int oldne; if (multithread.terminate) @@ -282,6 +283,7 @@ namespace netgen mesh.CalcSurfacesOfNode(); mesh.FindOpenElements(domain); + cout << "\thave close open elements " << mesh.GetNOpenElements() << endl; if (!mesh.GetNOpenElements()) return; @@ -292,6 +294,7 @@ namespace netgen if (mesh.HasOpenQuads()) { + cout << "\thave close open quads" << endl; string rulefile = ngdir; const char ** rulep = NULL; @@ -353,6 +356,7 @@ namespace netgen << "mesh has " << mesh.GetNE() << " prism/pyramid elements" << endl; mesh.FindOpenElements(domain); + cout << "\thave close open elements " << mesh.GetNOpenElements() << endl; } } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 1a50a250..298a2951 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1298,7 +1298,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) CreateMPfromKwargs(mp, kwargs); } MeshVolume (mp, self); - OptimizeVolume (mp, self); + // OptimizeVolume (mp, self); }, py::arg("mp")=nullptr, meshingparameter_description.c_str(), py::call_guard())