From 3fd1e82624e158fd7717c4bf584cc0a0ecd4d246 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:23:51 -0400 Subject: [PATCH 1/7] Migrate face/edge attributes off tets that remove_tets_by_ids deletes Attributes live only at their canonical slot (lowest incident tet id); deleting that tet re-canonicalizes surviving faces and edges onto cells the operations never maintained, so later reads see stale garbage. Migrate every such attribute to the simplex's new owner before touching the connectivity, mirroring the invariant consolidate_mesh maintains when renumbering. Latent and independent of the feature work; found while chasing a filter artifact that turned out to be the documented winding wrinkle (this was measured NOT to be that artifact's cause). Only the end-of-run filters call remove_tets_by_ids, and no suite config filters, so no reference output moves. Co-Authored-By: Claude Fable 5 --- src/wmtk/TetMesh.h | 90 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/src/wmtk/TetMesh.h b/src/wmtk/TetMesh.h index 3a5a69d199..9510d1d614 100644 --- a/src/wmtk/TetMesh.h +++ b/src/wmtk/TetMesh.h @@ -1196,6 +1196,96 @@ class TetMesh */ void remove_tets_by_ids(const std::vector& tids) { + // A face or edge attribute lives at ONE canonical slot, (lowest incident tet) * 4 + // or 6 + local index. Deleting a tet that is a simplex's canonical owner silently + // re-canonicalizes the simplex onto a surviving tet whose cell was never + // maintained -- the operations only ever write the canonical slot -- so the + // surviving cell holds stale garbage. Measured: after filtering, interior faces + // read back m_is_surface_fs = true 0.12 (eps 0.09) from the real surface, and the + // extracted "tracked surface" grows phantom triangles. Migrate every attribute + // whose owner dies to the simplex's new owner BEFORE touching the connectivity. + if (p_face_attrs != nullptr || p_edge_attrs != nullptr) { + std::vector removed(tet_capacity(), false); + for (const size_t tid : tids) { + removed[tid] = true; + } + for (const size_t tid : tids) { + const auto& tconn = m_tet_connectivity[tid]; + if (p_face_attrs != nullptr) { + for (int lf = 0; lf < 4; ++lf) { + std::array vids = { + {tconn[m_local_faces[lf][0]], + tconn[m_local_faces[lf][1]], + tconn[m_local_faces[lf][2]]}}; + // The face's other incident tet, if any. + auto common = set_intersection( + m_vertex_connectivity[vids[0]].m_conn_tets, + m_vertex_connectivity[vids[1]].m_conn_tets); + common = set_intersection( + common, + m_vertex_connectivity[vids[2]].m_conn_tets); + size_t other = std::numeric_limits::max(); + for (const size_t t : common) { + if (t != tid) { + other = t; + } + } + // Migrate only when this tet was the canonical owner (tid < other) + // and the other side survives; other < tid means the attribute + // already lives on the survivor. + if (other == std::numeric_limits::max() || removed[other] || + other < tid) { + continue; + } + std::array sorted = vids; + std::sort(sorted.begin(), sorted.end()); + for (int lo = 0; lo < 4; ++lo) { + std::array ovids = { + {m_tet_connectivity[other][m_local_faces[lo][0]], + m_tet_connectivity[other][m_local_faces[lo][1]], + m_tet_connectivity[other][m_local_faces[lo][2]]}}; + std::sort(ovids.begin(), ovids.end()); + if (ovids == sorted) { + p_face_attrs->move(tid * 4 + lf, other * 4 + lo); + break; + } + } + } + } + if (p_edge_attrs != nullptr) { + for (int le = 0; le < 6; ++le) { + const size_t v0 = tconn[m_local_edges[le][0]]; + const size_t v1 = tconn[m_local_edges[le][1]]; + const auto fan = set_intersection( + m_vertex_connectivity[v0].m_conn_tets, + m_vertex_connectivity[v1].m_conn_tets); + // Only the canonical owner migrates, once, to the lowest survivor. + if (fan.empty() || fan.front() != tid) { + continue; + } + size_t new_owner = std::numeric_limits::max(); + for (const size_t t : fan) { + if (!removed[t]) { + new_owner = t; // fan is sorted ascending; first kept wins + break; + } + } + if (new_owner == std::numeric_limits::max()) { + continue; // the edge dies with its whole fan + } + for (int lo = 0; lo < 6; ++lo) { + const size_t w0 = m_tet_connectivity[new_owner][m_local_edges[lo][0]]; + const size_t w1 = m_tet_connectivity[new_owner][m_local_edges[lo][1]]; + if ((w0 == v0 && w1 == v1) || (w0 == v1 && w1 == v0)) { + p_edge_attrs->move(tid * 6 + le, new_owner * 6 + lo); + break; + } + } + } + } + } + } + for (size_t tid : tids) { m_tet_connectivity[tid].m_is_removed = true; for (int j = 0; j < 4; j++) From ab00a94563a7501f5342761e1d45979e1df12d5f Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:24:01 -0400 Subject: [PATCH 2/7] 2D: anchor input free points, and keep feature outputs across filtering embed_segments returns its point provenance (computed and previously discarded); isolated input vertices and the new input_points files are anchored through the existing feature-point machinery -- same ball, same collapse policy, same retention audit. Free points are detected on the ORIGINAL input, so a sub-eps loop degenerating during simplification is not mistaken for a user point. Explicit features are collected before any filter deletes triangles: the mesh loses a discarded region, the new _features.obj (tracked curves + anchors, written only when free points were supplied) and the retention audit never lose a user-supplied feature. Measured: 12/12 retained across eps_rel and stop_energy sweeps (0/12 with preserve_feature_points off, worst 117x eps); the flood case with points in the discarded region reports 12/12 where it lost them before. Co-Authored-By: Claude Fable 5 --- .../wmtk/components/triwild/TriWildMesh.cpp | 29 ++- .../wmtk/components/triwild/TriWildMesh.h | 12 +- .../components/triwild/tests/CMakeLists.txt | 1 + .../triwild/tests/test_free_points.cpp | 179 ++++++++++++++++++ .../wmtk/components/triwild/triwild.cpp | 162 +++++++++++++++- .../wmtk/components/triwild/triwild_spec.json | 12 ++ src/wmtk/utils/EmbedSegments.cpp | 17 +- src/wmtk/utils/EmbedSegments.hpp | 9 +- 8 files changed, 409 insertions(+), 12 deletions(-) create mode 100644 components/triwild/wmtk/components/triwild/tests/test_free_points.cpp diff --git a/components/triwild/wmtk/components/triwild/TriWildMesh.cpp b/components/triwild/wmtk/components/triwild/TriWildMesh.cpp index 1495cc712e..c8c9649fe3 100644 --- a/components/triwild/wmtk/components/triwild/TriWildMesh.cpp +++ b/components/triwild/wmtk/components/triwild/TriWildMesh.cpp @@ -46,7 +46,8 @@ void TriWildMesh::init_mesh( const MatrixXi& E, const std::vector& tag_names, const MatrixXd& V_env, - const MatrixXi& E_env) + const MatrixXi& E_env, + const std::vector& free_point_vids) { assert(V.cols() == 2); assert(F.cols() == 3); @@ -224,6 +225,32 @@ void TriWildMesh::init_mesh( m_feature_points.size(), m_envelope_eps); } + + // Input free points: anchored unconditionally, junctions' cleanup exemption + // included. A junction is anchored (or not) by policy, because it is DERIVED from + // the curve network; a free point IS the input, so there is no policy question. A + // vertex that already carries a feature id (a free point coinciding with a polyline + // endpoint, merged by the arrangement's dedup) keeps the id it has -- one anchor at + // that position is enough, the retention audit is geometric. + size_t n_free = 0; + for (const size_t v : free_point_vids) { + if (v >= vert_capacity()) { + log_and_throw_error("Free-point vertex id {} out of range", v); + } + if (m_vertex_extra[v].m_feature_id != NO_FEATURE) { + continue; + } + m_vertex_extra[v].m_feature_id = m_feature_points.size(); + m_feature_points.push_back(m_vertex_attribute[v].m_posf); + ++n_free; + } + if (n_free > 0) { + logger().info( + "feature points: {} input free points anchored within {:.6} of their input " + "positions", + n_free, + m_envelope_eps); + } } // init envelope diff --git a/components/triwild/wmtk/components/triwild/TriWildMesh.h b/components/triwild/wmtk/components/triwild/TriWildMesh.h index 3d175dc054..a7d86844f5 100644 --- a/components/triwild/wmtk/components/triwild/TriWildMesh.h +++ b/components/triwild/wmtk/components/triwild/TriWildMesh.h @@ -172,6 +172,11 @@ class TriWildMesh : public wmtk::TriOptimizerMesh * stay near what the user gave us, not near the simplified version of it. Same * arrangement as tetwild, which hands its optimizer the envelope built on the * unsimplified input surface. + * @param free_point_vids rows of V that stand for input FREE POINTS -- input vertices + * with no incident segment. Each is anchored in m_feature_points exactly like a + * polyline endpoint: same ball, same collapse policy, same retention audit. They + * cannot be derived here the way endpoints and junctions are, because their + * valence in E is 0 -- the same valence as every background-grid vertex. */ void init_mesh( const MatrixXd& V, @@ -180,7 +185,8 @@ class TriWildMesh : public wmtk::TriOptimizerMesh const MatrixXi& E, const std::vector& tag_names, const MatrixXd& V_env, - const MatrixXi& E_env); + const MatrixXi& E_env, + const std::vector& free_point_vids = {}); void init_surfaces_and_boundaries(); @@ -235,8 +241,8 @@ class TriWildMesh : public wmtk::TriOptimizerMesh { if (const size_t n = m_feature_rejects.load(); n > 0) { logger().info( - "[feature] {} collapses refused to keep a polyline endpoint or junction " - "within {:.6} of its input position", + "[feature] {} collapses refused to keep a feature point (polyline endpoint, " + "junction, or input free point) within {:.6} of its input position", n, m_envelope_eps); } diff --git a/components/triwild/wmtk/components/triwild/tests/CMakeLists.txt b/components/triwild/wmtk/components/triwild/tests/CMakeLists.txt index 02835c18c8..a7e1ace433 100644 --- a/components/triwild/wmtk/components/triwild/tests/CMakeLists.txt +++ b/components/triwild/wmtk/components/triwild/tests/CMakeLists.txt @@ -12,6 +12,7 @@ set(SRC_FILES test_simplify_segments.cpp test_winding_tags.cpp test_feature_points.cpp + test_free_points.cpp test_coarsen_pass.cpp ) diff --git a/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp b/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp new file mode 100644 index 0000000000..cfcca97aa4 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp @@ -0,0 +1,179 @@ +#include +#include +#include + +#include + +#include +#include + +using namespace wmtk; +using namespace wmtk::components::triwild; + +// A free point is an input vertex with no incident segment. The arrangement triangulates +// every point it is handed, so a free point becomes an ordinary mesh vertex; what makes it +// a FEATURE is the anchor registered at init, which reuses the endpoint/junction machinery +// (test_feature_points.cpp) unchanged. These tests cover the two new pieces: the point map +// out of embed_segments, and the registration in init_mesh. + +TEST_CASE("embed-segments-point-map", "[triwild_operation][free_points]") +{ + // One segment plus two free points: one strictly inside the future triangulation, one + // exactly ON the segment's interior -- the point the arrangement must MERGE with a + // constrained-edge vertex rather than keep separate. + MatrixXd V(4, 2); + V << 0, 0, // + 2, 0, // segment (0,0)-(2,0) + 1, 1, // free point off the segment + 1, 0; // free point ON the segment interior + MatrixXi E(1, 2); + E << 0, 1; + + MatrixXd V_out; + std::vector V_rational; + MatrixXi F_out, E_out; + std::vector point_map; + utils::embed_segments(V, E, V_out, V_rational, F_out, E_out, nullptr, &point_map); + + REQUIRE(point_map.size() == 4); + for (int i = 0; i < 4; ++i) { + REQUIRE(point_map[i] >= 0); + REQUIRE(point_map[i] < V_out.rows()); + // Every input point survives at its exact coordinates. + CHECK(V_out(point_map[i], 0) == V(i, 0)); + CHECK(V_out(point_map[i], 1) == V(i, 1)); + } + + // The on-segment point splits the constrained edge: its output vertex must be an + // endpoint of some constrained edge, and the segment must now be tiled by two of them. + const int on_seg = point_map[3]; + int incident = 0; + for (int e = 0; e < E_out.rows(); ++e) { + if (E_out(e, 0) == on_seg || E_out(e, 1) == on_seg) { + ++incident; + } + } + CHECK(incident == 2); + CHECK(E_out.rows() == 2); + + // The off-segment point is on no constrained edge. + const int off_seg = point_map[2]; + for (int e = 0; e < E_out.rows(); ++e) { + CHECK(E_out(e, 0) != off_seg); + CHECK(E_out(e, 1) != off_seg); + } +} + +TEST_CASE("embed-segments-point-map-duplicates", "[triwild_operation][free_points]") +{ + // Exact duplicates merge to the same output vertex; the map reports it for both. + MatrixXd V(4, 2); + V << 0, 0, // + 1, 0, // + 0.5, 0.5, // + 0.5, 0.5; + MatrixXi E(1, 2); + E << 0, 1; + + MatrixXd V_out; + std::vector V_rational; + MatrixXi F_out, E_out; + std::vector point_map; + utils::embed_segments(V, E, V_out, V_rational, F_out, E_out, nullptr, &point_map); + + REQUIRE(point_map.size() == 4); + CHECK(point_map[2] >= 0); + CHECK(point_map[2] == point_map[3]); +} + +TEST_CASE("free-point-coincident-with-endpoint", "[triwild_operation][free_points]") +{ + // A free point EXACTLY at an open polyline's endpoint. The arrangement merges the two + // input points into one vertex; at init that vertex is first anchored as an endpoint, + // and the free-point registration must then KEEP that id rather than stack a second + // anchor at the same position -- the retention audit is geometric, one anchor covers + // both readings of the point. + MatrixXd V(3, 2); + V << 0, 0, // + 2, 0, // open polyline (0,0)-(2,0) + 0, 0; // free point, exactly the first endpoint + MatrixXi E(1, 2); + E << 0, 1; + + MatrixXd V_arr; + std::vector V_rational; + MatrixXi F_arr, E_arr; + std::vector point_map; + utils::embed_segments(V, E, V_arr, V_rational, F_arr, E_arr, nullptr, &point_map); + + // The duplicate merged: the free point maps to the same vertex as the endpoint. + REQUIRE(point_map[2] >= 0); + CHECK(point_map[2] == point_map[0]); + + Parameters params; + params.init(Vector2d(0, 0), Vector2d(2, 0)); + TriWildMesh mesh(params, /*envelope_eps=*/0.1); + mesh.init_mesh( + V_arr, + V_rational, + F_arr, + E_arr, + /*tag_names=*/{}, + V, + E, + /*free_point_vids=*/{size_t(point_map[2])}); + + // Exactly the polyline's two endpoint anchors -- no third anchor for the free point. + REQUIRE(mesh.m_feature_points.size() == 2); + // The shared vertex carries the endpoint's id (the first registered), and both anchors + // are covered geometrically. + CHECK(mesh.m_vertex_extra[size_t(point_map[0])].m_feature_id != NO_FEATURE); + const auto [kept, total] = mesh.feature_retention(); + CHECK(kept == 2); + CHECK(total == 2); +} + +TEST_CASE("init-mesh-registers-free-points", "[triwild_operation][free_points]") +{ + // Run the real pipeline front end -- arrangement, then init_mesh -- on a square curve + // with one free point inside, and check the anchor comes out the other side. + MatrixXd V(5, 2); + V << 0, 0, // + 4, 0, // + 4, 4, // + 0, 4, // + 2, 2; // free point + MatrixXi E(4, 2); + E << 0, 1, 1, 2, 2, 3, 3, 0; + + MatrixXd V_arr; + std::vector V_rational; + MatrixXi F_arr, E_arr; + std::vector point_map; + utils::embed_segments(V, E, V_arr, V_rational, F_arr, E_arr, nullptr, &point_map); + REQUIRE(point_map[4] >= 0); + + Parameters params; + params.init(Vector2d(0, 0), Vector2d(4, 4)); + TriWildMesh mesh(params, /*envelope_eps=*/0.1); + mesh.init_mesh( + V_arr, + V_rational, + F_arr, + E_arr, + /*tag_names=*/{}, + V, + E, + /*free_point_vids=*/{size_t(point_map[4])}); + + // The closed square has no endpoints and no junctions, so the free point is the only + // anchor. + REQUIRE(mesh.m_feature_points.size() == 1); + CHECK(mesh.m_feature_points[0] == Vector2d(2, 2)); + CHECK(mesh.m_vertex_extra[size_t(point_map[4])].m_feature_id == 0); + + // And the retention audit sees it. + const auto [kept, total] = mesh.feature_retention(); + CHECK(kept == 1); + CHECK(total == 1); +} diff --git a/components/triwild/wmtk/components/triwild/triwild.cpp b/components/triwild/wmtk/components/triwild/triwild.cpp index 974e335c39..65f0813fd4 100644 --- a/components/triwild/wmtk/components/triwild/triwild.cpp +++ b/components/triwild/wmtk/components/triwild/triwild.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -258,6 +259,55 @@ void triwild(nlohmann::json json_params) params.preserve_topology ? 0.0 : double(json_params["remove_duplicate_eps"]); wmtk::utils::read_input_curves(input_paths, remove_duplicate_eps, V_in, E_in, Vs, Es); + // Input point files: every vertex is a feature point, appended to V_in as an isolated + // vertex (referenced by no segment). Appended BEFORE the bounding box is taken, so a + // point outside the curves' box still ends up inside the triangulated domain. Not added + // to Vs/Es: points carry no winding number. + { + std::vector point_paths = json_params["input_points"]; + for (std::string& p : point_paths) { + p = resolve_path(root, p).string(); + } + for (const std::string& path : point_paths) { + MatrixXd Vp; + MatrixXi Ep; + io::read_edge_mesh(path, Vp, Ep, remove_duplicate_eps); + if (Ep.rows() > 0) { + logger().warn( + "input_points file {} has {} edges; only its {} vertices are used", + path, + Ep.rows(), + Vp.rows()); + } + logger().info("Read point file {}: #P = {}", path, Vp.rows()); + const int base = V_in.rows(); + V_in.conservativeResize(base + Vp.rows(), 2); + V_in.block(base, 0, Vp.rows(), 2) = Vp.block(0, 0, Vp.rows(), 2); + } + } + + // Free points to preserve: every vertex of V_in with no incident segment -- the point + // files above, plus any isolated vertex the curve files contained (read_edge_mesh keeps + // them for exactly this). Detected on the ORIGINAL input, not the simplified network: a + // sub-eps loop can degenerate to an isolated vertex during simplification, and such a + // remnant is the simplification doing its job, not an input point to pin. + std::vector free_points; + { + std::vector valence(V_in.rows(), 0); + for (int i = 0; i < E_in.rows(); ++i) { + ++valence[E_in(i, 0)]; + ++valence[E_in(i, 1)]; + } + for (int v = 0; v < V_in.rows(); ++v) { + if (valence[v] == 0) { + free_points.emplace_back(V_in.row(v)); + } + } + if (!free_points.empty()) { + logger().info("input free points: {}", free_points.size()); + } + } + // Informational input-topology report; gated behind DEBUG_euler because it is only // meaningful next to the matching computations later in the run. const bool debug_euler = json_params["DEBUG_euler"]; @@ -358,7 +408,47 @@ void triwild(nlohmann::json json_params) std::vector V_rational; // the same vertices, exact MatrixXi F; MatrixXi E; // constraint edges in the arrangement - wmtk::utils::embed_segments(V_simp, E_simp, V, V_rational, F, E); + std::vector point_map; + wmtk::utils::embed_segments( + V_simp, + E_simp, + V, + V_rational, + F, + E, + nullptr, + free_points.empty() ? nullptr : &point_map); + + // Find each free point's vertex in the arrangement. By POSITION, not by index: the + // simplification compacts V, so input row ids do not survive it -- but an isolated + // vertex has no incident segment, so no collapse ever moves it and its coordinates in + // V_simp are bit-identical to the input's. The arrangement's point provenance then maps + // that row to the output vertex. + std::vector free_point_vids; + if (!free_points.empty()) { + std::map, int> simp_row_of; + for (int v = 0; v < V_simp.rows(); ++v) { + simp_row_of.emplace(std::make_pair(V_simp(v, 0), V_simp(v, 1)), v); + } + for (const Vector2d& p : free_points) { + const auto it = simp_row_of.find({p[0], p[1]}); + if (it == simp_row_of.end()) { + log_and_throw_error( + "Input free point ({}, {}) not found after simplification; isolated " + "vertices must survive it unmoved", + p[0], + p[1]); + } + const int vid = point_map[it->second]; + if (vid < 0) { + log_and_throw_error( + "Input free point ({}, {}) was dropped by the arrangement", + p[0], + p[1]); + } + free_point_vids.push_back(size_t(vid)); + } + } // The arrangement is the baseline for everything after it. It is EXPECTED to differ // from the input: resolving a crossing inserts a vertex shared by both curves, which @@ -440,7 +530,7 @@ void triwild(nlohmann::json json_params) TriWildMesh mesh(params, opt_eps, NUM_THREADS); wmtk::set_preallocation_factor_from_json(mesh, json_params); - mesh.init_mesh(V, V_rational, F, E, tag_names, V_env, E_env); + mesh.init_mesh(V, V_rational, F, E, tag_names, V_env, E_env, free_point_vids); // After init_mesh, which is what builds the envelope, and after the simplification, which // uses its own object -- so this only disables the checks the optimizer makes. @@ -493,6 +583,31 @@ void triwild(nlohmann::json json_params) "output groups will be empty."); } + // Feature collections, taken BEFORE any filter deletes triangles: user-supplied + // features must not vanish from the feature outputs and audits when the region they + // live in is discarded (they do leave the triangle mesh itself). Same rule as 3D. + std::vector> collected_curve_segs; + std::vector collected_anchor_pts; + std::pair collected_retention{0, 0}; + double collected_retention_worst = 0; + { + for (const auto& e : mesh.get_edges()) { + if (!mesh.m_edge_attribute[e.eid(mesh)].m_is_surface_fs) { + continue; + } + collected_curve_segs.push_back( + {{mesh.m_vertex_attribute[e.vid(mesh)].m_posf, + mesh.m_vertex_attribute[e.switch_vertex(mesh).vid(mesh)].m_posf}}); + } + for (const auto& v : mesh.get_vertices()) { + const size_t vid = v.vid(mesh); + if (mesh.m_vertex_extra[vid].m_feature_id != NO_FEATURE) { + collected_anchor_pts.push_back(mesh.m_vertex_attribute[vid].m_posf); + } + } + collected_retention = mesh.feature_retention(&collected_retention_worst); + } + if (filter_option == "input") { mesh.filter_with_input_winding_number(); mesh.consolidate_mesh(); @@ -608,16 +723,20 @@ void triwild(nlohmann::json json_params) double feat_worst_ratio = 0; size_t feat_kept = 0, feat_total = 0; if (json_params["DEBUG_feature_retention"]) { - std::tie(feat_kept, feat_total) = mesh.feature_retention(&feat_worst_ratio); + // Pre-filter numbers (see the collection above): a feature in a discarded region + // left the triangle mesh with its region, but was preserved up to extraction and + // survives in _features.obj. + std::tie(feat_kept, feat_total) = collected_retention; + feat_worst_ratio = collected_retention_worst; } if (feat_total > 0) { if (feat_kept == feat_total) { logger().info("feature points retained: {}/{}", feat_kept, feat_total); } else { logger().warn( - "feature points retained: {}/{} -- {} polyline endpoints or junctions are no " - "longer represented within eps; the worst is {:.2f} x eps from the nearest " - "vertex", + "feature points retained: {}/{} -- {} feature points (polyline endpoints, " + "junctions, or input free points) are no longer represented within eps; the " + "worst is {:.2f} x eps from the nearest vertex", feat_kept, feat_total, feat_total - feat_kept, @@ -714,6 +833,37 @@ void triwild(nlohmann::json json_params) } mesh.write_msh_groups(output_path + ".msh"); + // The tracked curves and feature anchors, as an edge mesh with point records -- the 2D + // counterpart of tetwild's _features.obj, from the same pre-filter collections. Only + // written when the input supplied free points, mirroring 3D's features-present gate -- + // and keeping featureless runs byte-identical, output file set included. + if (!free_points.empty()) { + std::ofstream fout(output_path + "_features.obj"); + std::map, size_t> vid_of; + const auto obj_vertex = [&](const Vector2d& p) { + const std::array key = {{p[0], p[1]}}; + const auto [it, inserted] = vid_of.emplace(key, vid_of.size() + 1); + if (inserted) { + fout << "v " << p[0] << " " << p[1] << " 0\n"; + } + return it->second; + }; + std::vector> obj_edges; + for (const auto& seg : collected_curve_segs) { + obj_edges.push_back({{obj_vertex(seg[0]), obj_vertex(seg[1])}}); + } + std::vector obj_points; + for (const Vector2d& p : collected_anchor_pts) { + obj_points.push_back(obj_vertex(p)); + } + for (const auto& e : obj_edges) { + fout << "l " << e[0] << " " << e[1] << "\n"; + } + for (const size_t pid : obj_points) { + fout << "p " << pid << "\n"; + } + } + logger().info("======= finish ========="); } diff --git a/components/triwild/wmtk/components/triwild/triwild_spec.json b/components/triwild/wmtk/components/triwild/triwild_spec.json index 6e9be5223f..afc84fc857 100644 --- a/components/triwild/wmtk/components/triwild/triwild_spec.json +++ b/components/triwild/wmtk/components/triwild/triwild_spec.json @@ -6,6 +6,7 @@ "required": ["application", "input"], "optional": [ "output", + "input_points", "input_names", "input_dir", "num_threads", @@ -81,6 +82,17 @@ "type": "string", "doc": "Triangular input mesh." }, + { + "pointer": "/input_points", + "type": "list", + "default": [], + "doc": "List of input POINT files: every vertex of each file is a feature point to preserve, whether or not the file also has edges (they are ignored). Each point is inserted as a vertex of the triangulation and anchored like a polyline endpoint: some output vertex stays within eps of it, collapses that would break that are refused, and smoothing may move the carrying vertex only within the eps ball. Free points can also be supplied inline: a vertex of an /input curve file with no incident edge is anchored the same way. A point lying in a region the /filter discards leaves the triangle mesh with that region, but survives in the _features.obj feature output, and DEBUG_feature_retention reports the PRE-filter state -- a feature preserved through the optimization is never reported lost because extraction discarded its region." + }, + { + "pointer": "/input_points/*", + "type": "string", + "doc": "One input point file (any edge-mesh format; only the vertices are used)." + }, { "pointer": "/output", "type": "string", diff --git a/src/wmtk/utils/EmbedSegments.cpp b/src/wmtk/utils/EmbedSegments.cpp index 242f4c4e42..fb2988d0e8 100644 --- a/src/wmtk/utils/EmbedSegments.cpp +++ b/src/wmtk/utils/EmbedSegments.cpp @@ -115,7 +115,8 @@ void embed_segments( std::vector& V_rational, MatrixXi& F_out, MatrixXi& E_out, - std::vector>* E_out_sources) + std::vector>* E_out_sources, + std::vector* point_map) { assert(V.cols() == 2); assert(E.cols() == 2); @@ -227,6 +228,20 @@ void embed_segments( } } + // Where did each input point end up? The remesher reports {triangle, vertex} per input + // point, in input order; the input points went in first, so the first V.rows() entries + // are ours and the background-grid entries after them are dropped. + if (point_map != nullptr) { + assert(point_provenance.size() >= size_t(V.rows())); + point_map->assign(V.rows(), -1); + for (int i = 0; i < V.rows(); ++i) { + const uint32_t v = point_provenance[i][1]; + if (v != UINT32_MAX) { + (*point_map)[i] = int(v); + } + } + } + logger().info( "2D arrangement: #V = {}, #F = {}, #E_constrained = {} ({} vertices have no exact " "double representation)", diff --git a/src/wmtk/utils/EmbedSegments.hpp b/src/wmtk/utils/EmbedSegments.hpp index 460bc5d9ab..8d2b43e8b9 100644 --- a/src/wmtk/utils/EmbedSegments.hpp +++ b/src/wmtk/utils/EmbedSegments.hpp @@ -31,6 +31,12 @@ namespace wmtk::utils { * @param[out] E_out_sources optional: for each row of E_out, the input edges (rows of E) it * tiles, ascending. Usually one; more where input edges overlap, which is * exactly the case a geometric look-up on E_out alone cannot tell apart. + * @param[out] point_map optional: for each row of V, the row of V_out holding that exact + * point. The arrangement triangulates every input point, referenced by a + * segment or not, so this is how a free point (no incident segment) is found + * again in the output; exact duplicates in V map to the same output row. -1 + * only if the remesher dropped the point, which it does not do for points -- + * the entry exists so a caller can assert that rather than assume it. */ void embed_segments( const MatrixXd& V, @@ -39,7 +45,8 @@ void embed_segments( std::vector& V_rational, MatrixXi& F_out, MatrixXi& E_out, - std::vector>* E_out_sources = nullptr); + std::vector>* E_out_sources = nullptr, + std::vector* point_map = nullptr); /** * @brief Read every input edge mesh and concatenate them into one segment network. From 728b4936115a03d04ceab83ab9209ccbd30725fd Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:24:20 -0400 Subject: [PATCH 3/7] 3D core: carry feature-edge tags through the operations, guarded by a tube A per-edge union bool (m_is_feature_edge), mirroring how the tracked surface is one m_is_surface_fs for every input. An edge attribute's slot is (lowest incident tet)*6+local, so the tags are cached by vertex pair around every operation and rewritten after: region-driven over the created tets for split and the swaps (which also scrubs stale values from reused slots), and cache-driven with v1 renamed to v2 for collapse -- whose deleted tets' link edges contain neither collapse endpoint, a case a tet FACE cannot exhibit, and can live on in tets outside any enumerable region. Swaps that would delete a tagged edge (3-2, 4-4, 5-6) are refused. Geometric guards: m_feature_envelope, the tube a tagged edge must stay in. The collapse guard checks every affected tagged edge post-collapse against it -- curves coarsen along themselves and cannot be dragged sideways -- and smoothing pulls feature vertices (interior ones included, which the surface gate used to smooth as free) toward the tube and hard-checks their tagged edges. A smoothing_position_is_allowed hook joins the 3D path (the 2D one existed), for the anchor balls of 0-dimensional features. Under preserve_topology the substructure link condition consults overridable orders (substructure_order_of_*), so a tagged edge counts as 1D substructure -- including interior feature edges the surface-face walk cannot see -- without changing what get_order_of_* means anywhere else. Everything is behind m_track_feature_edges -- never behind p_edge_attrs, which topological_offset registers with its own attribute type. MshData gains 0D point elements (gmsh type 15) for the mixed output. Co-Authored-By: Claude Fable 5 --- src/wmtk/TetMesh.h | 29 ++++- src/wmtk/TetMeshSubstructure.cpp | 38 ++++-- src/wmtk/TetOptimizerMesh.cpp | 170 +++++++++++++++++++++++++ src/wmtk/TetOptimizerMesh.h | 119 +++++++++++++++++ src/wmtk/TetOptimizerMeshCollapse.cpp | 45 +++++++ src/wmtk/TetOptimizerMeshSplit.cpp | 25 ++++ src/wmtk/TetOptimizerMeshSwaps.cpp | 44 +++++++ src/wmtk/optimization/SmoothVertex.hpp | 21 ++- src/wmtk/utils/io.hpp | 24 +++- 9 files changed, 498 insertions(+), 17 deletions(-) diff --git a/src/wmtk/TetMesh.h b/src/wmtk/TetMesh.h index 9510d1d614..bd5c8c2d2a 100644 --- a/src/wmtk/TetMesh.h +++ b/src/wmtk/TetMesh.h @@ -22,14 +22,17 @@ namespace wmtk { class TetMesh { -private: +protected: /** * @brief local edges within a tet * + * Protected rather than private: TetOptimizerMesh's feature-edge tracker enumerates a + * tet's six edges by local id. */ static constexpr std::array, 6> m_local_edges = { {{{0, 1}}, {{1, 2}}, {{0, 2}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}}; +private: static constexpr std::array m_map_vertex2edge = {{0, 0, 1, 3}}; static constexpr std::array m_map_vertex2oppo_face = {{3, 1, 2, 0}}; static constexpr std::array m_map_edge2face = {{0, 0, 0, 1, 2, 1}}; @@ -1623,6 +1626,30 @@ class TetMesh */ size_t get_order_of_edge(const std::array& vids) const; + /** + * @brief The orders substructure_link_condition consults, overridable so an application + * can widen the 1D substructure beyond what the surface complex derives. + * + * The paper's condition is written in terms of simplex orders; TetOptimizerMesh overrides + * these so an input FEATURE edge counts as order 2 (and its endpoints accordingly) without + * touching the meaning of get_order_of_* anywhere else -- the feature tags are + * deliberately not part of the order classification. + */ + virtual size_t substructure_order_of_edge(const std::array& vids) const + { + return get_order_of_edge(vids); + } + virtual size_t substructure_order_of_vertex(const size_t vid) const + { + return get_order_of_vertex(vid); + } + /** + * @brief One-ring vertices x of `vid` where (vid, x) is 1D substructure NOT derivable + * from the surface complex (feature edges through the interior, for instance). The link + * condition's order-2 edge collection walks surface faces and would never see them. + */ + virtual void substructure_feature_neighbors(size_t, std::vector&) const {} + /** * @brief Link condition that also considers substructures. * diff --git a/src/wmtk/TetMeshSubstructure.cpp b/src/wmtk/TetMeshSubstructure.cpp index bf20ffcbd7..6dc4d58546 100644 --- a/src/wmtk/TetMeshSubstructure.cpp +++ b/src/wmtk/TetMeshSubstructure.cpp @@ -222,9 +222,9 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const using namespace simplex; - const size_t edge_order = get_order_of_edge({{u_id, v_id}}); - const size_t u_order = get_order_of_vertex(u_id); - const size_t v_order = get_order_of_vertex(v_id); + const size_t edge_order = substructure_order_of_edge({{u_id, v_id}}); + const size_t u_order = substructure_order_of_vertex(u_id); + const size_t v_order = substructure_order_of_vertex(v_id); // If the edge is lower order than both vertices, we know for sure that this edge must not // be collapsed. Example: edge in space (order 0) connecting two surfaces (order 1). @@ -283,19 +283,28 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const auto& [ev0, ev1] = e_opp.vertices(); // collect order 2 edges - if (u_order > 1 || get_order_of_vertex(ev0) > 1) { + if (u_order > 1 || substructure_order_of_vertex(ev0) > 1) { const Edge e0(u_id, ev0); - if (get_order_of_edge(e0.vertices()) > 1) { + if (substructure_order_of_edge(e0.vertices()) > 1) { order2_edges.add(e0); } } - if (u_order > 1 || get_order_of_vertex(ev1) > 1) { + if (u_order > 1 || substructure_order_of_vertex(ev1) > 1) { const Edge e1(u_id, ev1); - if (get_order_of_edge(e1.vertices()) > 1) { + if (substructure_order_of_edge(e1.vertices()) > 1) { order2_edges.add(e1); } } } + // Feature edges are 1D substructure whether or not any surface face touches them; + // the loop above walks surface faces and cannot see an interior feature edge. + { + std::vector feat; + substructure_feature_neighbors(u_id, feat); + for (const size_t x : feat) { + order2_edges.add(Edge(u_id, x)); + } + } order2_edges.sort_and_clean(); for (const Edge& e : order2_edges.edges()) { const Face fw(e, w_id); @@ -338,19 +347,26 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const auto& [ev0, ev1] = e_opp.vertices(); // collect order 2 edges - if (v_order > 1 || get_order_of_vertex(ev0) > 1) { + if (v_order > 1 || substructure_order_of_vertex(ev0) > 1) { const Edge e0(v_id, ev0); - if (get_order_of_edge(e0.vertices()) > 1) { + if (substructure_order_of_edge(e0.vertices()) > 1) { order2_edges.add(e0); } } - if (v_order > 1 || get_order_of_vertex(ev1) > 1) { + if (v_order > 1 || substructure_order_of_vertex(ev1) > 1) { const Edge e1(v_id, ev1); - if (get_order_of_edge(e1.vertices()) > 1) { + if (substructure_order_of_edge(e1.vertices()) > 1) { order2_edges.add(e1); } } } + { + std::vector feat; + substructure_feature_neighbors(v_id, feat); + for (const size_t x : feat) { + order2_edges.add(Edge(v_id, x)); + } + } order2_edges.sort_and_clean(); for (const Edge& e : order2_edges.edges()) { const Face fw(e, w_id); diff --git a/src/wmtk/TetOptimizerMesh.cpp b/src/wmtk/TetOptimizerMesh.cpp index db870e3c93..e483535974 100644 --- a/src/wmtk/TetOptimizerMesh.cpp +++ b/src/wmtk/TetOptimizerMesh.cpp @@ -692,4 +692,174 @@ double TetOptimizerMesh::get_length2(const Tuple& l) const return length; } +// ---- Feature-edge tag propagation helpers (see the header banner) ------------------------ + +bool TetOptimizerMesh::vertex_has_feature_edge(const size_t vid) const +{ + for (const size_t u : get_one_ring_vids_for_vertex(vid)) { + const Tuple e = tuple_from_edge({{vid, u}}); + if (e.is_valid(*this) && m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { + return true; + } + } + return false; +} + +size_t TetOptimizerMesh::substructure_order_of_edge(const std::array& vids) const +{ + size_t order = TetMesh::substructure_order_of_edge(vids); + if (order < 2 && m_track_feature_edges) { + const Tuple e = tuple_from_edge(vids); + if (e.is_valid(*this) && m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { + order = 2; + } + } + return order; +} + +size_t TetOptimizerMesh::substructure_order_of_vertex(const size_t vid) const +{ + size_t order = TetMesh::substructure_order_of_vertex(vid); + if (order < 3 && m_track_feature_edges) { + size_t n_tagged = 0; + for (const size_t u : get_one_ring_vids_for_vertex(vid)) { + const Tuple e = tuple_from_edge({{vid, u}}); + if (e.is_valid(*this) && m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { + ++n_tagged; + } + } + if (n_tagged == 2) { + order = std::max(order, 2); + } else if (n_tagged == 1 || n_tagged >= 3) { + order = 3; // curve endpoint or junction + } + } + return order; +} + +void TetOptimizerMesh::substructure_feature_neighbors(const size_t vid, std::vector& out) + const +{ + if (!m_track_feature_edges) { + return; + } + for (const size_t u : get_one_ring_vids_for_vertex(vid)) { + const Tuple e = tuple_from_edge({{vid, u}}); + if (e.is_valid(*this) && m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { + out.push_back(u); + } + } +} + +bool TetOptimizerMesh::feature_edges_at_vertex_inside(const size_t vid) const +{ + if (!m_track_feature_edges || !m_feature_envelope) { + return true; + } + for (const size_t u : get_one_ring_vids_for_vertex(vid)) { + const Tuple e = tuple_from_edge({{vid, u}}); + if (!e.is_valid(*this) || !m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { + continue; + } + if (m_feature_envelope->is_outside(std::array{ + {m_vertex_attribute[vid].m_posf, m_vertex_attribute[u].m_posf}})) { + return false; + } + } + return true; +} + + +void TetOptimizerMesh::feature_edges_cache( + const std::vector& tids, + std::map, bool>& cache) +{ + cache.clear(); + for (const size_t tid : tids) { + const auto vs = oriented_tet_vids(tid); + for (int local_eid = 0; local_eid < 6; ++local_eid) { + const auto [l0, l1] = m_local_edges[local_eid]; + std::array pair = {{vs[l0], vs[l1]}}; + if (pair[0] > pair[1]) { + std::swap(pair[0], pair[1]); + } + if (cache.count(pair) != 0) { + continue; + } + const size_t eid = tuple_from_edge(tid, local_eid).eid(*this); + cache.emplace(pair, m_feature_edge_attribute[eid].m_is_feature_edge); + } + } +} + +void TetOptimizerMesh::feature_edges_restore_region( + const std::vector& tids, + const std::map, bool>& cache) +{ + std::vector tets; + tets.reserve(tids.size()); + for (const size_t tid : tids) { + tets.push_back(tuple_from_tet(tid)); + } + feature_edges_restore_region(tets, cache); +} + +void TetOptimizerMesh::feature_edges_restore_region( + const std::vector& tets, + const std::map, bool>& cache) +{ + for (const Tuple& t : tets) { + const size_t tid = t.tid(*this); + const auto vs = oriented_tet_vids(tid); + for (int local_eid = 0; local_eid < 6; ++local_eid) { + const auto [l0, l1] = m_local_edges[local_eid]; + std::array pair = {{vs[l0], vs[l1]}}; + if (pair[0] > pair[1]) { + std::swap(pair[0], pair[1]); + } + const auto it = cache.find(pair); + // A pair not in the cache is a genuinely new edge: default, never inherit -- + // the slot may be reused and hold a stale value. + const bool tag = it != cache.end() && it->second; + const size_t eid = tuple_from_edge(tid, local_eid).eid(*this); + m_feature_edge_attribute[eid].m_is_feature_edge = tag; + } + } +} + +bool TetOptimizerMesh::feature_edges_restore_remap( + const std::map, bool>& cache, + const size_t v_old, + const size_t v_new) +{ + // Remap and OR-merge first: (v_old,x) and (v_new,x) become the same pair. + std::map, bool> merged; + for (const auto& [pair, tag] : cache) { + std::array p = pair; + for (size_t& v : p) { + if (v == v_old) { + v = v_new; + } + } + if (p[0] == p[1]) { + continue; // the collapsed edge itself; its tag dies with it (the chain shortens) + } + if (p[0] > p[1]) { + std::swap(p[0], p[1]); + } + merged[p] = merged[p] || tag; + } + for (const auto& [pair, tag] : merged) { + const Tuple t = tuple_from_edge(pair); + if (!t.is_valid(*this)) { + if (tag) { + return false; // a tagged edge vanished; abort, as the face path does + } + continue; + } + m_feature_edge_attribute[t.eid(*this)].m_is_feature_edge = tag; + } + return true; +} + } // namespace wmtk diff --git a/src/wmtk/TetOptimizerMesh.h b/src/wmtk/TetOptimizerMesh.h index 33d44960fc..3805198d9f 100644 --- a/src/wmtk/TetOptimizerMesh.h +++ b/src/wmtk/TetOptimizerMesh.h @@ -106,6 +106,76 @@ class TetOptimizerMesh : public wmtk::TetMesh, public wmtk::RationalPositions */ AttributeContainerGroup m_face_attr_group; + /** + * @brief Per-edge feature tag: does this tet edge tile an input feature curve. + * + * A single union bool, deliberately, mirroring how the tracked surface is one + * m_is_surface_fs for every input file and sheet. Which curve an edge belongs to is not + * carried through the optimization anywhere in this codebase; identity questions are + * answered at the boundaries -- provenance at insertion, combinatorial and geometric + * audits at the end. + */ + struct FeatureEdgeAttributes + { + bool m_is_feature_edge = false; + }; + using EdgeAttCol = AttributeCollection; + /// Sized, protected and rolled back only once enable_feature_edge_tracking() ran. + EdgeAttCol m_feature_edge_attribute; + + /** + * @brief Whether the shared operations maintain m_feature_edge_attribute. + * + * The operations are guarded by THIS flag, never by `p_edge_attrs != nullptr`: + * topological_offset also registers an edge collection there, with its own attribute + * type, and reinterpreting it as feature tags would corrupt it. + */ + bool m_track_feature_edges = false; + + /** + * @brief Envelope around the input feature curves: the tube a tagged edge must stay in. + * + * The collapse guard's veto (see collapse_edge_before). Kept SEPARATE from the surface + * envelope and the order-2 envelope: merging would let feature curves vouch for open + * boundaries (and vice versa) in guard decisions. Null => tagged edges are propagated + * but not geometrically constrained. + */ + std::shared_ptr m_feature_envelope; + + /// Turn on feature-edge tracking. Must run before init() so the connectivity init sizes + /// the collection; refuses to share p_edge_attrs with an application that already uses it. + void enable_feature_edge_tracking() + { + if (p_edge_attrs != nullptr && p_edge_attrs != &m_feature_edge_attribute) { + log_and_throw_error( + "enable_feature_edge_tracking: p_edge_attrs is already registered by the " + "application; feature tracking cannot share it"); + } + p_edge_attrs = &m_feature_edge_attribute; + m_track_feature_edges = true; + } + + /// Does any edge incident to `vid` carry the feature tag? Derived from the edge tags + /// rather than from a vertex flag, so it cannot go stale. + bool vertex_has_feature_edge(size_t vid) const; + + // Substructure widening for the link condition (preserve_topology): a tagged feature + // edge is order-2 substructure, a vertex with two incident tagged edges is on a curve + // (order 2), with one or three-plus it is a curve endpoint or junction (order 3). The + // meaning of get_order_of_* is untouched. + size_t substructure_order_of_edge(const std::array& vids) const override; + size_t substructure_order_of_vertex(size_t vid) const override; + void substructure_feature_neighbors(size_t vid, std::vector& out) const override; + /// Are all tagged edges at `vid` inside the feature tube, at the CURRENT positions? + /// True when tracking or the tube is off. Smoothing calls this with the candidate + /// position already written into the vertex attribute. + bool feature_edges_at_vertex_inside(size_t vid) const; + + /// Per-vertex positional constraint, on top of the envelopes -- the 3D counterpart of + /// TriOptimizerMesh's hook of the same name. An application uses this to pin a vertex to + /// a 0-dimensional feature it stands for, within a ball. Default: no constraint. + virtual bool smoothing_position_is_allowed(size_t, const Vector3d&) const { return true; } + /** * @brief The sentinel get_quality returns for an element AMIPS cannot score. * @@ -488,6 +558,8 @@ class TetOptimizerMesh : public wmtk::TetMesh, public wmtk::RationalPositions { double max_energy; std::map, FaceAttributes> changed_faces; + /// Feature-edge tags of the affected tets, by sorted vid pair (m_track_feature_edges). + std::map, bool> changed_edges; bool is_surface_flip = false; size_t sf_a = 0, sf_b = 0, sf_c = 0, sf_d = 0; @@ -501,9 +573,12 @@ class TetOptimizerMesh : public wmtk::TetMesh, public wmtk::RationalPositions size_t v2_id = 0; bool is_edge_on_surface = false; bool is_edge_open_boundary = false; + bool is_edge_on_feature = false; size_t edge_order = 0; double max_quality_before = 0.; std::vector>> changed_faces; + /// Feature-edge tags of the affected tets, by sorted vid pair (m_track_feature_edges). + std::map, bool> changed_edges; }; wmtk::threading::enumerable_thread_specific split_cache; @@ -518,12 +593,56 @@ class TetOptimizerMesh : public wmtk::TetMesh, public wmtk::RationalPositions std::vector> boundary_edges; std::vector changed_tids; std::vector changed_energies; + /// Feature-edge tags of the affected tets, by sorted vid pair (m_track_feature_edges). + std::map, bool> changed_edges; /// Coarsening pass only: the worst relative quality in the region the composite may /// disturb, measured before the collapse. See collapse_edge_after. double region_max_rel_before = 0.; }; wmtk::threading::enumerable_thread_specific collapse_cache; + // ---- Feature-edge tag propagation (m_track_feature_edges) ---------------------------- + // + // An edge attribute's slot is (lowest incident tet id) * 6 + local index, so it moves + // whenever that tet dies. Any edge whose slot can move is an edge of an affected tet + // (deleted tets are affected; created tets only affect their own six edges), so caching + // every edge of every affected tet before the operation and rewriting afterwards covers + // every slot that can move -- the same argument the face tags rely on. + // + // Restore comes in two forms because the operations differ in what they create: + // * split and the swaps CREATE tets, whose slots may be reused and hold stale values; + // restore_region walks every edge of every created tet and writes the cached value or + // the default. The created region provably contains every surviving cached edge. + // * collapse creates NO tets, and its link edges (la,lb) -- which contain neither + // endpoint, a case a tet FACE cannot exhibit -- can live on in tets outside any + // enumerable region. restore_remap walks the CACHE instead: remaps v1 -> v2, OR-merges + // colliding pairs, and writes each surviving edge wherever its slot now is. Returns + // false if a TAGGED cached edge no longer exists (the caller aborts, mirroring the + // face path's try_tuple_from_face bailout). + + /// Cache the feature tag of every edge of every tet in `tids`. + void feature_edges_cache( + const std::vector& tids, + std::map, bool>& cache); + /// Write back tags over every edge of every tet in `tets` (created regions). + void feature_edges_restore_region( + const std::vector& tets, + const std::map, bool>& cache); + void feature_edges_restore_region( + const std::vector& tids, + const std::map, bool>& cache); + /// Write back tags for every cached edge, with `v_old` renamed to `v_new` (collapse). + bool feature_edges_restore_remap( + const std::map, bool>& cache, + size_t v_old, + size_t v_new); + + /// Hook: the split vertex was created on a feature edge (or not). Runs only when + /// m_track_feature_edges is set; the counterpart of split_after_vertex for the feature + /// flag, separate so existing overrides keep their signature. + virtual void split_after_vertex_feature(size_t, bool) {} + + /// Set for the duration of coarsen_mesh(); read-only while a pass is running. bool m_coarsen_mode = false; diff --git a/src/wmtk/TetOptimizerMeshCollapse.cpp b/src/wmtk/TetOptimizerMeshCollapse.cpp index 071370b62e..4659fcd699 100644 --- a/src/wmtk/TetOptimizerMeshCollapse.cpp +++ b/src/wmtk/TetOptimizerMeshCollapse.cpp @@ -229,6 +229,42 @@ bool TetOptimizerMesh::collapse_edge_before(const Tuple& loc) // input is an edg cache.changed_faces.push_back(std::make_pair(f_attr, f_vids)); } + // Feature-edge tags: cache every edge of every tet incident to v1 -- the tets that get + // rewired or deleted. This includes the deleted tets' link edges (la,lb), which contain + // neither endpoint and can live on in tets outside any region enumerable after the + // collapse; that is why the restore walks the cache, not a region. + if (m_track_feature_edges) { + feature_edges_cache(n1_locs, cache.changed_edges); + + // The geometric guard: every tagged edge, as the collapse would leave it (v1 renamed + // to v2, which does not move), must stay inside the feature-curve tube. This is what + // lets a curve coarsen -- a collapse ALONG the curve keeps the merged edge in the + // tube and passes -- while a collapse that would drag the curve sideways, or weld + // two curves farther apart than the tube width, is refused. The exact counterpart of + // the surface-split's surface_triangle_is_outside check, and of the open-boundary + // collapse guard. + if (m_feature_envelope) { + for (const auto& [pair, tag] : cache.changed_edges) { + if (!tag) { + continue; + } + std::array p = pair; + for (size_t& v : p) { + if (v == v1_id) { + v = v2_id; + } + } + if (p[0] == p[1]) { + continue; // the collapsed edge itself; it disappears + } + if (m_feature_envelope->is_outside(std::array{ + {m_vertex_attribute[p[0]].m_posf, m_vertex_attribute[p[1]].m_posf}})) { + return false; + } + } + } + } + if (VA[v1_id].m_is_on_surface) { // this code must check if a face is tagged as surface face // only checking the vertices is not enough @@ -431,6 +467,15 @@ bool TetOptimizerMesh::collapse_edge_after(const Tuple& loc) m_face_attribute[std::get<1>(found.value())] = f_attr; } + // Feature-edge tags: cache-driven restore with v1 renamed to v2 (see the cache comment + // in collapse_edge_before). False means a tagged edge no longer exists -- abort, the + // rollback undoes everything. + if (m_track_feature_edges) { + if (!feature_edges_restore_remap(cache.changed_edges, v1_id, v2_id)) { + return false; + } + } + if (!m_coarsen_mode) { return true; } diff --git a/src/wmtk/TetOptimizerMeshSplit.cpp b/src/wmtk/TetOptimizerMeshSplit.cpp index b1422e3dcb..0055690ae7 100644 --- a/src/wmtk/TetOptimizerMeshSplit.cpp +++ b/src/wmtk/TetOptimizerMeshSplit.cpp @@ -179,6 +179,17 @@ bool TetOptimizerMesh::split_edge_before(const Tuple& loc0) } wmtk::vector_unique(cache.changed_faces, comp, is_equal); + if (m_track_feature_edges) { + std::vector tids; + tids.reserve(tets.size()); + for (const Tuple& t : tets) { + tids.push_back(t.tid(*this)); + } + feature_edges_cache(tids, cache.changed_edges); + cache.is_edge_on_feature = + m_feature_edge_attribute[loc0.eid(*this)].m_is_feature_edge; + } + return split_before_cells(loc0, tets); } @@ -363,6 +374,20 @@ bool TetOptimizerMesh::split_edge_after(const Tuple& loc) } } + // Feature-edge tags: rewrite every edge of the children (locs is exactly the one-ring of + // the new vertex, i.e. the created tets), then the split edge's two halves inherit its + // tag -- they are new pairs, so the region restore defaulted them. + if (m_track_feature_edges) { + feature_edges_restore_region(locs, cache.changed_edges); + if (cache.is_edge_on_feature) { + const Tuple e1 = tuple_from_edge({{v1_id, v_id}}); + const Tuple e2 = tuple_from_edge({{v2_id, v_id}}); + m_feature_edge_attribute[e1.eid(*this)].m_is_feature_edge = true; + m_feature_edge_attribute[e2.eid(*this)].m_is_feature_edge = true; + } + split_after_vertex_feature(v_id, cache.is_edge_on_feature); + } + m_vertex_attribute[v_id].partition_id = m_vertex_attribute[v1_id].partition_id; m_vertex_attribute[v_id].m_sizing_scalar = (m_vertex_attribute[v1_id].m_sizing_scalar + m_vertex_attribute[v2_id].m_sizing_scalar) / 2; diff --git a/src/wmtk/TetOptimizerMeshSwaps.cpp b/src/wmtk/TetOptimizerMeshSwaps.cpp index 2de486908f..c3cd25b05b 100644 --- a/src/wmtk/TetOptimizerMeshSwaps.cpp +++ b/src/wmtk/TetOptimizerMeshSwaps.cpp @@ -130,6 +130,10 @@ bool TetOptimizerMesh::swap_edge_before(const Tuple& t) if (is_edge_on_bbox(t)) { return false; } + // A tagged feature edge IS the operated edge here, and this swap deletes it outright. + if (m_track_feature_edges && m_feature_edge_attribute[t.eid(*this)].m_is_feature_edge) { + return false; + } // Surface edges are allowed only as a topology-preserving surface diagonal flip (see // prepare_surface_flip). If disabled, keep the old behavior of rejecting all surface-edge // swaps. Route on the direct incident-surface-face count so a genuine surface edge is never @@ -149,6 +153,10 @@ bool TetOptimizerMesh::swap_edge_before(const Tuple& t) cache.max_energy = max_energy; face_attribute_tracker(*this, incident_tets, m_face_attribute, cache.changed_faces); + if (m_track_feature_edges) { + feature_edges_cache(incident_tets, cache.changed_edges); + } + return true; } @@ -327,6 +335,10 @@ bool TetOptimizerMesh::swap_edge_after(const Tuple& t) } tracker_assign_after(*this, twotets, cache.changed_faces, m_face_attribute); + if (m_track_feature_edges) { + feature_edges_restore_region(twotets, cache.changed_edges); + } + if (cache.is_surface_flip) { // The generic tracker copied the old (interior) attributes onto the new @@ -436,6 +448,10 @@ bool TetOptimizerMesh::swap_face_before(const Tuple& t) } face_attribute_tracker(*this, twotets, m_face_attribute, cache.changed_faces); + if (m_track_feature_edges) { + feature_edges_cache(twotets, cache.changed_edges); + } + return true; } @@ -456,6 +472,10 @@ bool TetOptimizerMesh::swap_face_after(const Tuple& t) if (!swap_after_cells(new_tids, false)) return false; tracker_assign_after(*this, incident_tets, swap_cache.local().changed_faces, m_face_attribute); + if (m_track_feature_edges) { + feature_edges_restore_region(incident_tets, swap_cache.local().changed_edges); + } + cnt_swap++; return true; @@ -558,6 +578,10 @@ bool TetOptimizerMesh::swap_edge_44_before(const Tuple& t) if (is_edge_on_bbox(t)) { return false; } + // A tagged feature edge IS the operated edge here, and this swap deletes it outright. + if (m_track_feature_edges && m_feature_edge_attribute[t.eid(*this)].m_is_feature_edge) { + return false; + } // Surface edges are allowed only as a topology-preserving surface diagonal flip. The base 4-4 // swap is steered to the case that creates the new surface edge (c,d) by // swap_edge_44_accept_case; if no 4-4 diagonal yields (c,d) the swap is rejected. Route on the @@ -577,6 +601,10 @@ bool TetOptimizerMesh::swap_edge_44_before(const Tuple& t) cache.max_energy = max_energy; face_attribute_tracker(*this, incident_tets, m_face_attribute, cache.changed_faces); + if (m_track_feature_edges) { + feature_edges_cache(incident_tets, cache.changed_edges); + } + return true; } @@ -613,6 +641,10 @@ bool TetOptimizerMesh::swap_edge_44_after(const Tuple& t) } tracker_assign_after(*this, incident_tets, cache.changed_faces, m_face_attribute); + if (m_track_feature_edges) { + feature_edges_restore_region(incident_tets, cache.changed_edges); + } + if (cache.is_surface_flip) { // Re-tag the two new surface faces (the generic tracker reset them to interior). Net @@ -679,6 +711,10 @@ bool TetOptimizerMesh::swap_edge_56_before(const Tuple& t) if (is_edge_on_bbox(t)) { return false; } + // A tagged feature edge IS the operated edge here, and this swap deletes it outright. + if (m_track_feature_edges && m_feature_edge_attribute[t.eid(*this)].m_is_feature_edge) { + return false; + } // Surface edges are allowed only as a topology-preserving surface diagonal flip. The base 5-6 // swap is steered to the fan that creates the new surface edge (c,d) by // swap_edge_56_accept_case; if no fan yields (c,d) the swap is rejected. Route on the direct @@ -698,6 +734,10 @@ bool TetOptimizerMesh::swap_edge_56_before(const Tuple& t) cache.max_energy = max_energy; face_attribute_tracker(*this, incident_tets, m_face_attribute, cache.changed_faces); + if (m_track_feature_edges) { + feature_edges_cache(incident_tets, cache.changed_edges); + } + return true; } @@ -735,6 +775,10 @@ bool TetOptimizerMesh::swap_edge_56_after(const Tuple& t) } tracker_assign_after(*this, tids, cache.changed_faces, m_face_attribute); + if (m_track_feature_edges) { + feature_edges_restore_region(tids, cache.changed_edges); + } + if (cache.is_surface_flip) { // Re-tag the two new surface faces (the generic tracker reset them to interior). Net diff --git a/src/wmtk/optimization/SmoothVertex.hpp b/src/wmtk/optimization/SmoothVertex.hpp index d0efd980bd..0bf600b16f 100644 --- a/src/wmtk/optimization/SmoothVertex.hpp +++ b/src/wmtk/optimization/SmoothVertex.hpp @@ -225,8 +225,13 @@ bool smooth_vertex_3d( VA[vid].m_posf = x; }; + // A feature-curve vertex needs pulling and containment even when it is INTERIOR: the + // surface flag alone would smooth it as free and let it walk off its curve. Derived from + // the edge tags, so it cannot go stale. + const bool on_feature_curve = m.m_track_feature_edges && m.vertex_has_feature_edge(vid); const std::shared_ptr pull_env = - VA[vid].m_is_on_surface ? m.smoothing_energy_envelope(vid) : nullptr; + (VA[vid].m_is_on_surface || on_feature_curve) ? m.smoothing_energy_envelope(vid) + : nullptr; if (pull_env && opts.smoothing_mode == SmoothVertexOptions::SmoothingMode::Projected) { // Smooth as if the vertex were interior, then walk back onto the input. @@ -324,6 +329,20 @@ bool smooth_vertex_3d( solve(); } + // Per-vertex positional constraint (a 0-dimensional feature anchor's ball). Same hook + // and same placement as the 2D path. + if (!m.smoothing_position_is_allowed(vid, VA[vid].m_posf)) { + if (counters) ++counters->envelope; + return false; + } + + // Feature containment: every tagged edge at this vertex must stay inside the feature + // tube, at the candidate position (already written into the vertex attribute here). + if (on_feature_curve && !m.feature_edges_at_vertex_inside(vid)) { + if (counters) ++counters->envelope; + return false; + } + // Containment: every surface triangle at this vertex must still be inside. Checked // against the containment envelope, which is not necessarily the one it was pulled to. const std::shared_ptr check_env = diff --git a/src/wmtk/utils/io.hpp b/src/wmtk/utils/io.hpp index ba4323df02..e5d4255e77 100644 --- a/src/wmtk/utils/io.hpp +++ b/src/wmtk/utils/io.hpp @@ -22,6 +22,20 @@ namespace wmtk { class MshData { public: + template + void add_point_vertices(size_t num_vertices, const Fn& get_vertex_cb) + { + add_vertices<0>(num_vertices, get_vertex_cb); + } + + /// 0D point elements (gmsh element type 15), one per vertex of the last-added + /// (point) vertex block. The callback returns std::array. + template + void add_points(size_t num_points, const Fn& get_point_cb) + { + add_simplex_elements<0>(num_points, get_point_cb); + } + template void add_edge_vertices(size_t num_vertices, const Fn& get_vertex_cb) { @@ -398,7 +412,7 @@ class MshData template void add_vertices(size_t num_vertices, const Fn& get_vertex_cb) { - static_assert(DIM >= 1 && DIM <= 3, "Only 1,2,3D elements are supported!"); + static_assert(DIM >= 0 && DIM <= 3, "Only 0,1,2,3D elements are supported!"); if (num_vertices == 0) { logger().trace("Adding empty vertex block."); } @@ -429,8 +443,8 @@ class MshData void add_simplex_elements(size_t num_elements, const Fn& get_element_cb) { static_assert( - DIM == 1 || DIM == 2 || DIM == 3, - "Only 1,2,3D simplex elements are supported"); + DIM >= 0 && DIM <= 3, + "Only 0,1,2,3D simplex elements are supported"); if (num_elements == 0) return; if (m_spec.nodes.num_nodes == 0) { @@ -447,7 +461,9 @@ class MshData mshio::ElementBlock block; block.entity_dim = DIM; block.entity_tag = vertex_block.entity_tag; - if constexpr (DIM == 1) { + if constexpr (DIM == 0) { + block.element_type = 15; // 1-node point. + } else if constexpr (DIM == 1) { block.element_type = 1; // 2-node line. } else if constexpr (DIM == 2) { block.element_type = 2; // 3-node triangle. From c99e93c090e16c1da2fdd28558883ec1b9a8fd2e Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:24:20 -0400 Subject: [PATCH 4/7] Insertion: force feature edges and points through the arrangement, exactly The remesher's edge/point inputs (VolumeRemesher commit e09a229) were declared as empty locals and their provenance discarded. embed_triangles_in_tets now takes an EmbedFeaturesInput and returns an EmbedFeaturesResult with the edge tilings and point vertices remapped into the COMPACTED output numbering (the same v_map the tets go through), with contract checks that no feature vertex was compacted away. Passing features without the result struct throws: forcing features in and dropping where they went would leave them untaggable. Co-Authored-By: Claude Fable 5 --- src/wmtk/utils/EmbedTriangles.cpp | 48 +++++++++++++++++++++++++++++-- src/wmtk/utils/EmbedTriangles.hpp | 43 ++++++++++++++++++++++++++- 2 files changed, 87 insertions(+), 4 deletions(-) diff --git a/src/wmtk/utils/EmbedTriangles.cpp b/src/wmtk/utils/EmbedTriangles.cpp index 649b13a115..13ba60c2ce 100644 --- a/src/wmtk/utils/EmbedTriangles.cpp +++ b/src/wmtk/utils/EmbedTriangles.cpp @@ -30,8 +30,15 @@ void embed_triangles_in_tets( std::vector>& tets_after, std::vector& tet_face_on_input_surface, const EmbedTrianglesOptions& opts, - EmbedTrianglesProvenance* provenance) + EmbedTrianglesProvenance* provenance, + const EmbedFeaturesInput* features, + EmbedFeaturesResult* features_out) { + if (features != nullptr && features_out == nullptr) { + log_and_throw_error( + "embed_triangles_in_tets: features passed without features_out; forcing features " + "in and dropping where they went would leave them untaggable"); + } // Remesher outputs. The tet-based ones are what this consumes: out_tets (the // remesher's tetrahedra), final_tets_parent (parent polyhedral cell of each // tet), cells_with_faces_on_input (per-cell flag) and final_tets_parent_faces @@ -79,8 +86,11 @@ void embed_triangles_in_tets( // Step 3: run the exact arrangement. // volumeremesher embed - std::vector vr_edge_coords, vr_point_coords; - std::vector vr_edge_indexes; + static const std::vector no_coords; + static const std::vector no_indices; + const std::vector& vr_edge_coords = features ? features->edge_vrt_coord : no_coords; + const std::vector& vr_edge_indexes = features ? features->edge_indices : no_indices; + const std::vector& vr_point_coords = features ? features->point_coord : no_coords; std::vector>> vr_tri_provenance; std::vector vr_tri_group; std::vector>> vr_edge_provenance; @@ -380,6 +390,38 @@ void embed_triangles_in_tets( t[i] = v_map[t[i]]; } } + + // Feature provenance, into the compacted numbering. Every vertex on a forced edge or + // point is a vertex of some output tet -- that is what "forced into the output" means -- + // so v_map cannot be -1 for any of them; the throw is a remesher-contract check, not a + // reachable branch. + if (features_out != nullptr) { + features_out->edge_tiling.assign(vr_edge_provenance.size(), {}); + for (size_t e = 0; e < vr_edge_provenance.size(); ++e) { + auto& tiling = features_out->edge_tiling[e]; + tiling.reserve(vr_edge_provenance[e].size()); + for (const auto& seg : vr_edge_provenance[e]) { // {tet, v0, v1} + const int64_t a = v_map[seg[1]]; + const int64_t b = v_map[seg[2]]; + if (a < 0 || b < 0) { + log_and_throw_error( + "Feature edge {}: an output edge vertex was compacted away", e); + } + tiling.push_back({{size_t(a), size_t(b)}}); + } + } + features_out->point_vertex.assign(vr_point_provenance.size(), -1); + for (size_t p = 0; p < vr_point_provenance.size(); ++p) { + const uint32_t v = vr_point_provenance[p][1]; + if (v == UINT32_MAX) { + continue; + } + if (v_map[v] < 0) { + log_and_throw_error("Feature point {}: its output vertex was compacted away", p); + } + features_out->point_vertex[p] = v_map[v]; + } + } logger().info("done"); // Step 5: publish the tets. makeTetrahedra already emits WMTK-positively diff --git a/src/wmtk/utils/EmbedTriangles.hpp b/src/wmtk/utils/EmbedTriangles.hpp index 486809c27d..347a11168a 100644 --- a/src/wmtk/utils/EmbedTriangles.hpp +++ b/src/wmtk/utils/EmbedTriangles.hpp @@ -48,6 +48,40 @@ struct EmbedTrianglesProvenance std::vector triangle_group; }; +/** + * @brief Extra 1D/0D features to force into the output mesh, exactly. + * + * The remesher pins each edge as the shared crease of two forcing triangles and each point + * as the apex of a corner of three, then tracks them through the arrangement like the + * surface. Flat arrays, same conventions as the surface inputs. Everything -- the features + * AND the forcing-triangle apexes, which stick out by up to 0.1x the edge length -- must lie + * inside the background tet mesh's domain. + */ +struct EmbedFeaturesInput +{ + /// feature-edge vertices, 3 doubles each, xyz-interleaved + std::vector edge_vrt_coord; + /// feature edges, 2 vertex ids each, into edge_vrt_coord + std::vector edge_indices; + /// feature points, 3 doubles each, xyz-interleaved + std::vector point_coord; +}; + +/** + * @brief Where each feature ended up, in the COMPACTED output numbering (the same ids + * `tets_after` uses). + */ +struct EmbedFeaturesResult +{ + /// Per input edge: the output tet edges tiling it, as vertex pairs. Empty only if the + /// edge was degenerate. + std::vector>> edge_tiling; + /// Per input point: the output vertex exactly equal to it, or -1 if it was dropped + /// (which the remesher does not do for a valid point -- the -1 exists so a caller can + /// assert that rather than assume it). + std::vector point_vertex; +}; + /** * @brief Conformally insert a triangle soup into a background tet mesh, exactly. * @@ -82,6 +116,11 @@ struct EmbedTrianglesProvenance * @param[out] tets_after output tets * @param[out] tet_face_on_input_surface 4 flags per tet, in WMTK local face order * @param[out] provenance optional: which input triangles each surface face came from + * @param features optional: feature edges and points to force into the output + * @param[out] features_out optional: where each feature ended up, in the compacted output + * numbering. Required when `features` is non-null -- forcing + * features in and dropping where they went would leave the caller + * unable to tag them. */ void embed_triangles_in_tets( const std::vector& tri_vrt_coord, @@ -95,6 +134,8 @@ void embed_triangles_in_tets( std::vector>& tets_after, std::vector& tet_face_on_input_surface, const EmbedTrianglesOptions& opts = {}, - EmbedTrianglesProvenance* provenance = nullptr); + EmbedTrianglesProvenance* provenance = nullptr, + const EmbedFeaturesInput* features = nullptr, + EmbedFeaturesResult* features_out = nullptr); } // namespace wmtk::utils From 66a57db5c178cce95a9aef4e654f0ac5d3a266db Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:24:38 -0400 Subject: [PATCH 5/7] tetwild: insert, preserve, audit and extract feature edges and points New inputs input_edges / input_points: curve networks forced into the tetrahedralization as tet-edge chains and points as tet vertices, exactly; a vertex of an edge file with no incident edge is a free point, as in 2D. The input bounding box grows over the features before eps derives from it, and the insertion checks the forcing-triangle apexes fit the background box. Preservation: tiling edges tagged and endpoints flagged at init; the feature tube (feature_envelope_ratio, 0.5 with the order-2 rationale) guards collapses and smoothing; anchors (every input point, every open curve endpoint, junctions under allow_junction_cleanup -- triwild's model) give the 0-dimensional features the coverage that containment cannot: an open curve otherwise erodes from its own tips, every step invisibly inside the tube. Measured on the acceptance scene: coverage 0.0965 -> 0.0014 (< tube 0.0017), retention 5/5 vs 0/5 with the guard off, final energy unchanged. Audits at finalize, the curves' counterpart of the surface trio: two-sided feature deviation, per-component Euler characteristic of the feature network, DEBUG_feature_retention -- all fed from collections taken BEFORE any filter deletes tets, so a discarded region removes features from the tet mesh but never from the feature outputs (_features.obj) or the audits. filter='hybrid' with per-file input_roles (volume / surface): volume-role insides keep their tets by per-input winding, surface-role sheets survive as a triangle block, feature curves as a line block and anchors as 0D point elements in one mixed _hybrid.msh. Inherits filter='input''s documented barycenter-winding wrinkle; adopting tracked-winding semantics is listed in the PR's TODO. Tests: exact-rational tiling checker for the insertion (collinear, contiguous, endpoint to endpoint); a propagation fuzz whose discriminating power is itself verified (a deliberately disabled restore is caught; the uniform-pick version that could not catch it was strengthened until it did). Co-Authored-By: Claude Fable 5 --- .../wmtk/components/tetwild/Parameters.h | 11 + .../wmtk/components/tetwild/TetWildMesh.cpp | 132 +++++ .../wmtk/components/tetwild/TetWildMesh.h | 123 ++++- .../tetwild/VolumemesherInsertion.cpp | 160 ++++++- .../components/tetwild/tests/CMakeLists.txt | 1 + .../tetwild/tests/test_feature_tags.cpp | 295 ++++++++++++ .../tetwild/tests/test_insertion.cpp | 198 +++++++- .../wmtk/components/tetwild/tetwild.cpp | 452 +++++++++++++++++- .../wmtk/components/tetwild/tetwild_spec.json | 69 ++- 9 files changed, 1430 insertions(+), 11 deletions(-) create mode 100644 components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp diff --git a/components/tetwild/wmtk/components/tetwild/Parameters.h b/components/tetwild/wmtk/components/tetwild/Parameters.h index 9089bd596c..230042ad11 100644 --- a/components/tetwild/wmtk/components/tetwild/Parameters.h +++ b/components/tetwild/wmtk/components/tetwild/Parameters.h @@ -10,6 +10,17 @@ struct Parameters : public wmtk::OptimizerParameters /// surface envelope's. Deliberately below 1 where the surface envelope uses the full eps; /// see the doc on /order2_envelope_ratio in the spec for the measurement behind 0.5. double order2_envelope_ratio = 0.5; + /// Envelope thickness for the feature-curve tube, as a fraction of eps. Same value and + /// same rationale as order2_envelope_ratio: a curve has no collapse-blockage a wider + /// envelope would relieve, so the extra room is only geometry to resolve. + double feature_envelope_ratio = 0.5; + /// Keep the 0-dimensional features -- input points, and feature-curve endpoints (and + /// junctions, see allow_junction_cleanup) -- within m_feature_eps of where the input + /// put them. Mirrors triwild's parameter of the same name. + bool preserve_feature_points = true; + /// Anchor only the endpoints of open feature curves, letting junctions merge. Mirrors + /// triwild: the erosion the anchor exists for is an ENDPOINT property. + bool allow_junction_cleanup = true; Vector3d min = Vector3d::Zero(); Vector3d max = Vector3d::Ones(); Vector3d box_min = Vector3d::Zero(); diff --git a/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp b/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp index 31eb238c89..deb5596d3d 100644 --- a/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp +++ b/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp @@ -65,8 +65,50 @@ void TetWildMesh::optimization_sanity_checks_extra() } } +std::pair TetWildMesh::feature_retention(double* worst_ratio) const +{ + if (worst_ratio) { + *worst_ratio = 0; + } + if (m_feature_points.empty()) { + return {0, 0}; + } + + // One nearest-neighbour query per anchor against a kd-tree of the live vertices; the + // 2D version learned the hard way not to scan every vertex per feature. + std::vector pts; + pts.reserve(vert_capacity()); + for (const Tuple& v : get_vertices()) { + pts.push_back(m_vertex_attribute[v.vid(*this)].m_posf); + } + if (pts.empty()) { + return {0, m_feature_points.size()}; + } + const KNN knn(pts); + + size_t kept = 0; + for (const Vector3d& anchor : m_feature_points) { + uint32_t idx = 0; + double sq = 0; + knn.nearest_neighbor(anchor, idx, sq); + if (sq <= m_feature_eps * m_feature_eps) { + ++kept; + } else if (worst_ratio) { + *worst_ratio = std::max(*worst_ratio, std::sqrt(sq) / m_feature_eps); + } + } + return {kept, m_feature_points.size()}; +} + std::shared_ptr TetWildMesh::smoothing_energy_envelope(const size_t vid) const { + // A vertex on an input feature curve is pulled toward the feature tube. Checked FIRST: + // the feature is user-supplied where order 2 is derived, so where a vertex is both, the + // explicit constraint wins. Both are soft pulls; the hard vetoes (the collapse guard and + // the tagged-edge containment in smoothing) are unaffected by this order. + if (m_track_feature_edges && m_feature_envelope && vertex_has_feature_edge(vid)) { + return m_feature_envelope; + } // Order 2 means the vertex is on a surface boundary or a non-manifold edge. This is // broader than the old m_is_on_open_boundary flag, which covered only open boundaries. if (get_order_of_vertex(vid) >= 2 && m_order2_envelope && m_order2_envelope->initialized()) { @@ -626,6 +668,96 @@ void TetWildMesh::filter_with_input_surface_winding_number() remove_tets_by_ids(rm_tids); } +void TetWildMesh::filter_with_roles(const std::vector& volume_inputs) +{ + std::vector rm_tids; + for (const Tuple& t : get_tets()) { + const size_t tid = t.tid(*this); + const auto& wn = m_tet_attribute[tid].m_winding_number_per_input; + bool inside = false; + for (const size_t k : volume_inputs) { + if (k < wn.size() && wn[k] > 0.5) { + inside = true; + break; + } + } + if (!inside) { + rm_tids.emplace_back(tid); + } + } + remove_tets_by_ids(rm_tids); +} + +void TetWildMesh::output_hybrid_mesh( + const std::string& file, + const std::vector>& sheet_triangles, + const std::vector>& curve_segments, + const std::vector& anchor_points) +{ + wmtk::MshData msh; + + // Tet block: the filtered mesh as it stands (the caller consolidated it). + const auto& vtx = get_vertices(); + msh.add_tet_vertices(vtx.size(), [&](size_t k) { + return m_vertex_attribute[vtx[k].vid(*this)].m_posf; + }); + const auto& tets = get_tets(); + msh.add_tets(tets.size(), [&](size_t k) { + const auto vs = oriented_tet_vertices(tets[k]); + std::array data; + for (int j = 0; j < 4; j++) { + data[j] = vs[j].vid(*this); + } + return data; + }); + + // Face and edge blocks carry their own vertex arrays (welded by exact position). + const auto weld = [](const auto& elements, auto& verts, auto& indexed) { + std::map, size_t> vid_of; + for (const auto& el : elements) { + auto& out = indexed.emplace_back(); + for (size_t j = 0; j < el.size(); ++j) { + const std::array key = {{el[j][0], el[j][1], el[j][2]}}; + const auto [it, inserted] = vid_of.emplace(key, verts.size()); + if (inserted) { + verts.push_back(el[j]); + } + out[j] = it->second; + } + } + }; + if (!sheet_triangles.empty()) { + std::vector fv; + std::vector> ft; + weld(sheet_triangles, fv, ft); + msh.add_face_vertices(fv.size(), [&](size_t k) { return fv[k]; }); + msh.add_faces(ft.size(), [&](size_t k) { return ft[k]; }); + } + if (!curve_segments.empty()) { + std::vector ev; + std::vector> ee; + weld(curve_segments, ev, ee); + msh.add_edge_vertices(ev.size(), [&](size_t k) { return ev[k]; }); + msh.add_edges(ee.size(), [&](size_t k) { return ee[k]; }); + } + + if (!anchor_points.empty()) { + msh.add_point_vertices(anchor_points.size(), [&](size_t k) { return anchor_points[k]; }); + msh.add_points(anchor_points.size(), [&](size_t k) { + return std::array{{k}}; + }); + } + + msh.save(file, /*binary=*/true); + logger().info( + "hybrid output: {} tets, {} surface triangles, {} curve edges, {} anchor points -> {}", + tets.size(), + sheet_triangles.size(), + curve_segments.size(), + anchor_points.size(), + file); +} + void TetWildMesh::filter_with_tracked_surface_winding_number() { std::vector rm_tids; diff --git a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h index f680cf1869..22731dfc71 100644 --- a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h +++ b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -62,9 +63,37 @@ class TetWildMesh : public wmtk::TetOptimizerMesh /// counterpart is TriWildMesh::VertexExtras::m_feature_id, which differs deliberately /// -- see the comment there. bool m_is_on_open_boundary = false; + /// Whether the vertex lies on an input feature curve (see EdgeAttributes). + bool m_is_on_feature_curve = false; + /** + * Index into TetWildMesh::m_feature_points, or NO_FEATURE. The port of + * TriWildMesh::VertexExtras::m_feature_id, and for the same reason: preserving a + * 0-dimensional feature -- an input point, or a feature curve's endpoint -- is a + * COVERAGE property, and a containment test passes trivially when one anchor + * vertex is collapsed onto another while the first quietly stops being + * represented. See the discussion there. + */ + size_t m_feature_point_id = std::numeric_limits::max(); }; wmtk::AttributeCollection m_vertex_extra; + + /// A vertex that stands for no input feature point. See VertexExtras::m_feature_point_id. + static constexpr size_t NO_FEATURE = std::numeric_limits::max(); + + /** + * Anchor positions of the 0-dimensional features, indexed by + * VertexExtras::m_feature_point_id: the input feature points, and the feature-curve + * endpoints (junctions too, when allow_junction_cleanup is off). A vertex carrying + * feature f may never end up further than m_feature_eps from m_feature_points[f]; the + * anchor never moves, so the id is a direct lookup. + */ + std::vector m_feature_points; + /// Radius of the anchor ball. Set by the driver alongside m_feature_envelope; kept a + /// separate name from the tube eps for the same reason triwild keeps m_feature_eps + /// separate from m_envelope_eps. + double m_feature_eps = -1; + /// The base holds only wmtk::OptimizerParameters; this is the same object, typed, for the /// tetwild-only fields. Parameters& m_tet_params; @@ -94,10 +123,39 @@ class TetWildMesh : public wmtk::TetOptimizerMesh } bool collapse_before_vertex(size_t v1, size_t v2, double edge_length) override { + if (collapse_breaks_feature_point(v1, v2)) return false; if (edge_length <= 0 || !m_vertex_extra[v1].m_is_on_open_boundary) return true; return m_vertex_extra[v2].m_is_on_open_boundary || !m_order2_envelope->is_outside(m_vertex_attribute[v2].m_posf); } + + /** + * True iff collapsing v1 into v2 would drop or displace a feature-point anchor. + * Port of TriWildMesh::collapse_breaks_feature, distance test included: two anchors + * within m_feature_eps of each other are allowed to merge -- the survivor covers both + * -- and forbidding that deadlocks the mesh (measured in 2D, see the triwild comment). + */ + bool collapse_breaks_feature_point(const size_t v1_id, const size_t v2_id) const + { + if (!m_tet_params.preserve_feature_points) return false; + const size_t f1 = m_vertex_extra[v1_id].m_feature_point_id; + if (f1 == NO_FEATURE) return false; + return (m_vertex_attribute[v2_id].m_posf - m_feature_points[f1]).squaredNorm() > + m_feature_eps * m_feature_eps; + } + + bool smoothing_position_is_allowed(const size_t vid, const Vector3d& p) const override + { + const size_t f = m_vertex_extra[vid].m_feature_point_id; + if (f == NO_FEATURE || !m_tet_params.preserve_feature_points) return true; + // A ball, not a pin -- the vertex may move anywhere within m_feature_eps of its + // anchor, exactly the freedom the envelope gives everywhere else. + return (p - m_feature_points[f]).squaredNorm() <= m_feature_eps * m_feature_eps; + } + + /// {anchors still represented within m_feature_eps, total}. Geometric, like triwild's: + /// "some vertex is within eps of this anchor", whether or not it carries the id. + std::pair feature_retention(double* worst_ratio = nullptr) const; bool collapse_is_order_2_edge(const std::array& e) override { return is_open_boundary_edge(e); @@ -109,15 +167,36 @@ class TetWildMesh : public wmtk::TetOptimizerMesh { m_vertex_extra[v2].m_is_on_open_boundary = m_vertex_extra[v1].m_is_on_open_boundary || m_vertex_extra[v2].m_is_on_open_boundary; + m_vertex_extra[v2].m_is_on_feature_curve = + m_vertex_extra[v1].m_is_on_feature_curve || m_vertex_extra[v2].m_is_on_feature_curve; return true; } - void collapse_after_vertex(size_t, size_t v2) override + void collapse_after_vertex(size_t v1, size_t v2) override { + // The survivor picks up the anchor the collapsed vertex carried, if it has none. + if (m_vertex_extra[v2].m_feature_point_id == NO_FEATURE) { + m_vertex_extra[v2].m_feature_point_id = m_vertex_extra[v1].m_feature_point_id; + } if (m_vertex_extra[v2].m_is_on_open_boundary && !is_vertex_on_boundary(v2)) { m_vertex_extra[v2].m_is_on_open_boundary = false; } + // Same cleanup for the feature flag: v1's flag was OR-ed in (see + // collapse_after_connectivity), but if the survivor has no incident tagged edge the + // curve rerouted around it and the flag is stale. + if (m_track_feature_edges && m_vertex_extra[v2].m_is_on_feature_curve && + !vertex_has_feature_edge(v2)) { + m_vertex_extra[v2].m_is_on_feature_curve = false; + } + } + + void split_after_vertex_feature(size_t vid, bool on_feature) override + { + m_vertex_extra[vid].m_is_on_feature_curve = on_feature; + // A fresh split vertex stands for no anchor; its slot may be reused and stale. + m_vertex_extra[vid].m_feature_point_id = NO_FEATURE; } + /// Envelope a vertex is pulled toward while smoothing. std::shared_ptr smoothing_energy_envelope(const size_t vid) const override; TetWildMesh( @@ -224,6 +303,26 @@ class TetWildMesh : public wmtk::TetOptimizerMesh void filter_with_tracked_surface_winding_number(); void filter_with_flood_fill(); + /** + * The 'hybrid' filter: keep only tets inside (winding > 0.5) at least one of the + * inputs listed in `volume_inputs` (indices into the input list, matching + * m_winding_number_per_input). Everything else is scaffolding; the driver collects + * the surface-role faces and the feature curves BEFORE calling this, because they + * live on the tets this deletes. + */ + void filter_with_roles(const std::vector& volume_inputs); + + /** + * Write the mixed-dimensional output of the hybrid filter: the (already filtered) + * tet mesh, plus the collected surface-role triangles as a face block and the + * feature curves as an edge block and the anchor points as 0D point elements. + */ + void output_hybrid_mesh( + const std::string& file, + const std::vector>& sheet_triangles, + const std::vector>& curve_segments, + const std::vector& anchor_points); + // debug use std::atomic cnt_split = 0, cnt_collapse = 0; @@ -250,6 +349,14 @@ class TetWildMesh : public wmtk::TetOptimizerMesh * * This is the insertion path. See the banner in VolumemesherInsertion.cpp. */ + /** + * Optional feature inputs: `feature_edge_vertices`/`feature_edges` is an edge mesh to + * force into the tetrahedralization as tet edges, `feature_points` are points to force + * in as tet vertices. Where they ended up comes back in `features_out` (compacted + * output ids, see utils::EmbedFeaturesResult), which is required when features are + * passed. Features and their forcing-triangle apexes must lie inside the background + * box; the caller grows the input bounding box over them (tetwild.cpp does). + */ void insertion_by_volumeremesher( const std::vector& vertices, const std::vector>& faces, @@ -257,14 +364,24 @@ class TetWildMesh : public wmtk::TetOptimizerMesh std::vector>& facets_after, std::vector& is_v_on_input, std::vector>& tets_after, - std::vector& tet_face_on_input_surface); + std::vector& tet_face_on_input_surface, + const std::vector& feature_edge_vertices = {}, + const std::vector>& feature_edges = {}, + const std::vector& feature_points = {}, + utils::EmbedFeaturesResult* features_out = nullptr); + /** + * `features` tags the feature-edge tilings on the initial mesh; passing it registers + * m_edge_attribute with p_edge_attrs. Null => no feature tracking, nothing allocated. + */ void init_from_Volumeremesher( const std::vector& v_rational, const std::vector>& facets, const std::vector& is_v_on_input, const std::vector>& tets, - const std::vector& tet_face_on_input_surface); + const std::vector& tet_face_on_input_surface, + const utils::EmbedFeaturesResult* features = nullptr, + const std::vector* feature_anchors = nullptr); void init_from_file(std::string input_dir); diff --git a/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp b/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp index ceb0a6a987..6d56ed7123 100644 --- a/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp +++ b/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp @@ -45,9 +45,20 @@ void TetWildMesh::insertion_by_volumeremesher( std::vector>& polygon_faces, // out: triangular facets std::vector& is_v_on_input, // out: vertex-on-input-surface flags std::vector>& tets_after, // out: output tets - std::vector& tet_face_on_input_surface) // out: 4 face-on-surface flags per tet + std::vector& tet_face_on_input_surface, // out: 4 face-on-surface flags per tet + const std::vector& feature_edge_vertices, // optional feature-edge vertices + const std::vector>& feature_edges, // optional feature edges + const std::vector& feature_points, // optional feature points + utils::EmbedFeaturesResult* features_out) // out: where the features ended up { logger().info("Insertion Surface: #V = {}, #F = {}", vertices.size(), faces.size()); + if (!feature_edges.empty() || !feature_points.empty()) { + logger().info( + "Insertion features: #E = {} ({} vertices), #P = {}", + feature_edges.size(), + feature_edge_vertices.size(), + feature_points.size()); + } // Step 1: build the Delaunay background mesh. // generate background mesh @@ -118,6 +129,67 @@ void TetWildMesh::insertion_by_volumeremesher( opts.check_collinear_input = m_params.perform_sanity_checks; opts.check_orientation = m_params.perform_sanity_checks; opts.check_surface_provenance = m_params.perform_sanity_checks; + + // Marshal the feature inputs into the remesher's flat arrays, and check the domain: the + // features AND their forcing-triangle apexes -- which stick out by up to 0.1x the edge + // length along an axis -- must lie inside the background box, or the remesher's + // non-Delaunay insertion cannot locate them. The caller grows the input bbox over the + // features, so a violation here means only one thing: a feature edge so long that a + // tenth of it exceeds the box padding (diag/15). Throw with that diagnosis rather than + // letting the remesher assert. + const bool has_features = !feature_edges.empty() || !feature_points.empty(); + utils::EmbedFeaturesInput vr_features; + if (has_features) { + const auto& lo = m_tet_params.box_min; + const auto& hi = m_tet_params.box_max; + const auto in_box = [&](const Vector3d& p, double margin) { + for (int d = 0; d < 3; ++d) { + if (p[d] - margin < lo[d] || p[d] + margin > hi[d]) { + return false; + } + } + return true; + }; + vr_features.edge_vrt_coord.reserve(3 * feature_edge_vertices.size()); + for (const Vector3d& p : feature_edge_vertices) { + vr_features.edge_vrt_coord.insert(vr_features.edge_vrt_coord.end(), p.data(), p.data() + 3); + } + vr_features.edge_indices.reserve(2 * feature_edges.size()); + for (const auto& e : feature_edges) { + const Vector3d& a = feature_edge_vertices.at(e[0]); + const Vector3d& b = feature_edge_vertices.at(e[1]); + const double apex_margin = 0.1 * (a - b).norm(); + if (!in_box(a, apex_margin) || !in_box(b, apex_margin)) { + log_and_throw_error( + "Feature edge ({}, {}) or its forcing-triangle apexes (offset 0.1 x its " + "length {}) leave the background box [{}, {}]; the edge is too long " + "relative to the box padding", + a.transpose(), + b.transpose(), + (a - b).norm(), + lo.transpose(), + hi.transpose()); + } + vr_features.edge_indices.push_back(uint32_t(e[0])); + vr_features.edge_indices.push_back(uint32_t(e[1])); + } + // Point apexes are offset by 0.1x the estimated point spacing, itself bounded by the + // point set's own bbox diagonal -- small against the box padding; checked with the + // conservative margin of the padding itself would over-reject, so points are checked + // with no margin and the remesher's own debug assert stays as the backstop. + vr_features.point_coord.reserve(3 * feature_points.size()); + for (const Vector3d& p : feature_points) { + if (!in_box(p, 0)) { + log_and_throw_error( + "Feature point ({}) lies outside the background box [{}, {}]", + p.transpose(), + lo.transpose(), + hi.transpose()); + } + vr_features.point_coord.insert(vr_features.point_coord.end(), p.data(), p.data() + 3); + } + } + utils::embed_triangles_in_tets( tri_vrt_coord, triangle_indices, @@ -129,7 +201,10 @@ void TetWildMesh::insertion_by_volumeremesher( is_v_on_input, tets_after, tet_face_on_input_surface, - opts); + opts, + nullptr, + has_features ? &vr_features : nullptr, + has_features ? features_out : nullptr); // TODO this is a sanity check, but it is checked all the time for now, until insertion is @@ -165,8 +240,17 @@ void TetWildMesh::init_from_Volumeremesher( const std::vector>& facets, const std::vector& is_v_on_input, const std::vector>& tets, - const std::vector& tet_face_on_input_surface) + const std::vector& tet_face_on_input_surface, + const utils::EmbedFeaturesResult* features, + const std::vector* feature_anchors) { + // Turn on feature tracking BEFORE init so the connectivity init sizes the edge + // attributes along with everything else. Only when features exist: with the flag off the + // shared operations skip every edge-tag branch and nothing is allocated. + if (features != nullptr) { + enable_feature_edge_tracking(); + } + init_with_isolated_vertices(v_rational.size(), tets); assert(check_mesh_connectivity_validity()); @@ -231,6 +315,76 @@ void TetWildMesh::init_from_Volumeremesher( }, NUM_THREADS); + // Feature tags: mark each tiling edge with its input file id and flag its vertices. + // The tiling pairs come out of the same arrangement as the connectivity, so a missing + // edge is a contract violation, not an input problem -- throw, don't warn. + if (features != nullptr) { + size_t n_tagged = 0; + for (size_t e = 0; e < features->edge_tiling.size(); ++e) { + for (const auto& seg : features->edge_tiling[e]) { + const Tuple t = tuple_from_edge({{seg[0], seg[1]}}); + if (!t.is_valid(*this)) { + log_and_throw_error( + "Feature tiling edge ({}, {}) is not an edge of the initial mesh", + seg[0], + seg[1]); + } + m_feature_edge_attribute[t.eid(*this)].m_is_feature_edge = true; + m_vertex_extra[seg[0]].m_is_on_feature_curve = true; + m_vertex_extra[seg[1]].m_is_on_feature_curve = true; + ++n_tagged; + } + } + logger().info( + "feature tags: {} tet edges tagged across {} input edges", + n_tagged, + features->edge_tiling.size()); + + // Anchor the 0-dimensional features. Anchor positions are input coordinates, and + // every one of them exists in the output exactly -- input points via the point + // provenance, curve endpoints as tiling vertices -- and input coordinates are + // explicit points, so their rationals ARE doubles and an exact double lookup finds + // them. A vertex already carrying an id keeps it: one anchor per position is + // enough, the retention audit is geometric. + if (feature_anchors != nullptr && !feature_anchors->empty()) { + std::map, size_t> vid_of; + const auto key = [](const Vector3d& p) { + return std::array{{p[0], p[1], p[2]}}; + }; + for (const int64_t v : features->point_vertex) { + if (v >= 0) { + vid_of.emplace(key(m_vertex_attribute[size_t(v)].m_posf), size_t(v)); + } + } + for (const auto& tiling : features->edge_tiling) { + for (const auto& seg : tiling) { + for (const size_t v : {seg[0], seg[1]}) { + vid_of.emplace(key(m_vertex_attribute[v].m_posf), v); + } + } + } + size_t n_anchored = 0; + for (const Vector3d& p : *feature_anchors) { + const auto it = vid_of.find(key(p)); + if (it == vid_of.end()) { + log_and_throw_error( + "Feature anchor ({}) is not a vertex of the initial mesh", + p.transpose()); + } + if (m_vertex_extra[it->second].m_feature_point_id != NO_FEATURE) { + continue; + } + m_vertex_extra[it->second].m_feature_point_id = m_feature_points.size(); + m_feature_points.push_back(p); + ++n_anchored; + } + logger().info( + "feature anchors: {} points anchored within {:.6} of their input positions", + n_anchored, + m_feature_eps); + } + } + // track bounding box (parallel). The per-face exact-rational corner test is the // cost; on-bbox faces are rare, so each chunk collects (vertex, bbox-side) pairs // locally and merges once, and the per-vertex vectors are appended serially. diff --git a/components/tetwild/wmtk/components/tetwild/tests/CMakeLists.txt b/components/tetwild/wmtk/components/tetwild/tests/CMakeLists.txt index e8c32e3921..c96608ce63 100644 --- a/components/tetwild/wmtk/components/tetwild/tests/CMakeLists.txt +++ b/components/tetwild/wmtk/components/tetwild/tests/CMakeLists.txt @@ -10,6 +10,7 @@ add_component_test(${COMPONENT_NAME}) set(SRC_FILES test_file_write.cpp test_insertion.cpp + test_feature_tags.cpp test_operation_smooth.cpp test_operations.cpp test_surface_swap.cpp diff --git a/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp b/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp new file mode 100644 index 0000000000..e6f5b1fe90 --- /dev/null +++ b/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp @@ -0,0 +1,295 @@ +#include +#include +#include +#include + +#include + +#include +#include + +using namespace wmtk; +using namespace components::tetwild; + +// Feature-edge tag PROPAGATION through the shared operations. +// +// The failure this hunts is silent: an edge attribute's slot is (lowest incident tet) * 6 + +// local index, so a neighbouring operation that deletes that tet moves the slot, and a +// propagation bug leaves the tag at the dead one -- nothing crashes, the curve just loses a +// link. So the assertions here are about the TAGGED SET as a whole, re-read from the live +// slots after every single operation: +// +// * splits + swaps only (phase 1): the tagged set must still tile every input edge exactly +// -- collinear, contiguous, endpoint to endpoint, in rationals. Splits along the curve +// subdivide the tiling legally; swaps that would delete a tagged edge are vetoed; nothing +// else may change it. A dropped tag breaks coverage, a stale slot breaks collinearity. +// * with collapses (phase 2): collapses legally erode the curve (no anchors yet), so the +// invariant weakens to: every tagged edge lies EXACTLY on some input segment, and both its +// endpoints carry the vertex flag. Tags may vanish; they may never leak off the curve. + +namespace { + +struct FeatureScene +{ + Parameters params; // must outlive the mesh (held by reference) + std::shared_ptr mesh; + std::vector fe_verts; + std::vector> fe; + std::shared_ptr surf; +}; + +FeatureScene build_scene() +{ + FeatureScene sc; + // Cube [0,2]^3 with an interior 2-segment feature chain. + MatrixXd V(8, 3); + V << 0, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, // + 0, 0, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2; + MatrixXi F(12, 3); + F << 0, 2, 1, 0, 3, 2, 4, 5, 6, 4, 6, 7, 0, 1, 5, 0, 5, 4, // + 2, 3, 7, 2, 7, 6, 1, 2, 6, 1, 6, 5, 0, 4, 7, 0, 7, 3; + + std::vector vertices; + std::vector> faces; + VF_to_vectors(V, F, vertices, faces); + sc.params.init(vertices, faces); + + sc.surf = std::make_shared( + vertices, + 0); + { + std::vector frozen_verts; + sc.surf->create_mesh(vertices.size(), faces, frozen_verts, 0.1); + } + auto* env_ptr = &sc.surf->m_envelope; + const std::shared_ptr env(env_ptr, [](SampleEnvelope*) {}); + + sc.fe_verts = {{0.5, 0.5, 0.5}, {1.0, 1.0, 1.0}, {1.5, 1.0, 1.5}}; + sc.fe = {{{0, 1}}, {{1, 2}}}; + + std::vector v_rational; + std::vector> facets; + std::vector is_v_on_input; + std::vector> tets; + std::vector tet_face_on_input_surface; + utils::EmbedFeaturesResult features_out; + { + TetWildMesh insertion_mesh(sc.params, env, 0); + insertion_mesh.insertion_by_volumeremesher( + vertices, + faces, + v_rational, + facets, + is_v_on_input, + tets, + tet_face_on_input_surface, + sc.fe_verts, + sc.fe, + {}, + &features_out); + } + sc.mesh = std::make_shared(sc.params, env, 0); + { + std::vector fe_env(sc.fe.size()); + for (size_t i = 0; i < sc.fe.size(); ++i) { + fe_env[i] = Eigen::Vector2i(int(sc.fe[i][0]), int(sc.fe[i][1])); + } + sc.mesh->m_feature_envelope = std::make_shared(); + sc.mesh->m_feature_envelope->init(sc.fe_verts, fe_env, 0.05); + sc.mesh->m_feature_eps = 0.05; + } + // Anchor the chain's endpoints, as the driver would (valence-1 vertices). + const std::vector anchors = {sc.fe_verts[0], sc.fe_verts[2]}; + sc.mesh->init_from_Volumeremesher( + v_rational, + facets, + is_v_on_input, + tets, + tet_face_on_input_surface, + &features_out, + &anchors); + return sc; +} + +/// All currently tagged edges, re-read from the live slots. +std::vector> tagged_edges(TetWildMesh& m) +{ + std::vector> out; + for (const auto& e : m.get_edges()) { + if (!m.m_feature_edge_attribute[e.eid(m)].m_is_feature_edge) { + continue; + } + size_t a = e.vid(m); + size_t b = e.switch_vertex(m).vid(m); + if (a > b) { + std::swap(a, b); + } + out.push_back({{a, b}}); + } + return out; +} + +/// Exact parameter of vertex `vid` on segment (A,B), or nullopt if not exactly on it. +std::optional param_on_segment( + TetWildMesh& m, + const size_t vid, + const Vector3d& A, + const Vector3d& B) +{ + const Vector3r a{Rational(A[0]), Rational(A[1]), Rational(A[2])}; + const Vector3r d{Rational(B[0] - A[0]), Rational(B[1] - A[1]), Rational(B[2] - A[2])}; + const Vector3r p = m.m_vertex_attribute[vid].m_pos - a; + const Vector3r c = p.cross(d); + if (c[0] != Rational(0) || c[1] != Rational(0) || c[2] != Rational(0)) { + return std::nullopt; + } + const Rational t = p.dot(d) / d.dot(d); + if (t < Rational(0) || t > Rational(1)) { + return std::nullopt; + } + return t; +} + +/// Strict invariant (no collapses have run): every tagged edge lies EXACTLY on some input +/// segment, with flagged endpoints. +void require_tags_on_curve(FeatureScene& sc) +{ + for (const auto& e : tagged_edges(*sc.mesh)) { + bool on_some_segment = false; + for (const auto& seg : sc.fe) { + const auto t0 = + param_on_segment(*sc.mesh, e[0], sc.fe_verts[seg[0]], sc.fe_verts[seg[1]]); + const auto t1 = + param_on_segment(*sc.mesh, e[1], sc.fe_verts[seg[0]], sc.fe_verts[seg[1]]); + if (t0.has_value() && t1.has_value()) { + on_some_segment = true; + break; + } + } + REQUIRE(on_some_segment); + REQUIRE(sc.mesh->m_vertex_extra[e[0]].m_is_on_feature_curve); + REQUIRE(sc.mesh->m_vertex_extra[e[1]].m_is_on_feature_curve); + } +} + +/// Phase-2 invariant: collapses may re-route the curve inside its tube, so exact +/// collinearity no longer holds -- containment in the feature envelope is the contract the +/// collapse guard enforces, and is what is asserted. The anchors add the coverage half: +/// every anchor keeps a live vertex within its ball, after every single collapse. +void require_tags_in_tube(FeatureScene& sc) +{ + TetWildMesh& m = *sc.mesh; + for (const auto& e : tagged_edges(m)) { + REQUIRE(!sc.mesh->m_feature_envelope->is_outside(std::array{ + {m.m_vertex_attribute[e[0]].m_posf, m.m_vertex_attribute[e[1]].m_posf}})); + REQUIRE(m.m_vertex_extra[e[0]].m_is_on_feature_curve); + REQUIRE(m.m_vertex_extra[e[1]].m_is_on_feature_curve); + } + const auto [kept, total] = m.feature_retention(); + REQUIRE(total == 2); + REQUIRE(kept == 2); +} + +/// Phase-1 invariant: the tagged set tiles every input segment exactly. +void require_full_tiling(FeatureScene& sc) +{ + require_tags_on_curve(sc); + const auto tagged = tagged_edges(*sc.mesh); + for (const auto& seg : sc.fe) { + const Vector3d& A = sc.fe_verts[seg[0]]; + const Vector3d& B = sc.fe_verts[seg[1]]; + std::vector> intervals; + for (const auto& e : tagged) { + const auto t0 = param_on_segment(*sc.mesh, e[0], A, B); + const auto t1 = param_on_segment(*sc.mesh, e[1], A, B); + if (!t0.has_value() || !t1.has_value()) { + continue; + } + auto lo = *t0, hi = *t1; + if (hi < lo) { + std::swap(lo, hi); + } + if (lo == hi) { + continue; // an endpoint shared with the neighbouring input segment + } + intervals.emplace_back(lo, hi); + } + REQUIRE(!intervals.empty()); + std::sort(intervals.begin(), intervals.end(), [](const auto& x, const auto& y) { + return x.first < y.first; + }); + REQUIRE(intervals.front().first == Rational(0)); + for (size_t i = 1; i < intervals.size(); ++i) { + REQUIRE(intervals[i].first == intervals[i - 1].second); + } + REQUIRE(intervals.back().second == Rational(1)); + } +} + +} // namespace + +TEST_CASE("feature-tag-propagation-fuzz", "[tetwild_operation][features]") +{ + FeatureScene sc = build_scene(); + TetWildMesh& m = *sc.mesh; + require_full_tiling(sc); + + std::mt19937 rng(7); + // Half the picks are edges incident to a curve vertex: slot churn only endangers tags + // near the curve, and uniform picks on a growing mesh rarely go near it. Verified by + // negative control -- with uniform picks, a disabled split restore survived the fuzz. + const auto random_edge = [&]() -> TetMesh::Tuple { + if (rng() % 2 == 0) { + const auto tagged = tagged_edges(m); + if (!tagged.empty()) { + const auto& te = tagged[rng() % tagged.size()]; + const auto ring = m.get_one_ring_tets_for_vertex(m.tuple_from_vertex(te[0])); + if (!ring.empty()) { + const auto& t = ring[rng() % ring.size()]; + return m.tuple_from_edge(t.tid(m), int(rng() % 6)); + } + } + } + const auto edges = m.get_edges(); + return edges[rng() % edges.size()]; + }; + + // Phase 1: splits and swaps only. The curve must stay fully tiled after every operation. + size_t done = 0; + for (int i = 0; i < 1500 && done < 250; ++i) { + const TetMesh::Tuple e = random_edge(); + if (!e.is_valid(m)) { + continue; + } + bool changed = false; + std::vector new_tets; + switch (rng() % 4) { + case 0: changed = m.split_edge(e, new_tets); break; + case 1: changed = m.swap_edge(e, new_tets); break; + case 2: changed = m.swap_edge_44(e, new_tets); break; + case 3: changed = m.swap_edge_56(e, new_tets); break; + } + if (!changed) { + continue; + } + ++done; + require_full_tiling(sc); + } + REQUIRE(done >= 100); // the fuzz must actually have exercised operations + + // Phase 2: collapses join. Tags may erode; they may never leak off the curve. + size_t collapsed = 0; + for (int i = 0; i < 400 && collapsed < 40; ++i) { + const TetMesh::Tuple e = random_edge(); + if (!e.is_valid(m)) { + continue; + } + std::vector new_tets; + if (!m.collapse_edge(e, new_tets)) { + continue; + } + ++collapsed; + require_tags_in_tube(sc); + } + REQUIRE(collapsed >= 10); +} diff --git a/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp b/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp index fc350a759b..0b015c11ce 100644 --- a/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp +++ b/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp @@ -273,4 +273,200 @@ TEST_CASE("vertex_order", "[tetwild]") } CHECK(nvo3 == vo3_count); -} \ No newline at end of file +} +// --------------------------------------------------------------------------- +// Feature edges and points through the insertion, checked in exact arithmetic. +// +// The contract (VolumeRemesher commit e09a229): each input edge is represented +// by output tet edges that tile it -- collinear with it, contiguous, and +// running endpoint to endpoint -- and each input point by an output vertex +// exactly equal to it. These checks are the wmtk-side counterpart of the +// remesher's own verify_tracking, run on the COMPACTED ids the caller gets. +// --------------------------------------------------------------------------- +namespace { + +// Exact check that `tiling` tiles segment AB: every tiling vertex on the line +// through AB with parameter t in [0,1], and the pieces cover [0,1] with no gap. +void require_exact_tiling( + const std::vector>& tiling, + const std::vector& v_rational, + const Vector3d& seg_a, + const Vector3d& seg_b) +{ + REQUIRE(!tiling.empty()); + const Vector3r a{Rational(seg_a[0]), Rational(seg_a[1]), Rational(seg_a[2])}; + const Vector3r d{Rational(seg_b[0] - seg_a[0]), Rational(seg_b[1] - seg_a[1]), Rational(seg_b[2] - seg_a[2])}; + const Rational dd = d.dot(d); + + // Parameter of a vertex along AB, after requiring it exactly on the line. + const auto param = [&](size_t vid) -> Rational { + const Vector3r p = v_rational.at(vid) - a; + const Vector3r c = p.cross(d); + REQUIRE(c[0] == Rational(0)); + REQUIRE(c[1] == Rational(0)); + REQUIRE(c[2] == Rational(0)); + return p.dot(d) / dd; + }; + + std::vector> intervals; + intervals.reserve(tiling.size()); + for (const auto& e : tiling) { + Rational t0 = param(e[0]); + Rational t1 = param(e[1]); + if (t1 < t0) { + std::swap(t0, t1); + } + REQUIRE(t0 >= Rational(0)); + REQUIRE(t1 <= Rational(1)); + REQUIRE(t0 < t1); // no degenerate pieces + intervals.emplace_back(t0, t1); + } + std::sort(intervals.begin(), intervals.end(), [](const auto& x, const auto& y) { + return x.first < y.first; + }); + REQUIRE(intervals.front().first == Rational(0)); + for (size_t i = 1; i < intervals.size(); ++i) { + REQUIRE(intervals[i].first == intervals[i - 1].second); // contiguous, no overlap + } + REQUIRE(intervals.back().second == Rational(1)); +} + +} // namespace + +TEST_CASE("insertion-feature-edges-and-points", "[tetwild_operation][features]") +{ + // A cube [0,2]^3, with features interior to it, on its surface, and touching a corner. + MatrixXd V(8, 3); + V << 0, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, // + 0, 0, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2; + MatrixXi F(12, 3); + F << 0, 2, 1, 0, 3, 2, // z = 0 + 4, 5, 6, 4, 6, 7, // z = 2 + 0, 1, 5, 0, 5, 4, // y = 0 + 2, 3, 7, 2, 7, 6, // y = 2 + 1, 2, 6, 1, 6, 5, // x = 2 + 0, 4, 7, 0, 7, 3; // x = 0 + + std::vector vertices; + std::vector> faces; + VF_to_vectors(V, F, vertices, faces); + + Parameters params; + params.init(vertices, faces); + + components::shortest_edge_collapse::ShortestEdgeCollapse surf_mesh(vertices, 0); + { + std::vector frozen_verts; + surf_mesh.create_mesh(vertices.size(), faces, frozen_verts, 0.1); + } + const std::shared_ptr env(&surf_mesh.m_envelope, [](SampleEnvelope*) {}); + + // Features: an interior diagonal edge chain (two edges sharing a vertex), an edge lying + // ON the y = 0 face, an interior point, a point on the surface, and a point exactly at a + // cube corner (already a background vertex -- the dedup case). + const std::vector fe_verts = { + {0.5, 0.5, 0.5}, + {1.0, 1.0, 1.0}, + {1.5, 1.0, 1.5}, // interior chain 0-1, 1-2 + {0.5, 0.0, 0.5}, + {1.5, 0.0, 1.5}, // on-surface edge 3-4 + }; + const std::vector> fe = {{{0, 1}}, {{1, 2}}, {{3, 4}}}; + const std::vector fp = { + {1.0, 0.5, 1.5}, // interior + {2.0, 1.0, 1.0}, // on the x = 2 face + {0.0, 0.0, 0.0}, // exactly a cube corner + }; + + std::vector v_rational; + std::vector> facets; + std::vector is_v_on_input; + std::vector> tets; + std::vector tet_face_on_input_surface; + utils::EmbedFeaturesResult features_out; + { + TetWildMesh mesh_insertion(params, env, 0); + mesh_insertion.insertion_by_volumeremesher( + vertices, + faces, + v_rational, + facets, + is_v_on_input, + tets, + tet_face_on_input_surface, + fe_verts, + fe, + fp, + &features_out); + } + + // Every feature edge is tiled exactly, endpoint to endpoint. + REQUIRE(features_out.edge_tiling.size() == fe.size()); + for (size_t e = 0; e < fe.size(); ++e) { + require_exact_tiling( + features_out.edge_tiling[e], + v_rational, + fe_verts[fe[e][0]], + fe_verts[fe[e][1]]); + } + + // Every feature point is an output vertex, exactly. + REQUIRE(features_out.point_vertex.size() == fp.size()); + for (size_t p = 0; p < fp.size(); ++p) { + const int64_t vid = features_out.point_vertex[p]; + REQUIRE(vid >= 0); + for (int k = 0; k < 3; ++k) { + CHECK(v_rational[size_t(vid)][k] == Rational(fp[p][k])); + } + } + + // And the tiling edges are real tet edges of the output connectivity. + std::set> tet_edges; + for (const auto& t : tets) { + for (int i = 0; i < 4; ++i) { + for (int j = i + 1; j < 4; ++j) { + tet_edges.emplace(std::min(t[i], t[j]), std::max(t[i], t[j])); + } + } + } + for (const auto& tiling : features_out.edge_tiling) { + for (const auto& e : tiling) { + CHECK(tet_edges.count({std::min(e[0], e[1]), std::max(e[0], e[1])}) == 1); + } + } + + // Init the mesh with the features and check the tags landed: every tiling edge's slot + // is tagged and its endpoints are flagged. + TetWildMesh mesh(params, env, 0); + mesh.init_from_Volumeremesher( + v_rational, + facets, + is_v_on_input, + tets, + tet_face_on_input_surface, + &features_out); + + for (size_t ie = 0; ie < features_out.edge_tiling.size(); ++ie) { + for (const auto& e : features_out.edge_tiling[ie]) { + const auto t = mesh.tuple_from_edge({{e[0], e[1]}}); + REQUIRE(t.is_valid(mesh)); + CHECK(mesh.m_feature_edge_attribute[t.eid(mesh)].m_is_feature_edge); + CHECK(mesh.m_vertex_extra[e[0]].m_is_on_feature_curve); + CHECK(mesh.m_vertex_extra[e[1]].m_is_on_feature_curve); + } + } + // A vertex NOT on any feature curve is not flagged: count flagged vertices and compare + // with the union of tiling vertices. + std::set tiling_verts; + for (const auto& tiling : features_out.edge_tiling) { + for (const auto& e : tiling) { + tiling_verts.insert(e[0]); + tiling_verts.insert(e[1]); + } + } + size_t flagged = 0; + for (size_t v = 0; v < mesh.vert_capacity(); ++v) { + flagged += mesh.m_vertex_extra[v].m_is_on_feature_curve ? 1 : 0; + } + CHECK(flagged == tiling_verts.size()); +} diff --git a/components/tetwild/wmtk/components/tetwild/tetwild.cpp b/components/tetwild/wmtk/components/tetwild/tetwild.cpp index 9dca24e0b5..594613cca0 100644 --- a/components/tetwild/wmtk/components/tetwild/tetwild.cpp +++ b/components/tetwild/wmtk/components/tetwild/tetwild.cpp @@ -9,9 +9,11 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -85,6 +87,50 @@ std::vector compute_euler_characteristics(const MatrixXi& F) return euler_characteristics; } +/** + * @brief Euler characteristic (V - E) of each connected component of an edge network, + * sorted. The 1D counterpart of compute_euler_characteristics above, and the same check + * triwild runs on its curves: per component, 1 for a tree/open polyline, 0 for one loop, + * and so on. Comparing the sorted per-component values of the input feature network against + * the output's tagged edges catches exactly what the distance audits cannot -- components + * merged, split, or lost outright. + */ +std::vector curve_euler_characteristics(const std::vector>& edges) +{ + std::map root; // union-find over the vertices that appear + std::function find = [&](size_t v) -> size_t { + while (root[v] != v) { + root[v] = root[root[v]]; + v = root[v]; + } + return v; + }; + for (const auto& e : edges) { + for (const size_t v : {e[0], e[1]}) { + if (root.count(v) == 0) { + root[v] = v; + } + } + root[find(e[0])] = find(e[1]); + } + std::map> vc_ec; // component root -> (#V, #E) + for (const auto& [v, r] : root) { + (void)r; + ++vc_ec[find(v)].first; + } + for (const auto& e : edges) { + ++vc_ec[find(e[0])].second; + } + std::vector ecs; + ecs.reserve(vc_ec.size()); + for (const auto& [r, ve] : vc_ec) { + (void)r; + ecs.push_back(ve.first - ve.second); + } + std::sort(ecs.begin(), ecs.end()); + return ecs; +} + TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) { using wmtk::utils::resolve_path; @@ -107,9 +153,45 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) int max_its = json_params["max_iterations"]; std::string filter_option = json_params["filter"]; + // Roles for the 'hybrid' filter: one entry per input file, "volume" (tet-fill its + // inside) or "surface" (keep as an embedded sheet). Curves and points already declare + // their dimension by arriving through input_edges / input_points. + std::vector input_roles = json_params["input_roles"]; + std::vector volume_input_ids; + if (filter_option == "hybrid") { + const size_t n_inputs = json_params["input"].size(); + if (input_roles.empty()) { + input_roles.assign(n_inputs, "volume"); + } + if (input_roles.size() != n_inputs) { + log_and_throw_error( + "input_roles has {} entries for {} inputs", + input_roles.size(), + n_inputs); + } + for (size_t k = 0; k < input_roles.size(); ++k) { + if (input_roles[k] == "volume") { + volume_input_ids.push_back(k); + } else if (input_roles[k] != "surface") { + log_and_throw_error( + "input_roles[{}] = '{}'; must be 'volume' or 'surface'", + k, + input_roles[k]); + } + } + if (volume_input_ids.empty()) { + logger().warn("filter='hybrid' with no volume-role input: no tets will be kept."); + } + } else if (!input_roles.empty()) { + logger().warn("input_roles is only used by filter='hybrid'; ignoring it."); + } + params.epsr = json_params["eps_rel"]; params.lr = json_params["length_rel"]; params.order2_envelope_ratio = json_params["order2_envelope_ratio"]; + params.feature_envelope_ratio = json_params["feature_envelope_ratio"]; + params.preserve_feature_points = json_params["preserve_feature_points"]; + params.allow_junction_cleanup = json_params["allow_junction_cleanup"]; params.stop_energy = json_params["stop_energy"]; params.split_high_valence_threshold = json_params["split_high_valence_threshold"]; params.num_smoothing_passes = json_params["num_smoothing_passes"]; @@ -173,6 +255,114 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) box_minmax.second = V.colwise().maxCoeff(); VF_to_vectors(V, F, verts, tris); } + // Feature inputs: an edge mesh to force into the tetrahedralization as tet edges, and + // points to force in as tet vertices. Both join the bounding box BEFORE eps and the + // background box derive from it -- they are input geometry, and a feature outside the + // surface's box must still land inside the triangulated domain. + std::vector feature_edge_vertices; + std::vector> feature_edges; + std::vector feature_points; + { + std::vector edge_paths = json_params["input_edges"]; + std::vector point_paths = json_params["input_points"]; + for (std::string& p : edge_paths) { + p = resolve_path(root, p).string(); + } + for (std::string& p : point_paths) { + p = resolve_path(root, p).string(); + } + for (size_t file = 0; file < edge_paths.size(); ++file) { + const std::string& path = edge_paths[file]; + MatrixXd Ve; + MatrixXi Ee; + io::read_edge_mesh(path, Ve, Ee); + logger().info("Read feature edge mesh {}: #V = {}, #E = {}", path, Ve.rows(), Ee.rows()); + const size_t base = feature_edge_vertices.size(); + for (int i = 0; i < Ve.rows(); ++i) { + feature_edge_vertices.emplace_back(Ve.row(i)); + } + for (int i = 0; i < Ee.rows(); ++i) { + feature_edges.push_back({{base + size_t(Ee(i, 0)), base + size_t(Ee(i, 1))}}); + } + // A vertex of an edge file with no incident edge is a free point, same rule as + // triwild's 2D inputs. + std::vector valence(Ve.rows(), 0); + for (int i = 0; i < Ee.rows(); ++i) { + ++valence[Ee(i, 0)]; + ++valence[Ee(i, 1)]; + } + for (int v = 0; v < Ve.rows(); ++v) { + if (valence[v] == 0) { + feature_points.emplace_back(Ve.row(v)); + } + } + } + for (const std::string& path : point_paths) { + MatrixXd Vp; + MatrixXi Ep; + io::read_edge_mesh(path, Vp, Ep); + if (Ep.rows() > 0) { + logger().warn( + "input_points file {} has {} edges; only its {} vertices are used", + path, + Ep.rows(), + Vp.rows()); + } + logger().info("Read feature point file {}: #P = {}", path, Vp.rows()); + for (int i = 0; i < Vp.rows(); ++i) { + feature_points.emplace_back(Vp.row(i)); + } + } + for (const Vector3d& p : feature_edge_vertices) { + box_minmax.first = box_minmax.first.cwiseMin(p); + box_minmax.second = box_minmax.second.cwiseMax(p); + } + for (const Vector3d& p : feature_points) { + box_minmax.first = box_minmax.first.cwiseMin(p); + box_minmax.second = box_minmax.second.cwiseMax(p); + } + } + + // The 0-dimensional features to anchor: every input point, and the feature network's + // endpoints (valence 1) -- junctions (valence >= 3) too when allow_junction_cleanup is + // off. Valence 2 is a curve interior: never anchored, the tube handles it. + // Anchors are REGISTERED unconditionally; preserve_feature_points gates only the + // collapse and smoothing policy. Same split as triwild, and it keeps the retention + // audit honest when the guard is off -- it reports 0/N instead of nothing to check. + std::vector feature_anchors; + { + feature_anchors = feature_points; + std::vector valence(feature_edge_vertices.size(), 0); + for (const auto& e : feature_edges) { + ++valence[e[0]]; + ++valence[e[1]]; + } + size_t n_endpoints = 0, n_junctions = 0; + for (size_t v = 0; v < feature_edge_vertices.size(); ++v) { + if (valence[v] == 0 || valence[v] == 2) { + continue; + } + if (valence[v] == 1) { + ++n_endpoints; + } else { + ++n_junctions; + if (params.allow_junction_cleanup) { + continue; + } + } + feature_anchors.push_back(feature_edge_vertices[v]); + } + if (n_endpoints + n_junctions > 0) { + logger().info( + "feature anchors: {} curve endpoints, {} junctions; anchoring {} points in " + "total (junction cleanup {})", + n_endpoints, + n_junctions, + feature_anchors.size(), + params.allow_junction_cleanup ? "on" : "off"); + } + } + t_load = phase_timer.getElapsedTime(); phase_timer.start(); // surface simplification begins @@ -438,6 +628,7 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) // Exact arrangement of the simplified surface against a Delaunay background // mesh; the remesher's own tets are used directly, so no Steiner points. + utils::EmbedFeaturesResult features_out; mesh.insertion_by_volumeremesher( vsimp, fsimp, @@ -445,7 +636,29 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) facets, is_v_on_input, tets, - tet_face_on_input_surface); + tet_face_on_input_surface, + feature_edge_vertices, + feature_edges, + feature_points, + &features_out); + + if (!feature_edges.empty() || !feature_points.empty()) { + size_t tiling_edges = 0; + for (const auto& t : features_out.edge_tiling) { + tiling_edges += t.size(); + } + size_t points_found = 0; + for (const int64_t v : features_out.point_vertex) { + points_found += v >= 0 ? 1 : 0; + } + logger().info( + "features after insertion: {} edges tiled by {} tet edges, {}/{} points are " + "output vertices", + features_out.edge_tiling.size(), + tiling_edges, + points_found, + features_out.point_vertex.size()); + } logger().info("=== finished insertion"); @@ -454,12 +667,31 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) wmtk::set_preallocation_factor_from_json(mesh_new, json_params); mesh_new.m_input_names = json_params["input_names"].get>(); + const bool has_features = !feature_edges.empty() || !feature_points.empty(); + // The feature-curve tube: what the collapse guard checks tagged edges against, and what + // curve-vertex smoothing will be pulled toward. Around the ORIGINAL input feature edges, + // like every other envelope. Follows the surface envelope's choice of predicate. + if (!feature_edges.empty()) { + std::vector fe_env(feature_edges.size()); + for (size_t i = 0; i < feature_edges.size(); ++i) { + fe_env[i] = Eigen::Vector2i(int(feature_edges[i][0]), int(feature_edges[i][1])); + } + mesh_new.m_feature_envelope = std::make_shared(!use_sample_envelope); + mesh_new.m_feature_envelope->init( + feature_edge_vertices, + fe_env, + params.epsr * params.diag_l * params.feature_envelope_ratio); + } + // The anchor ball has the tube's radius; input_points-only runs need it too. + mesh_new.m_feature_eps = params.epsr * params.diag_l * params.feature_envelope_ratio; mesh_new.init_from_Volumeremesher( v_rational, facets, is_v_on_input, tets, - tet_face_on_input_surface); + tet_face_on_input_surface, + has_features ? &features_out : nullptr, + has_features && !feature_anchors.empty() ? &feature_anchors : nullptr); double insertion_time = insertion_timer.getElapsedTime(); @@ -523,6 +755,35 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) // phase, so skip_winding_number lets a caller that does not filter and does not need // the annotations opt out of them. When filtering is requested the flag is ignored // (the winding number is needed), with a warning. + // Feature collections, taken BEFORE any filter deletes tets: explicit features are + // user input, and a filter discarding the region they live in must not erase them from + // the feature outputs and audits (they do disappear from the tet mesh itself). Same + // rule in 2D. With filter='none' nothing is deleted and this equals the final state. + std::vector> hybrid_sheet_tris; // hybrid filter only + std::vector> hybrid_curve_segs; + std::vector hybrid_anchor_pts; + std::pair hybrid_retention{0, 0}; + double hybrid_retention_worst = 0; + const bool features_present = !feature_edges.empty() || !feature_points.empty(); + if (features_present) { + for (const auto& e : mesh_new.get_edges()) { + if (!mesh_new.m_feature_edge_attribute[e.eid(mesh_new)].m_is_feature_edge) { + continue; + } + hybrid_curve_segs.push_back( + {{mesh_new.m_vertex_attribute[e.vid(mesh_new)].m_posf, + mesh_new.m_vertex_attribute[e.switch_vertex(mesh_new).vid(mesh_new)] + .m_posf}}); + } + for (const auto& v : mesh_new.get_vertices()) { + const size_t vid = v.vid(mesh_new); + if (mesh_new.m_vertex_extra[vid].m_feature_point_id != TetWildMesh::NO_FEATURE) { + hybrid_anchor_pts.push_back(mesh_new.m_vertex_attribute[vid].m_posf); + } + } + hybrid_retention = mesh_new.feature_retention(&hybrid_retention_worst); + } + const bool skip_winding = json_params["skip_winding_number"] && filter_option == "none"; if (json_params["skip_winding_number"] && filter_option != "none") { logger().warn( @@ -566,6 +827,44 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) const int num_parts = mesh_new.flood_fill(); logger().info("flood fill parts {}", num_parts); mesh_new.filter_with_flood_fill(); + } else if (filter_option == "hybrid") { + // Collect everything that lives on scaffolding tets BEFORE deleting them: the + // surface-role sheets (tracked faces none of whose incident tets stay), every + // tagged feature curve, the anchor positions, and the pre-filter retention -- the + // audit that means something is the one taken while the features still exist in + // the tet mesh. + std::vector keep(mesh_new.tet_capacity(), false); + for (const auto& t : mesh_new.get_tets()) { + const size_t tid = t.tid(mesh_new); + const auto& wn = mesh_new.m_tet_attribute[tid].m_winding_number_per_input; + for (const size_t k : volume_input_ids) { + if (k < wn.size() && wn[k] > 0.5) { + keep[tid] = true; + break; + } + } + } + for (const auto& f : mesh_new.get_faces()) { + if (!mesh_new.m_face_attribute[f.fid(mesh_new)].m_is_surface_fs) { + continue; + } + bool any_kept = keep[f.tid(mesh_new)]; + const auto oppo = f.switch_tetrahedron(mesh_new); + if (oppo.has_value()) { + any_kept = any_kept || keep[(*oppo).tid(mesh_new)]; + } + if (any_kept) { + continue; // a volume boundary (or interior) face; the tets represent it + } + const size_t v1 = f.vid(mesh_new); + const size_t v2 = f.switch_vertex(mesh_new).vid(mesh_new); + const size_t v3 = f.switch_edge(mesh_new).switch_vertex(mesh_new).vid(mesh_new); + hybrid_sheet_tris.push_back( + {{mesh_new.m_vertex_attribute[v1].m_posf, + mesh_new.m_vertex_attribute[v2].m_posf, + mesh_new.m_vertex_attribute[v3].m_posf}}); + } + mesh_new.filter_with_roles(volume_input_ids); } else if (filter_option != "none") { logger().error("Unknown filter option '{}'. No filtering performed.", filter_option); } @@ -669,6 +968,7 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) double hausdorff_distance = -1; // d(output -> input), the invariant double coverage_distance = -1; // d(input -> output), diagnostic only std::vector ecs_output; + std::vector ecs_curves_in, ecs_curves_out; { Eigen::MatrixXd V(verts.size(), 3); for (int i = 0; i < verts.size(); ++i) { @@ -731,6 +1031,65 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) coverage_distance); } + // Feature-curve deviation, the 1D counterpart of the surface block above: + // containment (tagged edges near the input curves -- the invariant the tube veto + // enforces) and coverage (input curves near the tagged edges -- diagnostic, the + // direction nothing enforces). Sampled point-to-segment on both sides; feature + // networks are small next to surfaces, so brute force over the segments is fine. + if (json_params["DEBUG_hausdorff"] && !feature_edges.empty()) { + // Collected pre-filter: under a filter the mesh no longer carries them. + const std::vector>& tagged = hybrid_curve_segs; + const auto point_to_segments = + [](const Vector3d& p, const std::vector>& segs) { + double best = std::numeric_limits::infinity(); + for (const auto& s : segs) { + const Vector3d d = s[1] - s[0]; + const double dd = d.squaredNorm(); + double t = dd > 0 ? (p - s[0]).dot(d) / dd : 0.0; + t = std::clamp(t, 0.0, 1.0); + best = std::min(best, (p - (s[0] + t * d)).squaredNorm()); + } + return best; + }; + std::vector> input_segs(feature_edges.size()); + for (size_t i = 0; i < feature_edges.size(); ++i) { + input_segs[i] = { + {feature_edge_vertices[feature_edges[i][0]], + feature_edge_vertices[feature_edges[i][1]]}}; + } + const auto sweep = [&](const std::vector>& from, + const std::vector>& to) { + double worst = -1; + for (const auto& s : from) { + const int n = 32; + for (int i = 0; i <= n; ++i) { + const Vector3d p = s[0] + (double(i) / n) * (s[1] - s[0]); + worst = std::max(worst, point_to_segments(p, to)); + } + } + return worst < 0 ? worst : std::sqrt(worst); + }; + const double feat_containment = tagged.empty() ? -1 : sweep(tagged, input_segs); + const double feat_coverage = tagged.empty() + ? std::numeric_limits::infinity() + : sweep(input_segs, tagged); + const double feat_eps = + params.epsr * params.diag_l * params.feature_envelope_ratio; + logger().info( + "feature deviation: {} tagged edges | containment d(tagged->input) = {:.4} | " + "tube = {:.4}", + tagged.size(), + feat_containment, + feat_eps); + if (feat_containment > feat_eps) { + logger().warn("Tagged feature edges left the tube; the veto was violated."); + } + logger().info( + "feature deviation: coverage d(input->tagged) = {:.4} (diagnostic; large " + "means part of an input curve is no longer represented)", + feat_coverage); + } + // The Euler-characteristic check is a topology sanity check. It is expensive on // meshes with many components (tens of seconds), so it is off by default and only // computed when explicitly requested (DEBUG_euler) or when it is actually needed @@ -743,6 +1102,53 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) logger().warn("Output topology is not the same as the input topology!"); } } + + // The same check for the feature curves: input network vs the tagged output edges. + // Always computed when features exist (the networks are tiny next to the surface); + // the preserve_topology throw is with the other throws below. + if (!feature_edges.empty()) { + std::vector> tagged_pairs; + std::map, size_t> vid_of; + for (const auto& seg : hybrid_curve_segs) { + std::array pair; + for (int j = 0; j < 2; ++j) { + const std::array key = {{seg[j][0], seg[j][1], seg[j][2]}}; + pair[j] = vid_of.emplace(key, vid_of.size()).first->second; + } + tagged_pairs.push_back(pair); + } + ecs_curves_in = curve_euler_characteristics(feature_edges); + ecs_curves_out = curve_euler_characteristics(tagged_pairs); + logger().info( + "Euler characteristic, feature curves: input {} | tagged output {}", + ecs_curves_in, + ecs_curves_out); + if (ecs_curves_in != ecs_curves_out) { + logger().warn("Feature-curve topology is not the same as the input's!"); + } + } + + // The anchor invariant, measured on the finished mesh. + if (json_params["DEBUG_feature_retention"]) { + double worst_ratio = 0; + // Measured pre-filter: the anchors' vertices may have been deleted with a + // discarded region, but the features were preserved up to that point and + // survive in the feature outputs. + auto [kept, total] = hybrid_retention; + worst_ratio = hybrid_retention_worst; + if (total > 0) { + if (kept == total) { + logger().info("feature points retained: {}/{}", kept, total); + } else { + logger().warn( + "feature points retained: {}/{} -- the worst is {:.2f} x eps from " + "the nearest vertex", + kept, + total, + worst_ratio); + } + } + } } /////////output @@ -810,6 +1216,9 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) if (params.preserve_topology && ecs_input != ecs_output) { log_and_throw_error("Input topology was not preserved."); } + if (params.preserve_topology && ecs_curves_in != ecs_curves_out) { + log_and_throw_error("Feature-curve topology was not preserved."); + } } @@ -819,8 +1228,47 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) mesh_new.output_mesh(output_path + "_final.msh"); + if (filter_option == "hybrid") { + mesh_new.output_hybrid_mesh( + output_path + "_hybrid.msh", + hybrid_sheet_tris, + hybrid_curve_segs, + hybrid_anchor_pts); + } + igl::write_triangle_mesh(output_path + "_surface.obj", matV, matF); + // The tracked feature curves, as an edge mesh -- the 1D counterpart of _surface.obj. + // Only written when features exist, so featureless runs are byte-identical. + if (features_present) { + // From the pre-filter collections: a filter must not erase user-supplied features + // from the feature outputs (they do leave the tet mesh with their region). + std::ofstream fout(output_path + "_features.obj"); + std::map, size_t> vid_of; + const auto obj_vertex = [&](const Vector3d& p) { + const std::array key = {{p[0], p[1], p[2]}}; + const auto [it, inserted] = vid_of.emplace(key, vid_of.size() + 1); + if (inserted) { + fout << "v " << p[0] << " " << p[1] << " " << p[2] << "\n"; + } + return it->second; + }; + std::vector> obj_edges; + for (const auto& seg : hybrid_curve_segs) { + obj_edges.push_back({{obj_vertex(seg[0]), obj_vertex(seg[1])}}); + } + std::vector obj_points; + for (const Vector3d& p : hybrid_anchor_pts) { + obj_points.push_back(obj_vertex(p)); + } + for (const auto& e : obj_edges) { + fout << "l " << e[0] << " " << e[1] << "\n"; + } + for (const size_t pid : obj_points) { + fout << "p " << pid << "\n"; + } + } + wmtk::logger().info("======= finish ========="); return mesh_new.export_mesh_data(); diff --git a/components/tetwild/wmtk/components/tetwild/tetwild_spec.json b/components/tetwild/wmtk/components/tetwild/tetwild_spec.json index 822cfd21aa..4e467fa5ff 100644 --- a/components/tetwild/wmtk/components/tetwild/tetwild_spec.json +++ b/components/tetwild/wmtk/components/tetwild/tetwild_spec.json @@ -6,6 +6,9 @@ "required": ["application", "input"], "optional": [ "output", + "input_edges", + "input_points", + "input_roles", "input_names", "input_dir", "write_vtu", @@ -16,6 +19,10 @@ "simplify_envelope_ratio", "optimize_envelope_around_simplified", "order2_envelope_ratio", + "feature_envelope_ratio", + "preserve_feature_points", + "allow_junction_cleanup", + "DEBUG_feature_retention", "use_sample_envelope", "use_legacy_code", "num_threads", @@ -81,6 +88,64 @@ "type": "string", "doc": "Triangular input mesh." }, + { + "pointer": "/input_edges", + "type": "list", + "default": [], + "doc": "List of input feature EDGE meshes: curve networks forced into the tetrahedralization exactly -- every input edge is represented by a chain of tet edges tiling it, endpoint to endpoint. A vertex of one of these files with no incident edge is treated as a feature point (see /input_points). The edges and points, plus the forcing-triangle apexes the insertion builds for them (offset by up to 0.1 x the edge length), must lie inside the background box; the input bounding box is grown over them automatically. NOTE: inserted and tracked only, for now -- the optimization phase does not yet preserve them." + }, + { + "pointer": "/input_edges/*", + "type": "string", + "doc": "One input feature edge mesh." + }, + { + "pointer": "/input_points", + "type": "list", + "default": [], + "doc": "List of input feature POINT files: every vertex of each file is forced into the tetrahedralization as a tet vertex, exactly (edges in these files are ignored, with a warning). NOTE: inserted and tracked only, for now -- the optimization phase does not yet preserve them." + }, + { + "pointer": "/input_points/*", + "type": "string", + "doc": "One input feature point file (only the vertices are used)." + }, + { + "pointer": "/input_roles", + "type": "list", + "default": [], + "doc": "For filter='hybrid': one role per /input entry, 'volume' (tetrahedralize its inside, by per-input winding number) or 'surface' (keep as an embedded triangle sheet in the hybrid output). Empty means every input is 'volume'. Curves and points declare their dimension by arriving through /input_edges and /input_points. Ignored (with a warning) by every other filter." + }, + { + "pointer": "/input_roles/*", + "type": "string", + "options": ["volume", "surface"], + "doc": "Role of one input mesh." + }, + { + "pointer": "/preserve_feature_points", + "type": "bool", + "default": true, + "doc": "Keep the 0-dimensional features -- every /input_points point, and the endpoints of open feature curves (junctions too, when allow_junction_cleanup is off) -- within the feature tube radius of where the input put them. Some output vertex stays within that ball of each anchor, collapses that would break this are refused, and smoothing may move an anchor-carrying vertex only inside its ball. Without this an open feature curve erodes from its own tips: every eroding collapse keeps the surviving chain inside the tube, so no containment test can see it. Mirrors triwild's parameter of the same name." + }, + { + "pointer": "/allow_junction_cleanup", + "type": "bool", + "default": true, + "doc": "Anchor only the endpoints of open feature curves, letting junctions (valence >= 3 in the feature network) merge and move within the tube. The erosion the anchor exists for is an ENDPOINT property -- a curve eats its own tip -- and does not apply to a junction, where every curve through it stays constrained. Mirrors triwild." + }, + { + "pointer": "/DEBUG_feature_retention", + "type": "bool", + "default": false, + "doc": "Sanity check: count how many feature anchors still have a mesh vertex within the anchor ball, measured BEFORE any /filter discards regions (a feature preserved through the optimization is never reported lost because extraction dropped its region; the feature outputs keep it either way), and report the worst miss in multiples of the ball radius. The invariant preserve_feature_points maintains, measured rather than assumed." + }, + { + "pointer": "/feature_envelope_ratio", + "type": "float", + "default": 0.5, + "doc": "Envelope thickness for the feature-curve tube -- the region a tet edge tagged as tiling an input feature edge must stay inside -- as a fraction of the surface envelope's eps. Kept below 1 for the same measured reason as order2_envelope_ratio: on a curve a wider envelope relieves no blockage, it is only freedom for the curve to wander, i.e. more geometry to resolve at the same final quality." + }, { "pointer": "/output", "type": "string", @@ -180,9 +245,9 @@ { "pointer": "/filter", "type": "string", - "options": ["flood", "input", "tracked", "none"], + "options": ["flood", "input", "tracked", "hybrid", "none"], "default": "none", - "doc": "Remove the outside region based on different criteria. 'flood': flood fill. 'input': winding number w.r.t. the input. 'tracked': winding number w.r.t. the tracked surface. 'none': Do not filter. Flood fill only works if the input is closed. Otherwise, it results in an empty mesh. Filtering w.r.t. the input might cause wrinkles along the surface as some tets might be falsely tagged. Filtering w.r.t. the tracked surface can lead to missing pieces if the input consists of multiple components." + "doc": "Remove the outside region based on different criteria. 'flood': flood fill. 'input': winding number w.r.t. the input. 'tracked': winding number w.r.t. the tracked surface. 'hybrid': mixed-dimensional output driven by /input_roles -- tets are kept only inside volume-role inputs, surface-role inputs survive as an embedded triangle sheet and the feature curves as an edge chain, all written together to _hybrid.msh (the plain _final.msh keeps only the volume tets). 'none': Do not filter. Flood fill only works if the input is closed. Otherwise, it results in an empty mesh. Filtering w.r.t. the input might cause wrinkles along the surface as some tets might be falsely tagged. Filtering w.r.t. the tracked surface can lead to missing pieces if the input consists of multiple components." }, { "pointer": "/skip_winding_number", From 62683addb491cd2a5ca7d0e79a9d4bc590027bb4 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 18 Aug 2026 13:34:26 -0400 Subject: [PATCH 6/7] clang-format the feature-preservation changes Co-Authored-By: Claude Fable 5 --- .../wmtk/components/tetwild/TetWildMesh.cpp | 4 +- .../wmtk/components/tetwild/TetWildMesh.h | 2 +- .../tetwild/VolumemesherInsertion.cpp | 5 ++- .../tetwild/tests/test_feature_tags.cpp | 17 ++++----- .../tetwild/tests/test_insertion.cpp | 5 ++- .../wmtk/components/tetwild/tetwild.cpp | 38 +++++++++---------- .../triwild/tests/test_free_points.cpp | 2 +- src/wmtk/TetMesh.h | 5 +-- src/wmtk/TetOptimizerMesh.cpp | 5 ++- src/wmtk/TetOptimizerMeshCollapse.cpp | 5 ++- src/wmtk/TetOptimizerMeshSplit.cpp | 3 +- src/wmtk/optimization/SmoothVertex.hpp | 3 +- src/wmtk/utils/EmbedTriangles.cpp | 3 +- src/wmtk/utils/io.hpp | 4 +- 14 files changed, 49 insertions(+), 52 deletions(-) diff --git a/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp b/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp index deb5596d3d..54719311b7 100644 --- a/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp +++ b/components/tetwild/wmtk/components/tetwild/TetWildMesh.cpp @@ -743,9 +743,7 @@ void TetWildMesh::output_hybrid_mesh( if (!anchor_points.empty()) { msh.add_point_vertices(anchor_points.size(), [&](size_t k) { return anchor_points[k]; }); - msh.add_points(anchor_points.size(), [&](size_t k) { - return std::array{{k}}; - }); + msh.add_points(anchor_points.size(), [&](size_t k) { return std::array{{k}}; }); } msh.save(file, /*binary=*/true); diff --git a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h index 22731dfc71..99e34a97ca 100644 --- a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h +++ b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h @@ -23,10 +23,10 @@ #include #include -#include #include #include #include +#include #include #include diff --git a/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp b/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp index 6d56ed7123..fb9e5d1ae9 100644 --- a/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp +++ b/components/tetwild/wmtk/components/tetwild/VolumemesherInsertion.cpp @@ -152,7 +152,10 @@ void TetWildMesh::insertion_by_volumeremesher( }; vr_features.edge_vrt_coord.reserve(3 * feature_edge_vertices.size()); for (const Vector3d& p : feature_edge_vertices) { - vr_features.edge_vrt_coord.insert(vr_features.edge_vrt_coord.end(), p.data(), p.data() + 3); + vr_features.edge_vrt_coord.insert( + vr_features.edge_vrt_coord.end(), + p.data(), + p.data() + 3); } vr_features.edge_indices.reserve(2 * feature_edges.size()); for (const auto& e : feature_edges) { diff --git a/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp b/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp index e6f5b1fe90..68ed52ef86 100644 --- a/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp +++ b/components/tetwild/wmtk/components/tetwild/tests/test_feature_tags.cpp @@ -54,9 +54,8 @@ FeatureScene build_scene() VF_to_vectors(V, F, vertices, faces); sc.params.init(vertices, faces); - sc.surf = std::make_shared( - vertices, - 0); + sc.surf = + std::make_shared(vertices, 0); { std::vector frozen_verts; sc.surf->create_mesh(vertices.size(), faces, frozen_verts, 0.1); @@ -130,11 +129,8 @@ std::vector> tagged_edges(TetWildMesh& m) } /// Exact parameter of vertex `vid` on segment (A,B), or nullopt if not exactly on it. -std::optional param_on_segment( - TetWildMesh& m, - const size_t vid, - const Vector3d& A, - const Vector3d& B) +std::optional +param_on_segment(TetWildMesh& m, const size_t vid, const Vector3d& A, const Vector3d& B) { const Vector3r a{Rational(A[0]), Rational(A[1]), Rational(A[2])}; const Vector3r d{Rational(B[0] - A[0]), Rational(B[1] - A[1]), Rational(B[2] - A[2])}; @@ -180,8 +176,9 @@ void require_tags_in_tube(FeatureScene& sc) { TetWildMesh& m = *sc.mesh; for (const auto& e : tagged_edges(m)) { - REQUIRE(!sc.mesh->m_feature_envelope->is_outside(std::array{ - {m.m_vertex_attribute[e[0]].m_posf, m.m_vertex_attribute[e[1]].m_posf}})); + REQUIRE(!sc.mesh->m_feature_envelope->is_outside( + std::array{ + {m.m_vertex_attribute[e[0]].m_posf, m.m_vertex_attribute[e[1]].m_posf}})); REQUIRE(m.m_vertex_extra[e[0]].m_is_on_feature_curve); REQUIRE(m.m_vertex_extra[e[1]].m_is_on_feature_curve); } diff --git a/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp b/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp index 0b015c11ce..180f2bef2a 100644 --- a/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp +++ b/components/tetwild/wmtk/components/tetwild/tests/test_insertion.cpp @@ -295,7 +295,10 @@ void require_exact_tiling( { REQUIRE(!tiling.empty()); const Vector3r a{Rational(seg_a[0]), Rational(seg_a[1]), Rational(seg_a[2])}; - const Vector3r d{Rational(seg_b[0] - seg_a[0]), Rational(seg_b[1] - seg_a[1]), Rational(seg_b[2] - seg_a[2])}; + const Vector3r d{ + Rational(seg_b[0] - seg_a[0]), + Rational(seg_b[1] - seg_a[1]), + Rational(seg_b[2] - seg_a[2])}; const Rational dd = d.dot(d); // Parameter of a vertex along AB, after requiring it exactly on the line. diff --git a/components/tetwild/wmtk/components/tetwild/tetwild.cpp b/components/tetwild/wmtk/components/tetwild/tetwild.cpp index 594613cca0..4b9b69aa39 100644 --- a/components/tetwild/wmtk/components/tetwild/tetwild.cpp +++ b/components/tetwild/wmtk/components/tetwild/tetwild.cpp @@ -276,7 +276,8 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) MatrixXd Ve; MatrixXi Ee; io::read_edge_mesh(path, Ve, Ee); - logger().info("Read feature edge mesh {}: #V = {}, #E = {}", path, Ve.rows(), Ee.rows()); + logger() + .info("Read feature edge mesh {}: #V = {}, #E = {}", path, Ve.rows(), Ee.rows()); const size_t base = feature_edge_vertices.size(); for (int i = 0; i < Ve.rows(); ++i) { feature_edge_vertices.emplace_back(Ve.row(i)); @@ -772,8 +773,7 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) } hybrid_curve_segs.push_back( {{mesh_new.m_vertex_attribute[e.vid(mesh_new)].m_posf, - mesh_new.m_vertex_attribute[e.switch_vertex(mesh_new).vid(mesh_new)] - .m_posf}}); + mesh_new.m_vertex_attribute[e.switch_vertex(mesh_new).vid(mesh_new)].m_posf}}); } for (const auto& v : mesh_new.get_vertices()) { const size_t vid = v.vid(mesh_new); @@ -1039,18 +1039,18 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) if (json_params["DEBUG_hausdorff"] && !feature_edges.empty()) { // Collected pre-filter: under a filter the mesh no longer carries them. const std::vector>& tagged = hybrid_curve_segs; - const auto point_to_segments = - [](const Vector3d& p, const std::vector>& segs) { - double best = std::numeric_limits::infinity(); - for (const auto& s : segs) { - const Vector3d d = s[1] - s[0]; - const double dd = d.squaredNorm(); - double t = dd > 0 ? (p - s[0]).dot(d) / dd : 0.0; - t = std::clamp(t, 0.0, 1.0); - best = std::min(best, (p - (s[0] + t * d)).squaredNorm()); - } - return best; - }; + const auto point_to_segments = [](const Vector3d& p, + const std::vector>& segs) { + double best = std::numeric_limits::infinity(); + for (const auto& s : segs) { + const Vector3d d = s[1] - s[0]; + const double dd = d.squaredNorm(); + double t = dd > 0 ? (p - s[0]).dot(d) / dd : 0.0; + t = std::clamp(t, 0.0, 1.0); + best = std::min(best, (p - (s[0] + t * d)).squaredNorm()); + } + return best; + }; std::vector> input_segs(feature_edges.size()); for (size_t i = 0; i < feature_edges.size(); ++i) { input_segs[i] = { @@ -1070,11 +1070,9 @@ TetWildMesh::ExportStruct tetwild_with_export(nlohmann::json json_params) return worst < 0 ? worst : std::sqrt(worst); }; const double feat_containment = tagged.empty() ? -1 : sweep(tagged, input_segs); - const double feat_coverage = tagged.empty() - ? std::numeric_limits::infinity() - : sweep(input_segs, tagged); - const double feat_eps = - params.epsr * params.diag_l * params.feature_envelope_ratio; + const double feat_coverage = tagged.empty() ? std::numeric_limits::infinity() + : sweep(input_segs, tagged); + const double feat_eps = params.epsr * params.diag_l * params.feature_envelope_ratio; logger().info( "feature deviation: {} tagged edges | containment d(tagged->input) = {:.4} | " "tube = {:.4}", diff --git a/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp b/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp index cfcca97aa4..d7b23aea6e 100644 --- a/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp +++ b/components/triwild/wmtk/components/triwild/tests/test_free_points.cpp @@ -1,5 +1,5 @@ -#include #include +#include #include #include diff --git a/src/wmtk/TetMesh.h b/src/wmtk/TetMesh.h index bd5c8c2d2a..4aa77ac4c5 100644 --- a/src/wmtk/TetMesh.h +++ b/src/wmtk/TetMesh.h @@ -1224,9 +1224,8 @@ class TetMesh auto common = set_intersection( m_vertex_connectivity[vids[0]].m_conn_tets, m_vertex_connectivity[vids[1]].m_conn_tets); - common = set_intersection( - common, - m_vertex_connectivity[vids[2]].m_conn_tets); + common = + set_intersection(common, m_vertex_connectivity[vids[2]].m_conn_tets); size_t other = std::numeric_limits::max(); for (const size_t t : common) { if (t != tid) { diff --git a/src/wmtk/TetOptimizerMesh.cpp b/src/wmtk/TetOptimizerMesh.cpp index e483535974..14484eb684 100644 --- a/src/wmtk/TetOptimizerMesh.cpp +++ b/src/wmtk/TetOptimizerMesh.cpp @@ -761,8 +761,9 @@ bool TetOptimizerMesh::feature_edges_at_vertex_inside(const size_t vid) const if (!e.is_valid(*this) || !m_feature_edge_attribute[e.eid(*this)].m_is_feature_edge) { continue; } - if (m_feature_envelope->is_outside(std::array{ - {m_vertex_attribute[vid].m_posf, m_vertex_attribute[u].m_posf}})) { + if (m_feature_envelope->is_outside( + std::array{ + {m_vertex_attribute[vid].m_posf, m_vertex_attribute[u].m_posf}})) { return false; } } diff --git a/src/wmtk/TetOptimizerMeshCollapse.cpp b/src/wmtk/TetOptimizerMeshCollapse.cpp index 4659fcd699..5ee660af3c 100644 --- a/src/wmtk/TetOptimizerMeshCollapse.cpp +++ b/src/wmtk/TetOptimizerMeshCollapse.cpp @@ -257,8 +257,9 @@ bool TetOptimizerMesh::collapse_edge_before(const Tuple& loc) // input is an edg if (p[0] == p[1]) { continue; // the collapsed edge itself; it disappears } - if (m_feature_envelope->is_outside(std::array{ - {m_vertex_attribute[p[0]].m_posf, m_vertex_attribute[p[1]].m_posf}})) { + if (m_feature_envelope->is_outside( + std::array{ + {m_vertex_attribute[p[0]].m_posf, m_vertex_attribute[p[1]].m_posf}})) { return false; } } diff --git a/src/wmtk/TetOptimizerMeshSplit.cpp b/src/wmtk/TetOptimizerMeshSplit.cpp index 0055690ae7..e4fd3549ef 100644 --- a/src/wmtk/TetOptimizerMeshSplit.cpp +++ b/src/wmtk/TetOptimizerMeshSplit.cpp @@ -186,8 +186,7 @@ bool TetOptimizerMesh::split_edge_before(const Tuple& loc0) tids.push_back(t.tid(*this)); } feature_edges_cache(tids, cache.changed_edges); - cache.is_edge_on_feature = - m_feature_edge_attribute[loc0.eid(*this)].m_is_feature_edge; + cache.is_edge_on_feature = m_feature_edge_attribute[loc0.eid(*this)].m_is_feature_edge; } return split_before_cells(loc0, tets); diff --git a/src/wmtk/optimization/SmoothVertex.hpp b/src/wmtk/optimization/SmoothVertex.hpp index 0bf600b16f..8e2370fe82 100644 --- a/src/wmtk/optimization/SmoothVertex.hpp +++ b/src/wmtk/optimization/SmoothVertex.hpp @@ -230,8 +230,7 @@ bool smooth_vertex_3d( // the edge tags, so it cannot go stale. const bool on_feature_curve = m.m_track_feature_edges && m.vertex_has_feature_edge(vid); const std::shared_ptr pull_env = - (VA[vid].m_is_on_surface || on_feature_curve) ? m.smoothing_energy_envelope(vid) - : nullptr; + (VA[vid].m_is_on_surface || on_feature_curve) ? m.smoothing_energy_envelope(vid) : nullptr; if (pull_env && opts.smoothing_mode == SmoothVertexOptions::SmoothingMode::Projected) { // Smooth as if the vertex were interior, then walk back onto the input. diff --git a/src/wmtk/utils/EmbedTriangles.cpp b/src/wmtk/utils/EmbedTriangles.cpp index 13ba60c2ce..1e2abb2eb3 100644 --- a/src/wmtk/utils/EmbedTriangles.cpp +++ b/src/wmtk/utils/EmbedTriangles.cpp @@ -405,7 +405,8 @@ void embed_triangles_in_tets( const int64_t b = v_map[seg[2]]; if (a < 0 || b < 0) { log_and_throw_error( - "Feature edge {}: an output edge vertex was compacted away", e); + "Feature edge {}: an output edge vertex was compacted away", + e); } tiling.push_back({{size_t(a), size_t(b)}}); } diff --git a/src/wmtk/utils/io.hpp b/src/wmtk/utils/io.hpp index e5d4255e77..9041dc321a 100644 --- a/src/wmtk/utils/io.hpp +++ b/src/wmtk/utils/io.hpp @@ -442,9 +442,7 @@ class MshData template void add_simplex_elements(size_t num_elements, const Fn& get_element_cb) { - static_assert( - DIM >= 0 && DIM <= 3, - "Only 0,1,2,3D simplex elements are supported"); + static_assert(DIM >= 0 && DIM <= 3, "Only 0,1,2,3D simplex elements are supported"); if (num_elements == 0) return; if (m_spec.nodes.num_nodes == 0) { From c8ce8df5739c1127f36afa1c7ca61f9f83f607e7 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Wed, 19 Aug 2026 08:39:17 -0400 Subject: [PATCH 7/7] Stop open sheet boundaries from eroding: three fixes in the order-2 machinery An open sheet boundary contracted along itself (measured: a corner walked 0.5, 14% of a flat sheet's area gone) through three compounding holes: 1. The collapse-side order-2 chord check -- live since the 2025-10-10 open boundary fixes -- was parked by the collapse refactoring (cb9b81a66c), which could not call tetwild's is_open_boundary_edge from the shared base. Restored through the collapse_is_order_2_edge virtual that same refactoring introduced for the caching half. 2. collapse_after_vertex decided whether the survivor keeps its open boundary flag BEFORE the merged face tags were written, reading a half updated state and sporadically stripping genuine boundary vertices of their flag -- after which every boundary guard went blind to them. Moved after the attribute write-back. 3. Smoothing had no boundary-edge check at all: a boundary vertex is only PULLED toward the boundary, so each smooth can walk it along the curve -- legal per step, invisible to the one-sided containment -- and the coarsening composite's per-collapse ring smoothing amplifies the walk by thousands. Added boundary_edges_at_vertex_inside, the open-boundary counterpart of the feature check: every open-boundary edge at the vertex must stay inside the order-2 tube. Membership is decided by flags plus the exact incident-face count, NOT by the tube -- a tube-based membership self-disables the moment geometry leaves the tube. Repro: hybrid-data/feature_demo4 (cube + two open sheets + a wire). Sheet area with no boundary declarations: 2.470 -> 2.876 of 2.880, corners within 0.002. Behavior changes exactly on open-boundary models: in the integration + challenging suites every diffing config has open boundaries, every boundary-free config is byte-identical, all 31 challenging models still converge. Co-Authored-By: Claude Fable 5 --- .../wmtk/components/tetwild/TetWildMesh.h | 24 +++++++++++ src/wmtk/TetOptimizerMesh.h | 7 ++++ src/wmtk/TetOptimizerMeshCollapse.cpp | 42 ++++++++++++------- src/wmtk/optimization/SmoothVertex.hpp | 10 +++++ 4 files changed, 67 insertions(+), 16 deletions(-) diff --git a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h index 99e34a97ca..c86b0ffc2e 100644 --- a/components/tetwild/wmtk/components/tetwild/TetWildMesh.h +++ b/components/tetwild/wmtk/components/tetwild/TetWildMesh.h @@ -121,6 +121,30 @@ class TetWildMesh : public wmtk::TetOptimizerMesh { m_vertex_extra[vid].m_is_on_open_boundary = is_open_boundary; } + bool boundary_edges_at_vertex_inside(size_t vid) const override + { + if (!m_vertex_extra[vid].m_is_on_open_boundary || !m_order2_envelope) { + return true; + } + for (const size_t u : get_one_ring_vids_for_vertex(vid)) { + if (!m_vertex_extra[u].m_is_on_open_boundary) { + continue; + } + // Exact membership: an open boundary edge has exactly one incident tracked + // face. Flags alone would also catch interior chords between two boundary + // vertices (e.g. across a corner), which must not be tube-tested. + if (get_num_surface_faces_for_edge({{vid, u}}) != 1) { + continue; + } + if (m_order2_envelope->is_outside( + std::array{ + {m_vertex_attribute[vid].m_posf, m_vertex_attribute[u].m_posf}})) { + return false; + } + } + return true; + } + bool collapse_before_vertex(size_t v1, size_t v2, double edge_length) override { if (collapse_breaks_feature_point(v1, v2)) return false; diff --git a/src/wmtk/TetOptimizerMesh.h b/src/wmtk/TetOptimizerMesh.h index 3805198d9f..29c87cae7a 100644 --- a/src/wmtk/TetOptimizerMesh.h +++ b/src/wmtk/TetOptimizerMesh.h @@ -171,6 +171,13 @@ class TetOptimizerMesh : public wmtk::TetMesh, public wmtk::RationalPositions /// position already written into the vertex attribute. bool feature_edges_at_vertex_inside(size_t vid) const; + /// The derived-boundary counterpart: are all OPEN-BOUNDARY edges at `vid` inside the + /// order-2 tube, at the current positions? Default true; tetwild overrides with its + /// open-boundary flags and envelope. Membership is decided by flags + the topological + /// face count, NOT by the tube itself -- a tube-based membership test self-disables the + /// moment geometry leaves the tube, which is this defect's signature. + virtual bool boundary_edges_at_vertex_inside(size_t) const { return true; } + /// Per-vertex positional constraint, on top of the envelopes -- the 3D counterpart of /// TriOptimizerMesh's hook of the same name. An application uses this to pin a vertex to /// a 0-dimensional feature it stands for, within a ball. Default: no constraint. diff --git a/src/wmtk/TetOptimizerMeshCollapse.cpp b/src/wmtk/TetOptimizerMeshCollapse.cpp index 5ee660af3c..4723fbe0b0 100644 --- a/src/wmtk/TetOptimizerMeshCollapse.cpp +++ b/src/wmtk/TetOptimizerMeshCollapse.cpp @@ -426,24 +426,30 @@ bool TetOptimizerMesh::collapse_edge_after(const Tuple& loc) // return false; // } } - // for (const auto& vids : cache.boundary_edges) { - // if (!is_open_boundary_edge(vids)) { - // // edge was an open boundary before (that is why it got cached) but is not anymore - // // after collapse - // return false; - // } - //} + // The order-2 chord gate: every cached order-2 edge (v1,x), remapped to (v2,x) at + // caching time, must still BE order-2 substructure -- flags and containment in the + // order-2 envelope included. Without it a sheet erodes IN-PLANE from its open + // boundary: boundary-to-boundary collapses cut corners with triangles lying exactly + // on the sheet's plane, invisible to the one-sided surface envelope (measured: + // 14% of a flat sheet's area gone, a corner receded 150x the order-2 tube). The + // check was live from the 2025-10-10 open-boundary fixes until the collapse + // refactoring (cb9b81a66c) parked it -- the shared base could not call tetwild's + // is_open_boundary_edge -- while giving its caching half the virtual below. This is + // that predicate, restored through the same virtual. + for (const auto& vids : cache.boundary_edges) { + if (!collapse_is_order_2_edge(vids)) { + // The edge was order-2 before the collapse (that is why it was cached) but + // is not anymore -- the collapse moved it off its curve. + return false; + } + } } - // Must run HERE, before the attribute updates below -- not at the end of the function. - // tetwild's override asks is_vertex_on_boundary(v2), which reads BOTH - // m_vertex_attribute[..].m_is_on_surface and m_face_attribute[..].m_is_surface_fs - // (TetWildMesh.cpp), and the two blocks below overwrite exactly those: the vertex flag is - // OR-ed from v1 just after round(), and the cached face attributes are written onto the - // post-collapse faces. Asking afterwards is asking a different question, and it changes - // which vertices keep their open-boundary flag -- and therefore which later collapses are - // allowed. - collapse_after_vertex(v1_id, v2_id); + // (collapse_after_vertex moved below: its is_vertex_on_boundary(v2) test must see the + // POST-collapse face tags, which the loops below are what write. Reading the + // half-updated state cleared open-boundary flags on genuine sheet corners -- measured + // on feature_demo4: three corners lost their flags, every boundary guard went blind to + // them, and the coarsening pass collapsed them away, 0.5 deep.) //// update attrs // tet attr @@ -477,6 +483,10 @@ bool TetOptimizerMesh::collapse_edge_after(const Tuple& loc) } } + // Now that the vertex flag OR (above) and the face tags (the loop above) are the + // post-collapse truth, decide whether the survivor keeps its open-boundary flag. + collapse_after_vertex(v1_id, v2_id); + if (!m_coarsen_mode) { return true; } diff --git a/src/wmtk/optimization/SmoothVertex.hpp b/src/wmtk/optimization/SmoothVertex.hpp index 8e2370fe82..1b03241998 100644 --- a/src/wmtk/optimization/SmoothVertex.hpp +++ b/src/wmtk/optimization/SmoothVertex.hpp @@ -342,6 +342,16 @@ bool smooth_vertex_3d( return false; } + // Open-boundary containment, the same test for the DERIVED order-2 curves: every open + // boundary edge at this vertex must stay inside the order-2 tube. Without it smoothing + // walks boundary vertices along the boundary -- each step legal for the vertex, no check + // on the edges -- and the coarsening composite's mass smoothing contracts an open sheet + // boundary far past its tube (measured: a corner walked 0.5, 14% of the sheet gone). + if (!m.boundary_edges_at_vertex_inside(vid)) { + if (counters) ++counters->envelope; + return false; + } + // Containment: every surface triangle at this vertex must still be inside. Checked // against the containment envelope, which is not necessarily the one it was pulled to. const std::shared_ptr check_env =