From d2f7c24a5e9cd0163d8995fb69aeea2d30bc1b92 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Fri, 14 Feb 2025 17:02:40 +0100 Subject: [PATCH] Utility function for debugging --- libsrc/meshing/debugging.cpp | 126 ++++++++++++++++++++++++++++++++++- libsrc/meshing/debugging.hpp | 9 ++- 2 files changed, 131 insertions(+), 4 deletions(-) diff --git a/libsrc/meshing/debugging.cpp b/libsrc/meshing/debugging.cpp index 43b837fd..feba17ab 100644 --- a/libsrc/meshing/debugging.cpp +++ b/libsrc/meshing/debugging.cpp @@ -1,8 +1,9 @@ #include +#include "debugging.hpp" namespace netgen { - unique_ptr GetOpenElements( const Mesh & m, int dom = 0, bool only_quads = false ) + unique_ptr GetOpenElements( const Mesh & m, int dom, bool only_quads ) { static Timer t("GetOpenElements"); RegionTimer rt(t); auto mesh = make_unique(); @@ -96,6 +97,127 @@ namespace netgen return mesh_ptr; } + void CheckMesh (const Mesh& mesh, MESHING_STEP step) + { + if (step == MESHCONST_OPTVOLUME) + { + bool have_error = false; + for (auto el : mesh.VolumeElements()) + { + double volume = el.Volume(mesh.Points()); + if (volume < 0) + { + have_error = true; + cout << "volume of element " << el << " is negative: " << volume << endl; + } + } + if (have_error) + throw Exception("Negative volume"); + CheckElementsAroundEdges(mesh); + } + } -} + void CheckElementsAroundEdges (const Mesh& mesh) + { + static Mesh last_good_mesh; + + Array> edges; + auto elementsonnode = mesh.CreatePoint2ElementTable(); + BuildEdgeList(mesh, elementsonnode, edges); + mesh.BoundaryEdge(1, 2); // trigger build of boundary edges + + ArrayMem hasbothpoints; + for (auto [pi0, pi1] : edges) + { + if (mesh.BoundaryEdge(pi0, pi1)) + continue; + + hasbothpoints.SetSize(0); + for (ElementIndex ei : elementsonnode[pi0]) + if (mesh[ei].PNums().Contains(pi1)) + hasbothpoints.Append(ei); + + bool skip = false; + for (ElementIndex ei : hasbothpoints) + { + if (mesh[ei].GetType() != TET) + { + skip = true; + break; + } + } + if (skip) + continue; + + int nsuround = hasbothpoints.Size(); + ArrayMem suroundpts(nsuround + 1); + suroundpts = PointIndex::INVALID; + ArrayMem tetused(nsuround); + tetused = false; + tetused[0] = true; + + auto el = mesh[hasbothpoints[0]]; + PointIndex pi2 = PointIndex::INVALID; + PointIndex pi3 = PointIndex::INVALID; + for (auto pi : el.PNums()) + if (pi != pi0 && pi != pi1) + { + pi3 = pi2; + pi2 = pi; + } + suroundpts[0] = pi2; + suroundpts[1] = pi3; + + for (auto i : Range(2, nsuround + 1)) + { + PointIndex oldpi = suroundpts[i - 1]; + PointIndex newpi = PointIndex::INVALID; + + for (int k = 0; k < nsuround && !newpi.IsValid(); k++) + if (!tetused[k]) + { + const Element& nel = mesh[hasbothpoints[k]]; + for (int k2 = 0; k2 < 4 && !newpi.IsValid(); k2++) + if (nel[k2] == oldpi) + { + newpi = nel[0] - pi0 + nel[1] - pi1 + nel[2] - oldpi + nel[3]; + tetused[k] = true; + suroundpts[i] = newpi; + + ArrayMem nelpts{nel[0], nel[1], nel[2], nel[3]}; + ArrayMem check_points{pi0, pi1, oldpi, newpi}; + QuickSort(check_points); + QuickSort(nelpts); + if (check_points != nelpts) + { + cout << __FILE__ << ":" << __LINE__ << "\tFound error" << endl; + cout << "i = " << i << endl; + cout << "oldpi = " << oldpi << endl; + cout << "newpi = " << newpi << endl; + cout << "Elements: " << endl; + cout << "nel " << nel << endl; + for (auto ei : hasbothpoints) + cout << mesh[ei] << endl; + cout << endl; + cout << "check_points: " << check_points << endl; + cout << "nelpts: " << nelpts << endl; + cout << "hasbothpoints: " << hasbothpoints << endl; + cout << "suroundpts: " << suroundpts << endl; + throw Exception("Found error"); + } + } + } + } + if (suroundpts.Last() != suroundpts[0]) + { + cout << __FILE__ << ":" << __LINE__ << "\tFound error" << endl; + cout << "hasbothpoints: " << hasbothpoints << endl; + cout << "suroundpts: " << suroundpts << endl; + for (auto ei : hasbothpoints) + cout << mesh[ei] << endl; + throw Exception("Found error"); + } + } + } +} // namespace netgen diff --git a/libsrc/meshing/debugging.hpp b/libsrc/meshing/debugging.hpp index b8c0ef7e..722fffa0 100644 --- a/libsrc/meshing/debugging.hpp +++ b/libsrc/meshing/debugging.hpp @@ -1,4 +1,4 @@ -#include "meshclass.hpp" +#include "meshfunc.hpp" namespace netgen @@ -7,6 +7,11 @@ namespace netgen unique_ptr FilterMesh( const Mesh & m, FlatArray points, FlatArray sels = Array{}, FlatArray els = Array{} ); + // Checks if the mesh is valid. This is called automatically on various places if debugparam.slowchecks is set + void CheckMesh( const Mesh & m, MESHING_STEP meshing_step ); - + // Sometimes during SwapImprove we discover topological errors in the mesh. For instance, an edge is adjacent to 8 tets around it, but + // the 8 "other" points of the tets don't form a closed path around the edge. Instead there are 2 sets of 4 points/tets each, which are not connected. + // This function checks for such errors and returns true if any are found. + void CheckElementsAroundEdges( const Mesh & m ); }