From 241f6eb9e54af6c3601aa95ed3a9435ac4adddba Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Sun, 24 Mar 2024 20:39:50 +0100 Subject: [PATCH] Some debug output cleanup --- libsrc/meshing/boundarylayer.cpp | 136 +++++++++++++++++++++++-------- libsrc/meshing/debugging.cpp | 58 ++++++++++--- libsrc/meshing/debugging.hpp | 24 +++--- libsrc/meshing/meshfunc.cpp | 2 +- 4 files changed, 159 insertions(+), 61 deletions(-) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 3202a002..a1be8f51 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -216,11 +216,11 @@ struct GrowthVectorLimiter { } auto GetSideTrig(SurfaceElementIndex sei, int index, double shift = 0.0, bool grow_first = true) { - cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; - auto trig = GetTrig(sei, 0.0); - cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; + // cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; + auto trig = GetMappedTrig(sei, 0.0); + // cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; auto sel = mesh[sei]; - cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; + // cout << "get side trig " << sei << "/" << mesh.SurfaceElements().Size() << endl; auto index1 = (index+1)%3; if(!grow_first) index1 = (index+2)%3; trig[index] = GetMappedPoint(sel[index1], shift); @@ -274,8 +274,7 @@ struct GrowthVectorLimiter { } dshift /= 0.9; intersection = isIntersectingTrig(pi_from, pi_to, sei, dshift); - if(dshift < 1) - cout << "final dshift " << dshift << "\t" << intersection.lam0 << endl; + // if(dshift < 1) cout << "final dshift " << dshift << "\t" << intersection.lam0 << endl; limits[pi_from] *= intersection.lam0; auto sel = mesh[sei]; for(auto i : Range(3)) @@ -308,14 +307,15 @@ struct GrowthVectorLimiter { void LimitSelfIntersection () { // check for self-intersection within new elements (prisms/hexes) + bool found_debug_element = false; auto isIntersecting = [&](SurfaceElementIndex sei, double shift) { if(sei>=tool.nse) return false; // ignore new surface elements, side trigs are only built from original surface elements const auto sel = mesh[sei]; - cout << "check " << sei << " " << mesh[sei] << " shift = " << shift << endl; + // cout << "check " << sei << " " << mesh[sei] << " shift = " << shift << endl; auto np = sel.GetNP(); for(auto i : Range(np)) { - if(sel[i] > tool.np) continue; + if(sel[i] > tool.np) return false; if(tool.mapto[sel[i]].Size() == 0) return false; } for(auto i : Range(np)) { @@ -327,17 +327,67 @@ struct GrowthVectorLimiter { // cout << "check " << seg[0] << seg[1] << endl; // cout << "\t" << trig[0] << trig[1] << trig[2] << endl; // } + if(false && sei==104) { + auto isect = isIntersectingPlane(seg, trig); + + Array> points; + if(i==0 && side) { + points.Append(mesh[sel[0]]); + points.Append(mesh[sel[1]]); + points.Append(mesh[sel[2]]); + // debug_gui.DrawTrigs(points, "sel_"+ToString(sei), "green"); + + points.Append(GetMappedPoint(sel[0], shift*limits[sel[0]])); + points.Append(GetMappedPoint(sel[1], shift*limits[sel[1]])); + points.Append(GetMappedPoint(sel[2], shift*limits[sel[2]])); + // debug_gui.DrawTrigs(points, "sel_"+ToString(sei), "green"); + } + + string name = "face_"+ToString(i); + + points.SetSize(0); + points.Append(trig[0]); + points.Append(trig[1]); + points.Append(trig[2]); + // debug_gui.DrawTrigs(points, name, "blue"); + + points.SetSize(0); + points.Append(seg[0]); + points.Append(seg[1]); + // debug_gui.DrawLines(points, name, "red"); + + if(isect) { + points.SetSize(0); + points.Append(isect.p); + // debug_gui.DrawPoints(points, name, "red"); + + + + + } + // debug_gui.DrawLines(points, "points_sel_"+ToString(sei)+"_"+ToString(i)+"_"+ToString(side), "blue"); + + + found_debug_element = true; + + } if(auto isect = isIntersectingPlane(seg, trig); isect) { // if(sel.PNums().Contains(1)) - cout << "found intersection " << shift << "\t"<< sei << " " << mesh[sei] << endl; - cout << isect.lam0 << ',' << isect.p << endl; - cout << "\t" << seg[0] << seg[1] << endl; - cout << "\t" << trig[0] << trig[1] << trig[2] << endl; + // cout << "found intersection " << shift << "\t"<< sei << " " << mesh[sei] << endl; + // cout << isect.lam0 << ',' << isect.p << endl; + // cout << "\t" << seg[0] << seg[1] << endl; + // cout << "\t" << trig[0] << trig[1] << trig[2] << endl; return true; } } } } + + if(false && sei==104) { + // int i_; + // cin >> i_; + // debug_gui.ResetObjects(); + } return false; }; @@ -380,18 +430,18 @@ struct GrowthVectorLimiter { const double step_factor = 0.9; while(shift>0.01 && isIntersecting(sei, shift*safety)) { shift *= step_factor; - cout << "reduce " << sel << "\tshift = " << shift << endl; + // cout << "reduce " << sel << "\tshift = " << shift << endl; double max_limit = 0; - cout << "\tlimits before " << limits[sel[0]] << '\t' << limits[sel[1]] << '\t' << limits[sel[2]] << endl; + // cout << "\tlimits before " << limits[sel[0]] << '\t' << limits[sel[1]] << '\t' << limits[sel[2]] << endl; for(auto i : Range(np)) max_limit = max(max_limit, limits[sel[i]]); for(auto i : Range(np)) if(max_limit == limits[sel[i]]) limits[sel[i]] *= step_factor; - cout << "\tlimits after " << limits[sel[0]] << '\t' << limits[sel[1]] << '\t' << limits[sel[2]] << endl; + // cout << "\tlimits after " << limits[sel[0]] << '\t' << limits[sel[1]] << '\t' << limits[sel[2]] << endl; if(max_limit < 0.01) { // if(sel.PNums().Contains(1)) - cout << "self intersection" << endl; + // cout << "self intersection" << endl; break; } } @@ -487,7 +537,7 @@ struct GrowthVectorLimiter { for(PointIndex pi : mesh.Points().Range()) { bbox.Add(mesh[pi]); - bbox.Add(GetPoint(pi, 1.0)); + bbox.Add(GetPoint(pi, 1.1)); // if(tool.mapto[pi].Size() >0) // bbox.Add(mesh[tool.mapto[pi].Last()]); } @@ -526,27 +576,34 @@ struct GrowthVectorLimiter { auto np_new = mesh.Points().Size(); // cout << "np_new " << np_new << endl; // cout << "too.np " << tool.np << endl; + int counter = 0; for(auto i : IntRange(tool.np, np_new)) { // cout << "handle point " << i << endl; PointIndex pi_to = i+PointIndex::BASE; + PointIndex pi_from = map_from[pi_to]; + if(!pi_from.IsValid()) + throw Exception("Point not mapped"); + // if(mesh[pi_to].Type() == INNERPOINT) // continue; // if(growthvectors[pi_to].Length2() == 0.0) // continue; Box<3> box(Box<3>::EMPTY_BOX); auto seg = GetSeg(pi_to, seg_shift); - box.Add(seg[0]); - box.Add(seg[1]); + + box.Add(GetPoint(pi_to, 0)); + box.Add(GetPoint(pi_to, limits[pi_from])); tree->GetFirstIntersecting(box.PMin(), box.PMax(), [&](SurfaceElementIndex sei) { const auto & sel = mesh[sei]; // cout << "got intersecting " << sei << endl; - auto pi_from = map_from[pi_to]; if(sel.PNums().Contains(pi_from)) return false; + counter++; f(pi_to, sei); return false; }); } + cout << "intersections counter " << counter << endl; } }; @@ -606,13 +663,13 @@ struct GrowthVectorLimiter { GrowthVectorLimiter limiter(*this); //, mesh, params, limits, growthvectors, total_height); // limit to not intersect with other (original) surface elements - // double trig_shift = 0; - // double seg_shift = 2.1; - // limiter.FindTreeIntersections(trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { - // // cout << "found intersection 1 " << pi_to << ", " << sei << endl; - // if(sei >= nse) return; // ignore new surface elements in first pass - // limiter.LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); - // }); + double trig_shift = 0; + double seg_shift = 2.1; + limiter.FindTreeIntersections(trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { + // cout << "found intersection 1 " << pi_to << ", " << sei << endl; + if(sei >= nse) return; // ignore new surface elements in first pass + limiter.LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); + }); limiter.LimitSelfIntersection(); @@ -620,14 +677,14 @@ struct GrowthVectorLimiter { // growthvectors[i] *= limits[i]; // limits = 1.0; - // // now limit again with shifted surace elements - // trig_shift = 1.0; - // seg_shift = 1.0; - // limiter.FindTreeIntersections(trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { - // // cout << "Should not have intersection with original surface element anymore" << endl; - // // cout << "found intersection 2 " << pi_to << ", " << sei << endl; - // limiter.LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); - // }); + // // now limit again with shifted surface elements + trig_shift = 1.1; + seg_shift = 1.1; + limiter.FindTreeIntersections(trig_shift, seg_shift, [&] (PointIndex pi_to, SurfaceElementIndex sei) { + // cout << "Should not have intersection with original surface element anymore" << endl; + // cout << "found intersection 2 " << pi_to << ", " << sei << endl; + limiter.LimitGrowthVector(pi_to, sei, trig_shift, seg_shift); + }); // for (auto i : Range(limits)) // if(limits[i] < 1.0) @@ -639,7 +696,7 @@ struct GrowthVectorLimiter { // limits[pi_from] = min(limits[pi_from], limits[pi_to]); // } - cout << "limits " << endl << limits << endl; + // cout << "limits " << endl << limits << endl; for(auto i : Range(growthvectors)) growthvectors[i] *= limits[i]; @@ -2329,7 +2386,14 @@ struct GrowthVectorLimiter { static Timer timer("Create Boundarylayers"); RegionTimer regt(timer); + cout << "generate boundary layers" << endl; + cout << "points before " << mesh.GetNP() << endl; + cout << "sels before " << mesh.GetNSE() << endl; + + // debug_gui.DrawMesh(mesh); + BoundaryLayerTool tool(mesh, blp); + cout << tool.nse << ',' << tool.np << endl; tool.Perform(); mesh.Save("layer.vol"); } diff --git a/libsrc/meshing/debugging.cpp b/libsrc/meshing/debugging.cpp index 4e3e31b8..0dd1f8d8 100644 --- a/libsrc/meshing/debugging.cpp +++ b/libsrc/meshing/debugging.cpp @@ -192,28 +192,60 @@ void DebuggingGUI::Stop() { cout << "joined thread" << endl; } -void DebuggingGUI::DrawMesh(const string& name, const Mesh& m) { +void DebuggingGUI::DrawMesh(const Mesh& m) { py::gil_scoped_acquire acquire; const auto webgui = py::module::import("netgen.webgui"); const auto dumps = py::module::import("json").attr("dumps"); - auto smesh = make_shared(); - *smesh = m; + mesh = make_shared(); + *mesh = m; const auto py_data = - webgui.attr("Draw")(smesh, py::arg("show") = false).attr("GetData")(); + webgui.attr("Draw")(mesh, py::arg("show") = false).attr("GetData")(); const string d = py::cast(dumps(py_data)); data = json::parse(d); Send(d); } -void DebuggingGUI::DrawPoints(const string& name, const Mesh& m, - FlatArray points) {} -void DebuggingGUI::DrawLines(const string& name, const Mesh& m, - FlatArray lines) {} -void DebuggingGUI::DrawTrigs(const string& name, const Mesh& m, - FlatArray trigs) {} -void DebuggingGUI::DrawTets(const string& name, const Mesh& m, - FlatArray tets) {} -void DebuggingGUI::AddComponent(const Mesh& m) {} +void DebuggingGUI::DrawPoints(FlatArray> position, string name, + string color) { + DrawObject(position, "points", name, color); +} + +void DebuggingGUI::DrawLines(FlatArray> position, string name, + string color) { + DrawObject(position, "lines", name, color); +} +void DebuggingGUI::DrawTrigs(FlatArray> position, string name, + string color) { + // Array> pnts; + // cout << "trigs " << position << endl; + // for (auto i : Range(position.Size()/3)) { + // pnts.Append(position[3*i+0]); + // pnts.Append(position[3*i+1]); + // pnts.Append(position[3*i+1]); + // pnts.Append(position[3*i+2]); + // pnts.Append(position[3*i+2]); + // pnts.Append(position[3*i+0]); + + // } + // cout << "pnts size " << pnts.Size() << endl; + DrawObject(position, "trigs", name, color); +} + +void DebuggingGUI::DrawObject(FlatArray> position, string type, string name, string color) { + json p = json::array_t(); + for (auto pnt : position) { + p.push_back(pnt[0]); + p.push_back(pnt[1]); + p.push_back(pnt[2]); + } + json d; + d["type"] = type; + d["name"] = name; + d["color"] = color; + d["position"] = p; + data["objects"].push_back(d); + Send(data); +} void DebuggingGUI::Send(const string& message) { if (loop) { diff --git a/libsrc/meshing/debugging.hpp b/libsrc/meshing/debugging.hpp index 1060c179..f4f1674d 100644 --- a/libsrc/meshing/debugging.hpp +++ b/libsrc/meshing/debugging.hpp @@ -1,6 +1,6 @@ #include "meshing.hpp" -#define NETGEN_DEBUGGING_GUI +// #define NETGEN_DEBUGGING_GUI #ifdef NETGEN_DEBUGGING_GUI #include "json.hpp" @@ -24,22 +24,24 @@ class DebuggingGUI { void Start(); void Stop(); - void DrawMesh(const string& name, const Mesh& m); - void DrawPoints(const string& name, const Mesh& m, - FlatArray points); - void DrawLines(const string& name, const Mesh& m, - FlatArray lines); - void DrawTrigs(const string& name, const Mesh& m, - FlatArray trigs); - void DrawTets(const string& name, const Mesh& m, - FlatArray tets); - void AddComponent(const Mesh& m); + void DrawMesh(const Mesh& m); + void DrawPoints(FlatArray> position, string name = "points", + string color = "black"); + void DrawLines(FlatArray> position, string name = "lines", + string color = "black"); + void DrawTrigs(FlatArray> position, string name = "trigs", + string color = "black"); + + void DrawObject(FlatArray> position, string type, string name, + string color); + void ResetObjects() { data["objects"] = json::array_t(); } private: thread gui_thread; void *app, *loop, *token; std::set websockets; json data; + shared_ptr mesh; void Send(const json& data) { Send(data.dump()); } void Send(const string& message); diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 3cfb8bf4..d42b36b1 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -739,7 +739,7 @@ namespace netgen MeshQuality3d (mesh3d); } - debug_gui.DrawMesh("Mesh", mesh3d); + // debug_gui.DrawMesh("Mesh", mesh3d); multithread.task = savetask; return MESHING3_OK;