diff --git a/cmake/wmtk_data.cmake b/cmake/wmtk_data.cmake index a76d3a1791..9a09ab0a59 100644 --- a/cmake/wmtk_data.cmake +++ b/cmake/wmtk_data.cmake @@ -16,7 +16,7 @@ ExternalProject_Add( SOURCE_DIR ${WMTK_DATA_ROOT} GIT_REPOSITORY https://github.com/wildmeshing/data2.git - GIT_TAG 36bb3968aa1e685d3da2e38b87bc6fa5ab328a13 + GIT_TAG fc5576bb4d6444776911f16b8b4bdde94c5c9b1a CONFIGURE_COMMAND "" BUILD_COMMAND "" diff --git a/components/CMakeLists.txt b/components/CMakeLists.txt index 1dc30bbe1c..7f9c67842b 100644 --- a/components/CMakeLists.txt +++ b/components/CMakeLists.txt @@ -23,6 +23,7 @@ add_subdirectory("shortest_edge_collapse/wmtk/components/shortest_edge_collapse" add_subdirectory("isotropic_remeshing/wmtk/components/isotropic_remeshing") add_subdirectory("qslim/wmtk/components/qslim") add_subdirectory("tetwild/wmtk/components/tetwild") +add_subdirectory("triwild/wmtk/components/triwild") add_subdirectory("simplicial_embedding/wmtk/components/simplicial_embedding") add_subdirectory("image_simulation/wmtk/components/image_simulation") add_subdirectory("c1_simplification/wmtk/components/c1_simplification") diff --git a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp index dcc5dbe40d..3cb1e2bba0 100644 --- a/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp @@ -1259,7 +1259,7 @@ bool EmbedSurface::embed_surface_tetgen() m_V_emb_r.resizeLike(m_V_emb); for (int i = 0; i < m_V_emb.rows(); ++i) { - m_V_emb_r.row(i) = to_rational(m_V_emb.row(i)); + m_V_emb_r.row(i) = to_rational(Vector3d(m_V_emb.row(i))); } // add tags diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 90f09fb60e..f4d707ea8e 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -1273,7 +1273,6 @@ bool ImageSimulationMeshTri::collapse_edge_before(const Tuple& loc) } bool ImageSimulationMeshTri::collapse_edge_after(const Tuple& loc) - { auto& VA = m_vertex_attribute; auto& cache = collapse_cache.local(); diff --git a/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp b/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp index 2daf564e7b..6970571027 100644 --- a/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/VolumemesherInsertion.cpp @@ -637,7 +637,7 @@ void ImageSimulationMesh::init_from_image( for (int i = 0; i < vert_capacity(); i++) { m_vertex_attribute[i].m_pos = V.row(i); - m_vertex_attribute[i].m_posf = to_double(V.row(i)); + m_vertex_attribute[i].m_posf = to_double(m_vertex_attribute[i].m_pos); } // sanity check @@ -704,8 +704,8 @@ void ImageSimulationMesh::init_from_image( m_face_attribute.m_attributes.resize(T.rows() * 4); for (int i = 0; i < vert_capacity(); i++) { - m_vertex_attribute[i].m_pos = to_rational(V.row(i)); m_vertex_attribute[i].m_posf = V.row(i); + m_vertex_attribute[i].m_pos = to_rational(m_vertex_attribute[i].m_posf); m_vertex_attribute[i].m_is_rounded = true; } diff --git a/components/triwild/wmtk/components/triwild/CMakeLists.txt b/components/triwild/wmtk/components/triwild/CMakeLists.txt new file mode 100644 index 0000000000..86d15731ce --- /dev/null +++ b/components/triwild/wmtk/components/triwild/CMakeLists.txt @@ -0,0 +1,45 @@ +# triwild + +set(COMPONENT_NAME triwild) +add_component(${COMPONENT_NAME}) + +if(NOT ${WMTK_ENABLE_COMPONENT_${COMPONENT_NAME}}) + return() +endif() + +set(SRC_FILES + TriWildMesh.cpp + TriWildMesh.h + + Parameters.h + + EdgeSplitting.cpp + EdgeCollapsing.cpp + EdgeSwapping.cpp + Smooth.cpp + + triwild.hpp + triwild.cpp + triwild_spec.hpp + + init_from_delaunay.hpp + init_from_delaunay.cpp +) + +target_sources(wmtk_${COMPONENT_NAME} PRIVATE ${SRC_FILES}) + +target_link_libraries(wmtk_${COMPONENT_NAME} PUBLIC + wmtk::toolkit + wmtk::data + igl::predicates + jse::jse +) + +###################### +## TESTS +###################### + +# if(${WMTK_COMPONENT_TESTS}) +# add_subdirectory(tests) +# endif() + diff --git a/components/triwild/wmtk/components/triwild/EdgeCollapsing.cpp b/components/triwild/wmtk/components/triwild/EdgeCollapsing.cpp new file mode 100644 index 0000000000..dfa1e818cd --- /dev/null +++ b/components/triwild/wmtk/components/triwild/EdgeCollapsing.cpp @@ -0,0 +1,300 @@ +#include "TriWildMesh.h" +#include "oneapi/tbb/concurrent_vector.h" +#include "wmtk/TriMesh.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace wmtk::components::triwild { + +void TriWildMesh::collapse_all_edges(bool is_limit_length) +{ + std::vector> all_ops; + + auto setup_and_execute = [&](auto& executor) { + executor.priority = [](const TriWildMesh& m, Op op, const Tuple& t) { + return -m.get_length2(t); + }; + executor.num_threads = NUM_THREADS; + executor.is_weight_up_to_date = [&](const TriWildMesh& m, + const std::tuple& ele) { + const auto& VA = m_vertex_attribute; + auto& [weight, op, tup] = ele; + const double length = m.get_length2(tup); + if (length != -weight) { + return false; + } + // + size_t v1_id = tup.vid(*this); + size_t v2_id = tup.switch_vertex(*this).vid(*this); + double sizing_ratio = (VA[v1_id].m_sizing_scalar + VA[v2_id].m_sizing_scalar) / 2; + if (is_limit_length && length > m_params.collapsing_l2 * sizing_ratio * sizing_ratio) + return false; + return true; + }; + + // Execute!! + do { + all_ops.clear(); + const auto all_edges = get_edges(); + logger().info("#E = {}", all_edges.size()); + for (const Tuple& loc : all_edges) { + // collect all edges. Filtering too long edges happens in `is_weight_up_to_date` + all_ops.emplace_back("edge_collapse", loc); + all_ops.emplace_back("edge_collapse", loc.switch_vertex(*this)); + } + executor(*this, all_ops); + logger().info( + "success: {}, failed: {}", + executor.get_cnt_success(), + executor.get_cnt_fail()); + } while (executor.get_cnt_success() > 0); + }; + + igl::Timer timer; + timer.start(); + if (NUM_THREADS > 0) { + auto executor = ExecutePass(ExecutionPolicy::kPartition); + executor.lock_vertices = [](TriWildMesh& m, const Tuple& e, int task_id) -> bool { + return m.try_set_edge_mutex_two_ring(e, task_id); + }; + setup_and_execute(executor); + logger().info("edge collapse time parallel: {:.4}s", timer.getElapsedTimeInSec()); + } else { + auto executor = ExecutePass(ExecutionPolicy::kSeq); + setup_and_execute(executor); + logger().info("edge collapse time serial: {:.4}s", timer.getElapsedTimeInSec()); + } +} + +bool TriWildMesh::collapse_edge_before(const Tuple& loc) // input is an edge +{ + const auto& VA = m_vertex_attribute; + auto& cache = collapse_cache.local(); + + cache.changed_edges.clear(); + cache.changed_fids.clear(); + cache.changed_energies.clear(); + cache.surface_edges.clear(); + + size_t v1_id = loc.vid(*this); + auto loc1 = switch_vertex(loc); + size_t v2_id = loc1.vid(*this); + + cache.v1_id = v1_id; + cache.v2_id = v2_id; + + cache.edge_length = (VA[v1_id].m_posf - VA[v2_id].m_posf).norm(); + + ///check if on bbox/surface/boundary + // bbox + if (!VA[v1_id].on_bbox_faces.empty()) { + if (VA[v2_id].on_bbox_faces.size() < VA[v1_id].on_bbox_faces.size()) { + return false; + } + for (int on_bbox : VA[v1_id].on_bbox_faces) + if (std::find( + VA[v2_id].on_bbox_faces.begin(), + VA[v2_id].on_bbox_faces.end(), + on_bbox) == VA[v2_id].on_bbox_faces.end()) { + return false; + } + } + + // surface + if (cache.edge_length > 0 && VA[v1_id].m_is_on_surface) { + if (VA[v1_id].m_is_rounded && !VA[v2_id].m_is_on_surface) { + return false; // do not collapse away from surface + } + + if (VA[v1_id].m_is_rounded && get_order_of_vertex(v1_id) > 1 && + get_order_of_vertex(v2_id) <= 1) { + /** + * In general, we don't want to collapse away from feature vertices. It is always fine + * to collapse into a feature vertex, though. However, we allow to feature vertices to + * be collapsed. + */ + return false; + } + + // if both vertices are on the surface, the collapsing edge should be inside the envelope + const size_t eid = loc.eid(*this); + if (VA[v2_id].m_is_on_surface && !m_edge_attribute.at(eid).m_is_surface_fs) { + const Vector2d& p1 = VA[v1_id].m_posf; + const Vector2d& p2 = VA[v2_id].m_posf; + if (m_envelope->is_outside(std::array{p1, p2})) { + return false; + } + } + } + + const auto& n1_locs = get_one_ring_fids_for_vertex(loc); + + cache.changed_fids.reserve(n1_locs.size()); + cache.max_energy = 0; + for (const size_t& tid : n1_locs) { + const double q = m_face_attribute.at(tid).m_quality; + cache.max_energy = std::max(cache.max_energy, q); + const auto vs = oriented_tri_vids(tid); + if (vs[0] != v2_id && vs[1] != v2_id && vs[2] != v2_id) { + cache.changed_fids.emplace_back(tid); + } + } + + // pre-compute after-collapse energies + cache.changed_energies.reserve(cache.changed_fids.size()); + for (const size_t tid : cache.changed_fids) { + std::array vs = oriented_tri_vids(tid); + for (size_t i = 0; i < 3; ++i) { + if (vs[i] == v1_id) { + vs[i] = v2_id; + break; + } + } + + if (is_inverted(vs)) { + return false; + } + double q = get_quality(vs); + // quality check only when v1 is rounded + if (VA[v1_id].m_is_rounded && q > m_params.stop_energy && q > cache.max_energy) { + return false; + } + cache.changed_energies.emplace_back(q); + } + assert(cache.changed_energies.size() == cache.changed_fids.size()); + + // + const auto& n12_locs = get_incident_fids_for_edge(loc); + for (const size_t& tid : n12_locs) { + auto vs = oriented_tri_vids(tid); + std::array e_vids = {{v1_id, 0}}; + int cnt = 1; + // get the vertex that is not v1/v2, i.e., the edge-link vertices. + for (int j = 0; j < 3; j++) { + if (vs[j] != v1_id && vs[j] != v2_id) { + e_vids[cnt] = vs[j]; + cnt++; + } + } + auto [_1, global_eid1] = tuple_from_edge(e_vids); + auto [_2, global_eid2] = tuple_from_edge({{v2_id, e_vids[1]}}); + auto e_attr = m_edge_attribute.at(global_eid1); + e_attr.merge(m_edge_attribute.at(global_eid2)); + cache.changed_edges.push_back(std::make_pair(e_attr, e_vids)); + } + + 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 + std::vector> fs; + for (const size_t& tid : n1_locs) { + const auto vs = oriented_tri_vids(tid); + + int j_v1 = -1; + auto skip = [&]() { + for (int j = 0; j < 3; j++) { + const size_t vid = vs[j]; + if (vid == v2_id) { + // ignore tets incident to the edge (v1,v2) + return true; // v1-v2 definitely not on surface. + } + if (vid == v1_id) j_v1 = j; + } + return false; + }; + if (skip()) continue; + + for (int k = 0; k < 2; k++) { + auto va = vs[(j_v1 + 1 + k) % 3]; + if (VA[va].m_is_on_surface) { + std::array f = {{v1_id, va}}; + const auto [f_tuple, fid] = tuple_from_edge(f); + if (!m_edge_attribute.at(fid).m_is_surface_fs) { + // check if this face is actually on the surface + continue; + } + std::sort(f.begin(), f.end()); + fs.push_back(f); + } + } + } + wmtk::vector_unique(fs); + + cache.surface_edges.reserve(fs.size()); + for (auto& f : fs) { + std::replace(f.begin(), f.end(), v1_id, v2_id); + cache.surface_edges.push_back(f); + } + } + + if (m_params.preserve_topology) { + const bool v1_surf = VA[v1_id].m_is_on_surface; + const bool v2_surf = VA[v2_id].m_is_on_surface; + const bool v1_bbox = !VA[v1_id].on_bbox_faces.empty(); + const bool v2_bbox = !VA[v2_id].on_bbox_faces.empty(); + + if ((v1_surf || v1_bbox) && (v2_surf || v2_bbox)) { + if (!substructure_link_condition(loc)) { + return false; + } + } + } + + return true; +} + +bool TriWildMesh::collapse_edge_after(const Tuple& loc) +{ + auto& VA = m_vertex_attribute; + auto& cache = collapse_cache.local(); + size_t v1_id = cache.v1_id; + size_t v2_id = cache.v2_id; + + if (!TriMesh::collapse_edge_after(loc)) { + return false; + } + + // surface + if (cache.edge_length > 0) { + for (auto& vids : cache.surface_edges) { + const Vector2d a = VA.at(vids[0]).m_posf; + const Vector2d b = VA.at(vids[1]).m_posf; + // surface envelope + bool is_out = m_envelope->is_outside(std::array{a, b}); + if (is_out) { + return false; + } + } + } + + //// update attrs + // tet attr + for (int i = 0; i < cache.changed_fids.size(); i++) { + m_face_attribute[cache.changed_fids[i]].m_quality = cache.changed_energies[i]; + } + // vertex attr + VA[v2_id].m_is_on_surface = VA.at(v1_id).m_is_on_surface || VA.at(v2_id).m_is_on_surface; + + // no need to update on_bbox_faces + // face attr + for (auto& info : cache.changed_edges) { + auto& f_attr = info.first; + auto& old_vids = info.second; + // + auto [_, global_fid] = tuple_from_edge({{v2_id, old_vids[1]}}); + if (global_fid == -1) { + return false; + } + m_edge_attribute[global_fid] = f_attr; + } + + return true; +} + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/EdgeSplitting.cpp b/components/triwild/wmtk/components/triwild/EdgeSplitting.cpp new file mode 100644 index 0000000000..e04a891dfc --- /dev/null +++ b/components/triwild/wmtk/components/triwild/EdgeSplitting.cpp @@ -0,0 +1,204 @@ +#include "TriWildMesh.h" + +#include +#include +#include +#include + +namespace wmtk::components::triwild { + +void TriWildMesh::split_all_edges() +{ + igl::Timer timer; + double time; + timer.start(); + auto collect_all_ops = std::vector>(); + for (const Tuple& loc : get_edges()) { + collect_all_ops.emplace_back("edge_split", loc); + } + time = timer.getElapsedTime(); + logger().info("edge split prepare time: {:.4}s", time); + auto setup_and_execute = [&](auto& executor) { + executor.renew_neighbor_tuples = + [](const TriWildMesh& m, std::string op, const auto& newts) { + std::vector> op_tups; + for (const auto& t : newts) { + op_tups.emplace_back(op, t); + op_tups.emplace_back(op, t.switch_edge(m)); + op_tups.emplace_back(op, t.switch_vertex(m).switch_edge(m)); + } + return op_tups; + }; + + executor.priority = [&](const TriWildMesh& m, std::string op, const Tuple& t) { + return m.get_length2(t); + }; + executor.num_threads = NUM_THREADS; + executor.is_weight_up_to_date = [&](const TriWildMesh& m, const auto& ele) { + auto [weight, op, tup] = ele; + auto length = m.get_length2(tup); + if (length != weight) { + return false; + } + // + size_t v1_id = tup.vid(*this); + size_t v2_id = tup.switch_vertex(*this).vid(*this); + const auto& VA = m_vertex_attribute; + double sizing_ratio = 0.5 * (VA[v1_id].m_sizing_scalar + VA[v2_id].m_sizing_scalar); + if (length < m_params.splitting_l2 * sizing_ratio * sizing_ratio) { + return false; + } + return true; + }; + executor(*this, collect_all_ops); + }; + if (NUM_THREADS > 0) { + timer.start(); + auto executor = ExecutePass(ExecutionPolicy::kPartition); + executor.lock_vertices = [&](auto& m, const auto& e, int task_id) -> bool { + return m.try_set_edge_mutex_two_ring(e, task_id); + }; + setup_and_execute(executor); + time = timer.getElapsedTime(); + wmtk::logger().info("edge split operation time parallel: {:.4}s", time); + } else { + timer.start(); + auto executor = ExecutePass(ExecutionPolicy::kSeq); + setup_and_execute(executor); + time = timer.getElapsedTime(); + wmtk::logger().info("edge split operation time serial: {:.4}s", time); + } +} + +bool TriWildMesh::split_edge_before(const Tuple& loc0) +{ + auto& cache = split_cache.local(); + + cache.changed_edges.clear(); + cache.faces.clear(); + + cache.v1_id = loc0.vid(*this); + cache.v2_id = loc0.switch_vertex(*this).vid(*this); + + cache.old_e_attrs = m_edge_attribute[loc0.eid(*this)]; + + const simplex::Edge edge(cache.v1_id, cache.v2_id); + + const auto faces = get_incident_fids_for_edge(loc0); + for (const size_t fid : faces) { + auto vs = oriented_tri_vids(fid); + for (int j = 0; j < 3; j++) { + const simplex::Edge e(vs[j], vs[(j + 1) % 3]); + if (e == edge) { + continue; + } + if (cache.changed_edges.count(e) != 0) { + continue; + } + auto [_, eid] = tuple_from_edge(e.vertices()); + cache.changed_edges[e] = m_edge_attribute[eid]; + } + } + + // store tet attributes + for (const size_t fid : faces) { + const simplex::Face face = simplex_from_face(fid); + const size_t opp = face.opposite_vertex(edge).id(); + cache.faces[opp] = m_face_attribute.at(fid); + } + + return true; +} + +bool TriWildMesh::split_edge_after(const Tuple& loc) +{ + if (!TriMesh::split_edge_after( + loc)) // note: call from super class, cannot be done with pure virtual classes + return false; + + const std::vector locs = get_one_ring_tris_for_vertex(loc.switch_vertex(*this)); + const size_t v_id = loc.switch_vertex(*this).vid(*this); + + auto& cache = split_cache.local(); + + const size_t v1_id = cache.v1_id; + const size_t v2_id = cache.v2_id; + + /// check inversion & rounding + auto& p = m_vertex_attribute[v_id].m_posf; + p = (m_vertex_attribute[v1_id].m_posf + m_vertex_attribute[v2_id].m_posf) / 2; + m_vertex_attribute[v_id].m_is_rounded = true; + + // this has to be done before the inversion check + m_vertex_attribute[v_id].m_pos = to_rational(p); + + for (auto& loc : locs) { + if (is_inverted(loc)) { + m_vertex_attribute[v_id].m_is_rounded = false; + break; + } + } + if (!m_vertex_attribute[v_id].m_is_rounded) { + m_vertex_attribute[v_id].m_pos = + (m_vertex_attribute[v1_id].m_pos + m_vertex_attribute[v2_id].m_pos) / 2; + p = to_double(m_vertex_attribute[v_id].m_pos); + } + + // update face attributes + { + // v1 - v_new + const auto faces1 = get_incident_fids_for_edge(v1_id, v_id); + const simplex::Edge edge1(v1_id, v_id); + for (const size_t fid : faces1) { + const simplex::Face face = simplex_from_face(fid); + const size_t opp = face.opposite_vertex(edge1).id(); + m_face_attribute[fid] = cache.faces[opp]; + } + // v2 - v_new + const auto faces2 = get_incident_fids_for_edge(v2_id, v_id); + const simplex::Edge edge2(v2_id, v_id); + for (const size_t fid : faces2) { + const simplex::Face face = simplex_from_face(fid); + const size_t opp = face.opposite_vertex(edge2).id(); + m_face_attribute[fid] = cache.faces[opp]; + } + assert(faces1.size() + faces2.size() == locs.size()); + + const auto [_1, eid1] = tuple_from_edge(edge1.vertices()); + const auto [_2, eid2] = tuple_from_edge(edge2.vertices()); + + m_edge_attribute[eid1] = cache.old_e_attrs; + m_edge_attribute[eid2] = cache.old_e_attrs; + for (const auto& [vid, _] : cache.faces) { + const auto [_tup, eid] = tuple_from_edge({{v_id, vid}}); + m_edge_attribute[eid].reset(); + } + } + + /// update quality + for (const Tuple& loc : locs) { + m_face_attribute[loc.fid(*this)].m_quality = get_quality(loc); + } + + /// update vertex attribute + // bbox + m_vertex_attribute[v_id].on_bbox_faces = wmtk::set_intersection( + m_vertex_attribute[v1_id].on_bbox_faces, + m_vertex_attribute[v2_id].on_bbox_faces); + // surface + m_vertex_attribute[v_id].m_is_on_surface = cache.old_e_attrs.m_is_surface_fs; + + /// update edge attribute + for (const auto& [e, e_attr] : cache.changed_edges) { + auto [_, eid] = tuple_from_edge(e.vertices()); + m_edge_attribute[eid] = e_attr; + } + + 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; + + return true; +} + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/EdgeSwapping.cpp b/components/triwild/wmtk/components/triwild/EdgeSwapping.cpp new file mode 100644 index 0000000000..8f14f82414 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/EdgeSwapping.cpp @@ -0,0 +1,178 @@ +#include "TriWildMesh.h" + +#include +#include +#include +#include +#include +#include "wmtk/utils/TupleUtils.hpp" + +#include + +namespace wmtk::components::triwild { + +auto renew = [](const TriWildMesh& m, auto op, auto& tris) { + using Tuple = TriMesh::Tuple; + std::vector edges; + for (const auto& t : tris) { + for (auto j = 0; j < 3; j++) { + edges.push_back(m.tuple_from_edge(t.fid(m), j)); + } + } + wmtk::unique_edge_tuples(m, edges); + + std::vector> optup; + optup.reserve(edges.size()); + for (const Tuple& e : edges) { + optup.emplace_back(op, e); + } + return optup; +}; + +auto edge_locker = [](auto& m, const auto& e, int task_id) { + // TODO: this should not be here + return m.try_set_edge_mutex_two_ring(e, task_id); +}; + +size_t TriWildMesh::swap_all_edges() +{ + std::vector> collect_all_ops; + { + igl::Timer timer; + timer.start(); + const auto edges = get_edges(); + collect_all_ops.reserve(edges.size()); + for (const Tuple& t : edges) { + collect_all_ops.emplace_back("edge_swap", t); + } + timer.stop(); + logger().info("edge collapse prepare time: {:.4}s", timer.getElapsedTimeInSec()); + } + logger().info("#E = {}", collect_all_ops.size()); + + auto setup_and_execute = [&](auto& executor) { + executor.renew_neighbor_tuples = renew; + executor.num_threads = NUM_THREADS; + executor.priority = [](const TriWildMesh& m, std::string op, const Tuple& e) { + return m.swap_weight(e); + }; + executor.should_renew = [](auto val) { return (val > 0); }; + executor.is_weight_up_to_date = [](const TriWildMesh& m, auto& ele) { + auto& [val, _, e] = ele; + const double w = m.swap_weight(e); + return (w > 1e-5) && ((w - val) * (w - val) < 1e-8); + }; + executor(*this, collect_all_ops); + }; + if (NUM_THREADS > 0) { + auto executor = ExecutePass(ExecutionPolicy::kPartition); + executor.lock_vertices = edge_locker; + setup_and_execute(executor); + } else { + auto executor = ExecutePass(ExecutionPolicy::kSeq); + setup_and_execute(executor); + } + + return true; +} + +double TriWildMesh::swap_weight(const Tuple& t) const +{ + const SmartTuple tt(*this, t); + const auto t_opp = tt.switch_face(); + if (!t_opp) { + return std::numeric_limits::lowest(); + } + + if (is_edge_on_surface(t)) { + return std::numeric_limits::lowest(); + } + + const size_t v0 = tt.vid(); + const size_t v1 = tt.switch_vertex().vid(); + const size_t v2 = tt.switch_edge().switch_vertex().vid(); + const size_t v3 = t_opp.value().switch_edge().switch_vertex().vid(); + + // before swap + const double q012 = get_quality({{v0, v1, v2}}); + const double q031 = get_quality({{v0, v3, v1}}); + // after swap + const double q032 = get_quality({{v0, v3, v2}}); + const double q231 = get_quality({{v2, v3, v1}}); + + const double q_before = std::max(q012, q031); + const double q_after = std::max(q032, q231); + + return q_before - q_after; +} + +bool TriWildMesh::swap_edge_before(const Tuple& t) +{ + if (is_edge_on_surface(t)) { + return false; + } + + const auto& FA = m_face_attribute; + auto& cache = swap_cache.local(); + cache.changed_edges.clear(); + + const auto incident_faces = get_incident_fids_for_edge(t); + + cache.face_tags = FA[incident_faces[0]].tags; + + double max_energy = -1.0; + for (const size_t fid : incident_faces) { + max_energy = std::max(FA[fid].m_quality, max_energy); + } + cache.max_energy = max_energy; + + // cache edges + simplex::Edge edge = simplex_from_edge(t); + for (const size_t fid : incident_faces) { + for (int j = 0; j < 3; j++) { + const Tuple tup = tuple_from_edge(fid, j); + simplex::Edge e = simplex_from_edge(tup); + if (e == edge) { + continue; + } + cache.changed_edges.try_emplace(e, m_edge_attribute[tup.eid(*this)]); + } + } + + return true; +} + +bool TriWildMesh::swap_edge_after(const Tuple& t) +{ + auto& cache = swap_cache.local(); + const auto incident_faces = get_incident_fids_for_edge(t); + + auto& FA = m_face_attribute; + + double max_energy = -1.0; + for (const size_t fid : incident_faces) { + if (is_inverted(fid)) { + return false; + } + double q = get_quality(fid); + FA[fid].m_quality = q; + max_energy = std::max(q, max_energy); + + FA[fid].tags = cache.face_tags; + } + if (max_energy >= cache.max_energy) { + return false; + } + + // cached edges + for (const auto& [e, e_attrs] : cache.changed_edges) { + const auto [_, eid] = tuple_from_edge(e.vertices()); + m_edge_attribute[eid] = e_attrs; + } + m_edge_attribute[t.eid(*this)].reset(); + + return true; +} + + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/Parameters.h b/components/triwild/wmtk/components/triwild/Parameters.h new file mode 100644 index 0000000000..7d91913d49 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/Parameters.h @@ -0,0 +1,74 @@ +#pragma once + +namespace wmtk::components::triwild { +struct Parameters +{ + double epsr = 2e-3; // relative error bound (wrt diagonal) + double eps = -1.; // absolute error bound + double lr = 5e-2; // target edge length (relative) + double l = -1.; + double l_min = -1; + double diag_l = -1.; + Vector2d box_min = Vector2d::Zero(); + Vector2d box_max = Vector2d::Ones(); + bool preserve_topology = false; + std::string output_path; + + double splitting_l2 = -1.; // the lower bound length (squared) for edge split + double collapsing_l2 = + std::numeric_limits::max(); // the upper bound length (squared) for edge collapse + + double stop_energy = 10; + + bool debug_output = false; + bool perform_sanity_checks = false; + + double w_amips = 1e-4; + double w_envelope = 0; + + void init(const Vector2d& min_, const Vector2d& max_) + { + box_min = min_; + box_max = max_; + diag_l = (box_max - box_min).norm(); + if (l > 0) { + lr = l / diag_l; + } else { + l = lr * diag_l; + } + splitting_l2 = l * l * (16 / 9.); + collapsing_l2 = l * l * (16 / 25.); + + if (eps > 0) { + epsr = eps / diag_l; + } else { + eps = epsr * diag_l; + } + + l_min = eps; + } + void init( + const std::vector& vertices, + const std::vector>& faces) + { + Vector2d min_, max_; + for (size_t i = 0; i < vertices.size(); i++) { + if (i == 0) { + min_ = vertices[i]; + max_ = vertices[i]; + continue; + } + for (int j = 0; j < 2; j++) { + if (vertices[i][j] < min_[j]) { + min_[j] = vertices[i][j]; + } + if (vertices[i][j] > max_[j]) { + max_[j] = vertices[i][j]; + } + } + } + + init(min_, max_); + } +}; +} // namespace wmtk::components::triwild diff --git a/components/triwild/wmtk/components/triwild/Smooth.cpp b/components/triwild/wmtk/components/triwild/Smooth.cpp new file mode 100644 index 0000000000..e8affeec85 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/Smooth.cpp @@ -0,0 +1,262 @@ + +#include "TriWildMesh.h" +#include "wmtk/ExecutionScheduler.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace wmtk::components::triwild { + +bool TriWildMesh::smooth_before(const Tuple& t) +{ + // try to round. + const bool r = round(t); + + const size_t vid = t.vid(*this); + + if (!m_vertex_attribute[vid].on_bbox_faces.empty()) { + return false; + } + + if (m_vertex_attribute[vid].m_is_rounded) { + return true; + } + // Note: no need to roll back. + return r; +} + + +bool TriWildMesh::smooth_after(const Tuple& t) +{ + // Newton iterations are encapsulated here. + logger().trace("Newton iteration for vertex smoothing."); + const size_t vid = t.vid(*this); + + const auto& VA = m_vertex_attribute; + + const auto& locs = get_one_ring_fids_for_vertex(t); + assert(locs.size() > 0); + + double max_quality = 0.; + for (const size_t fid : locs) { + max_quality = std::max(max_quality, m_face_attribute[fid].m_quality); + if (is_inverted_f(fid)) { + // Neighbors that are not rounded could cause a face to be inverted in floats + return false; + } + } + + auto& solver = m_solver.local(); + if (!solver) { + solver = optimization::create_basic_solver(); + } + + std::vector> assembles = get_amips_assembles(t); + + const Vector2d old_pos = VA[vid].m_posf; + + auto amips_energy = get_amips_energy(t); + std::shared_ptr total_energy = amips_energy; + + auto solve = [&]() { + VectorXd x = VA[vid].m_posf; + try { + solver->minimize(*total_energy, x); + } catch (const std::exception&) { + // polysolve might throw errors that we want to ignore (e.g., line search failed) + } + m_vertex_attribute[vid].m_posf = x; + m_vertex_attribute[vid].m_pos = to_rational(Vector2d(x)); + }; + + const auto surf_assembles = get_surface_assembles(t); + if (VA[vid].m_is_on_surface) { + assert(!surf_assembles.empty()); + + auto energy_sum = std::make_shared(); + + auto envelope_energy = get_envelope_energy(t); + // do one solve without weights for amips and envelope + if (m_params.w_amips > 0) { + energy_sum->add_energy(amips_energy, 1. / m_params.w_amips); + } + if (m_params.w_envelope > 0) { + energy_sum->add_energy(envelope_energy, 1. / m_params.w_envelope); + } + total_energy = energy_sum; + solve(); + + // second solve (with proper weights) + energy_sum = std::make_shared(); + + if (m_params.w_amips > 0) { + energy_sum->add_energy(amips_energy); + } + if (m_params.w_envelope > 0) { + energy_sum->add_energy(envelope_energy); + } + total_energy = energy_sum; + solve(); + } else { + // not on the surface + solve(); + } + + logger().trace("old pos {} -> new pos {}", old_pos.transpose(), VA[vid].m_posf.transpose()); + + // check surface containment + if (VA[vid].m_is_on_surface) { + for (size_t i = 1; i < surf_assembles.size(); ++i) { + std::array edge; + edge[0] = VA[vid].m_posf; + edge[1] = surf_assembles[i]; + if (m_envelope->is_outside(edge)) { + return false; + } + } + } + + // quality (only check if not on surface) + auto max_after_quality = 0.; + for (const size_t fid : locs) { + if (is_inverted(fid)) { + return false; + } + const double q = get_quality(fid); + m_face_attribute[fid].m_quality = q; + max_after_quality = std::max(max_after_quality, q); + } + if (!VA[vid].m_is_on_surface) { + if (max_after_quality > max_quality) { + return false; + } + } + + return true; +} + +void TriWildMesh::smooth_all_vertices(const size_t n_iters) +{ + for (size_t i = 0; i < n_iters; ++i) { + // log_total_surface_energy(); + igl::Timer timer; + timer.start(); + std::vector> collect_all_ops; + for (const Tuple& t : get_vertices()) { + collect_all_ops.emplace_back("vertex_smooth", t); + } + logger().info("vertex smoothing prepare time: {:.4}s", timer.getElapsedTimeInSec()); + logger().info("#V = {}", collect_all_ops.size()); + if (NUM_THREADS > 0) { + timer.start(); + ExecutePass executor(ExecutionPolicy::kPartition); + executor.lock_vertices = [](auto& m, const auto& e, int task_id) -> bool { + return m.try_set_vertex_mutex_one_ring(e, task_id); + }; + executor.num_threads = NUM_THREADS; + executor(*this, collect_all_ops); + logger().info("vertex smoothing time parallel: {:.4}s", timer.getElapsedTimeInSec()); + } else { + timer.start(); + ExecutePass executor(ExecutionPolicy::kSeq); + executor(*this, collect_all_ops); + logger().info("vertex smoothing time serial: {:.4}s", timer.getElapsedTimeInSec()); + } + if (m_params.debug_output) { + write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); + } + } +} + +std::vector TriWildMesh::get_surface_assembles(const Tuple& t) const +{ + const auto& VA = m_vertex_attribute; + const size_t vid = t.vid(*this); + + std::vector surface_pts; + + if (!VA[vid].m_is_on_surface) { + return surface_pts; + } + + const auto es = get_surface_edges_for_vertex(vid); + + surface_pts.resize(es.size() + 1); + surface_pts[0] = VA[vid].m_posf; + + for (size_t i = 0; i < es.edges().size(); ++i) { + const auto& vs = es.edges()[i].vertices(); + size_t neighbor_id = vs[0] != vid ? vs[0] : vs[1]; + surface_pts[i + 1] = VA[neighbor_id].m_posf; + } + + return surface_pts; +} + +std::shared_ptr TriWildMesh::get_envelope_energy( + const Tuple& t) const +{ + const double w = m_s_envelope * m_params.w_envelope; + + auto envelope_energy = std::make_shared(m_envelope, w); + return envelope_energy; +} + +std::vector> TriWildMesh::get_amips_assembles(const Tuple& t) const +{ + const size_t vid = t.vid(*this); + const auto& locs = get_one_ring_fids_for_vertex(t); + + const auto& VA = m_vertex_attribute; + + std::vector> assembles; + assembles.reserve(locs.size()); + + for (const size_t fid : locs) { + if (is_inverted(fid)) { + log_and_throw_error("Inverted face in amips assemble!"); + } + std::array local_verts = oriented_tri_vids(fid); + { + size_t v_loc = 0; + for (size_t i = 0; i < 3; ++i) { + if (local_verts[i] == vid) { + v_loc = i; + break; + } + } + std::array buf = local_verts; + local_verts[0] = buf[v_loc]; + local_verts[1] = buf[(v_loc + 1) % 3]; + local_verts[2] = buf[(v_loc + 2) % 3]; + } + + std::array T; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 2; j++) { + T[i * 2 + j] = VA[local_verts[i]].m_posf[j]; + } + } + assembles.push_back(T); + } + + return assembles; +} + +std::shared_ptr TriWildMesh::get_amips_energy(const Tuple& t) const +{ + const double w = m_params.w_amips > 0 ? m_s_amips * m_params.w_amips : 1; + + const auto assembles = get_amips_assembles(t); + auto amips_energy = std::make_shared(assembles, w); + assert(amips_energy->initial_position() == m_vertex_attribute.at(t.vid(*this)).m_posf); + return amips_energy; +} + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/TriWildMesh.cpp b/components/triwild/wmtk/components/triwild/TriWildMesh.cpp new file mode 100644 index 0000000000..bb5eb91eeb --- /dev/null +++ b/components/triwild/wmtk/components/triwild/TriWildMesh.cpp @@ -0,0 +1,1086 @@ + +#include "TriWildMesh.h" + +#include "wmtk/utils/Rational.hpp" + +#include +#include +#include +#include +#include +#include +#include + +// clang-format off +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// clang-format on + +//#include +#include +#include +#include + +namespace wmtk::components::triwild { + +VertexAttributes::VertexAttributes(const Vector2r& p) +{ + m_pos = p; + m_posf = to_double(p); +} + +void TriWildMesh::mesh_improvement(int max_its) +{ + ////preprocessing + partition_mesh_morton(); + + logger().info("========it pre========"); + local_operations({{0, 1, 0, 0}}, false); + + ////operation loops + bool is_hit_min_edge_length = false; + const int M = 2; + int m = 0; + double pre_max_energy = 0., pre_avg_energy = 0.; + for (int it = 0; it < max_its; it++) { + ///ops + logger().info("\n========it {}========", it); + auto [max_energy, avg_energy] = local_operations({{1, 1, 1, 1}}); + + ///energy check + logger().info("max energy {:.6} | stop {:.6}", max_energy, m_params.stop_energy); + if (max_energy < m_params.stop_energy) { + break; + } + consolidate_mesh(); + + logger().info("#V = {}, #T = {}", vert_capacity(), tri_capacity()); + + auto cnt_round = 0, cnt_verts = 0; + TriMesh::for_each_vertex([&](auto& v) { + if (m_vertex_attribute[v.vid(*this)].m_is_rounded) cnt_round++; + cnt_verts++; + }); + if (cnt_round < cnt_verts) { + logger().info("rounded {}/{}", cnt_round, cnt_verts); + } else { + logger().info("All rounded!", cnt_round, cnt_verts); + } + + ///sizing field + if (it > 0 && pre_max_energy - max_energy < 5e-1 && + (pre_avg_energy - avg_energy) / avg_energy < 0.1) { + m++; + if (m == M) { + logger().info(">>>>adjust_sizing_field..."); + is_hit_min_edge_length = adjust_sizing_field_serial(max_energy); + // is_hit_min_edge_length = adjust_sizing_field(max_energy); + logger().info(">>>>adjust_sizing_field finished..."); + m = 0; + } + } else { + m = 0; + pre_max_energy = max_energy; + pre_avg_energy = avg_energy; + } + if (is_hit_min_edge_length) { + // todo: maybe to do sth + } + } + + logger().info("========it post========"); + local_operations({{0, 1, 0, 0}}); +} + +std::tuple TriWildMesh::local_operations( + const std::array& ops, + bool collapse_limit_length) +{ + igl::Timer timer; + + std::tuple energy; + + auto sanity_checks = [this]() { + if (!m_params.perform_sanity_checks) { + return; + } + logger().info("Perform sanity checks..."); + const auto faces = get_edges_by_condition([](auto& f) { return f.m_is_surface_fs; }); + for (const auto& verts : faces) { + const auto& p0 = m_vertex_attribute[verts[0]].m_posf; + const auto& p1 = m_vertex_attribute[verts[1]].m_posf; + if (m_envelope->is_outside(std::array{p0, p1})) { + logger().error("Edge {} is outside!", verts); + } + } + + // check for inverted faces + for (const Tuple& t : get_faces()) { + if (!is_inverted(t)) { + continue; + } + const auto vs = oriented_tri_vids(t); + logger().error("Face {} is inverted! Vertices = {}", t.fid(*this), vs); + } + logger().info("Sanity checks done."); + }; + + sanity_checks(); + + timer.start(); + for (int i = 0; i < ops.size(); i++) { + if (i == 0) { + for (int n = 0; n < ops[i]; n++) { + logger().info("==splitting {}==", n); + split_all_edges(); + logger().info( + "#V = {}, #F = {} after split", + get_vertices().size(), + get_faces().size()); + } + if (m_params.debug_output) { + write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); + } + auto [max_energy, avg_energy] = get_max_avg_energy(); + logger().info("split max energy = {:.6} avg = {:.6}", max_energy, avg_energy); + sanity_checks(); + } else if (i == 1) { + for (int n = 0; n < ops[i]; n++) { + logger().info("==collapsing {}==", n); + collapse_all_edges(collapse_limit_length); + logger().info( + "#V = {}, #F = {} after collapse", + get_vertices().size(), + get_faces().size()); + } + if (m_params.debug_output) { + write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); + } + auto [max_energy, avg_energy] = get_max_avg_energy(); + logger().info("collapse max energy = {:.6} avg = {:.6}", max_energy, avg_energy); + sanity_checks(); + } else if (i == 2) { + for (int n = 0; n < ops[i]; n++) { + logger().info("==swapping {}==", n); + size_t cnt_success = swap_all_edges(); + if (cnt_success == 0) { + break; + } + } + if (m_params.debug_output) { + write_vtu(fmt::format("debug_{}", m_debug_print_counter++)); + } + auto [max_energy, avg_energy] = get_max_avg_energy(); + logger().info("swap max energy = {:.6} avg = {:.6}", max_energy, avg_energy); + sanity_checks(); + } else if (i == 3) { + logger().info("==smoothing =="); + smooth_all_vertices(ops[i]); + auto [max_energy, avg_energy] = get_max_avg_energy(); + logger().info("smooth max energy = {:.6} avg = {:.6}", max_energy, avg_energy); + sanity_checks(); + } + } + energy = get_max_avg_energy(); + logger().info("max energy = {:.6}", std::get<0>(energy)); + logger().info("avg energy = {:.6}", std::get<1>(energy)); + logger().info("time = {:.4}s", timer.getElapsedTimeInSec()); + + + return energy; +} + +void TriWildMesh::init_mesh(const MatrixXd& V, const MatrixXi& F, const MatrixXi& E) +{ + assert(V.cols() == 2); + assert(F.cols() == 3); + assert(E.cols() == 2); + + init(F); + + assert(check_mesh_connectivity_validity()); + + m_vertex_attribute.m_attributes.resize(V.rows()); + m_edge_attribute.m_attributes.resize(F.rows() * 3); + + for (int i = 0; i < vert_capacity(); i++) { + m_vertex_attribute[i].m_pos = to_rational(Vector2d(V.row(i))); + m_vertex_attribute[i].m_posf = V.row(i); + } + + // init quality and check for inverted mesh + bool is_mesh_inverted = false; + for (const Tuple& t : get_faces()) { + if (is_mesh_inverted ^ is_inverted(t)) { + if (!is_mesh_inverted) { + is_mesh_inverted = true; + } else { + log_and_throw_error("Tets with different orientations in the input!"); + } + } + m_face_attribute[t.fid(*this)].m_quality = get_quality(t); + } + + if (is_mesh_inverted) { + log_and_throw_error( + "Input mesh is fully inverted! This should not happen... Might be a bug."); + } + + // mark edges as on surface if they are in E + for (int i = 0; i < E.rows(); i++) { + std::array vids = {(size_t)E(i, 0), (size_t)E(i, 1)}; + const auto [e, eid] = tuple_from_edge(vids); + if (!e.is_valid(*this)) { + log_and_throw_error("Edge {} in E is not found in the mesh!", vids); + } + m_edge_attribute[eid].m_is_surface_fs = true; + m_vertex_attribute[vids[0]].m_is_on_surface = true; + m_vertex_attribute[vids[1]].m_is_on_surface = true; + } + + // init envelope + if (m_envelope) { + log_and_throw_error("Envelope was already initialized once."); + } + assert(m_V_envelope.empty() && m_E_envelope.empty()); + + m_V_envelope.resize(V.rows()); + for (size_t i = 0; i < m_V_envelope.size(); ++i) { + m_V_envelope[i] = V.row(i); + } + m_E_envelope.resize(E.rows()); + for (size_t i = 0; i < m_E_envelope.size(); ++i) { + m_E_envelope[i] = E.row(i); + } + + m_envelope = std::make_shared(); + m_envelope->init(m_V_envelope, m_E_envelope, m_envelope_eps); + + // Sanity check: All surface edges must be inside the envelope + { + logger().info("Envelope sanity check"); + const auto surf_edges = get_edges_by_condition([](auto& f) { return f.m_is_surface_fs; }); + for (const auto& verts : surf_edges) { + std::array pp = { + m_vertex_attribute[verts[0]].m_posf, + m_vertex_attribute[verts[1]].m_posf}; + if (m_envelope->is_outside(pp)) { + log_and_throw_error("Edge {} is outside!", verts); + } + } + logger().info("Envelope sanity check done"); + } + + // track bounding box + const auto edges = get_edges(); + for (size_t i = 0; i < edges.size(); i++) { + const auto vids = get_edge_vids(edges[i]); + int on_bbox = -1; + for (int k = 0; k < 2; k++) { + if (m_vertex_attribute[vids[0]].m_pos[k] == m_params.box_min[k] && + m_vertex_attribute[vids[1]].m_pos[k] == m_params.box_min[k]) { + on_bbox = k * 2; + break; + } + if (m_vertex_attribute[vids[0]].m_pos[k] == m_params.box_max[k] && + m_vertex_attribute[vids[1]].m_pos[k] == m_params.box_max[k]) { + on_bbox = k * 2 + 1; + break; + } + } + if (on_bbox < 0) { + continue; + } + if (edges[i].switch_face(*this)) { + log_and_throw_error("Boundary edge {} is not on the boundary!", vids); + } + + const size_t eid = edges[i].eid(*this); + m_edge_attribute[eid].m_is_bbox_fs = on_bbox; + + for (const size_t vid : vids) { + m_vertex_attribute[vid].on_bbox_faces.push_back(on_bbox); + } + } + + for_each_vertex( + [&](auto& v) { wmtk::vector_unique(m_vertex_attribute[v.vid(*this)].on_bbox_faces); }); + + //// rounding + size_t cnt_round = 0; + + for (int i = 0; i < vert_capacity(); i++) { + Tuple v = tuple_from_vertex(i); + if (round(v)) { + cnt_round++; + } + } + + if (cnt_round < vert_capacity()) { + logger().info("Rounded {}/{}", cnt_round, vert_capacity()); + } else { + logger().info("All rounded!", cnt_round, vert_capacity()); + } +} + +bool TriWildMesh::adjust_sizing_field_serial(double max_energy) +{ + logger().info("#V {}, #F {}", vert_capacity(), tri_capacity()); + + const double stop_filter_energy = m_params.stop_energy * 0.8; + double filter_energy = std::max(max_energy / 100, stop_filter_energy); + filter_energy = std::min(filter_energy, 100.); + + const auto recover_scalar = 1.5; + const auto refine_scalar = 0.5; + const auto min_refine_scalar = m_params.l_min / m_params.l; + + // outputs scale_multipliers + std::vector scale_multipliers(vert_capacity(), recover_scalar); + + std::vector pts; + std::queue v_queue; + + for (int i = 0; i < tri_capacity(); i++) { + const Tuple t = tuple_from_tri(i); + if (!t.is_valid(*this)) { + continue; + } + const size_t fid = t.fid(*this); + if (m_face_attribute.at(fid).m_quality < filter_energy) { + continue; + } + const auto vs = oriented_tri_vids(t); + Vector2d c(0, 0); // center + for (int j = 0; j < 3; j++) { + c += m_vertex_attribute.at(vs[j]).m_posf; + v_queue.emplace(vs[j]); + } + c /= 3; + pts.emplace_back(Vector3d(c[0], c[1], 0)); + } + + logger().info("filter energy {} Low Quality Tets {}", filter_energy, pts.size()); + + const double R = m_params.l * 1.8; + + int sum = 0; + int adjcnt = 0; + + std::vector visited(vert_capacity(), false); + + KNN knn(pts); + + std::vector cache_one_ring; + // size_t vid; + while (!v_queue.empty()) { + sum++; + const size_t vid = v_queue.front(); + v_queue.pop(); + if (visited[vid]) continue; + visited[vid] = true; + adjcnt++; + + const Vector2d& pos_v = m_vertex_attribute.at(vid).m_posf; + const Vector3d p(pos_v[0], pos_v[1], 0); + double sq_dist = 0.; + uint32_t idx; + knn.nearest_neighbor(p, idx, sq_dist); + const double dist = std::sqrt(sq_dist); + + if (dist > R) { // outside R-ball, unmark. + continue; + } + + scale_multipliers[vid] = std::min( + scale_multipliers[vid], + dist / R * (1 - refine_scalar) + refine_scalar); // linear interpolate + + get_one_ring_vids_for_vertex_duplicate(vid, cache_one_ring); + for (size_t n_vid : cache_one_ring) { + if (visited[n_vid]) { + continue; + } + v_queue.push(n_vid); + } + } + + logger().info("sum = {}; adjacent = {}", sum, adjcnt); + + std::atomic_bool is_hit_min_edge_length = false; + + for (int i = 0; i < vert_capacity(); i++) { + const Tuple v = tuple_from_vertex(i); + if (!v.is_valid(*this)) { + continue; + } + const size_t vid = v.vid(*this); + auto& v_attr = m_vertex_attribute[vid]; + + auto new_scale = v_attr.m_sizing_scalar * scale_multipliers[vid]; + if (new_scale > 1) { + v_attr.m_sizing_scalar = 1; + } else if (new_scale < min_refine_scalar) { + is_hit_min_edge_length = true; + v_attr.m_sizing_scalar = min_refine_scalar; + } else { + v_attr.m_sizing_scalar = new_scale; + } + } + + return is_hit_min_edge_length.load(); +} + +void TriWildMesh::write_msh_groups(std::string file, const bool write_envelope) +{ + consolidate_mesh(); + + wmtk::MshData msh; + + const auto& vtx = get_vertices(); + msh.add_face_vertices(vtx.size(), [&](size_t k) { + auto i = vtx[k].vid(*this); + Vector2d p2 = m_vertex_attribute[i].m_posf; + return Vector3d(p2[0], p2[1], 0); + }); + + const auto& faces = get_faces(); + + int64_t max_tag = -1; + for (const Tuple& t : faces) { + const size_t fid = t.fid(*this); + const auto& tags = m_face_attribute[fid].tags; + if (tags.size() == 0) { + continue; + } + int64_t mt = *tags.rbegin(); + max_tag = std::max(max_tag, mt); + } + + if (m_tags_count < max_tag + 1) { + logger().warn( + "Max tag is {} but m_tags_count is {}. Adjusting m_tags_count.", + max_tag, + m_tags_count); + m_tags_count = max_tag + 1; + } + + std::vector faces_with_tag; + faces_with_tag.reserve(faces.size()); + + auto msh_add_faces = [&]() { + msh.add_faces(faces_with_tag.size(), [&](size_t k) { + auto vs = oriented_tri_vids(faces_with_tag[k]); + std::array data; + for (int j = 0; j < 3; j++) { + data[j] = vs[j]; + } + return data; + }); + }; + + // ambient mesh (no non-zero tags) + for (const Tuple& t : faces) { + const size_t fid = t.fid(*this); + if (m_face_attribute[fid].tags.empty()) { + faces_with_tag.push_back(t); + } + } + msh_add_faces(); + + msh.add_physical_group("ambient"); + + + // add a group for each tag + for (size_t tag_img = 0; tag_img < m_tags_count; ++tag_img) { + faces_with_tag.clear(); + for (const Tuple& t : faces) { + const size_t fid = t.fid(*this); + if (m_face_attribute[fid].tags.count(tag_img)) { + faces_with_tag.push_back(t); + } + } + + if (faces_with_tag.empty()) { + continue; + } + + msh.add_empty_vertices(2); + msh_add_faces(); + + std::string group_name; + if (m_tag_id_to_name.count(tag_img)) { + group_name = m_tag_id_to_name[tag_img]; + } else { + group_name = fmt::format("tag_{}", tag_img); + while (m_tag_name_to_id.count(group_name)) { + group_name += "_"; + } + m_tag_name_to_id[group_name] = tag_img; + m_tag_id_to_name[tag_img] = group_name; + logger().warn( + "Tag {} does not have a name. Assigning the name {}.", + tag_img, + group_name); + } + msh.add_physical_group(group_name); + } + + if (m_envelope && write_envelope) { + msh.add_edge_vertices(m_V_envelope.size(), [this](size_t k) { + return Vector3d(m_V_envelope[k][0], m_V_envelope[k][1], 0); + }); + msh.add_edges(m_E_envelope.size(), [this](size_t k) { return m_E_envelope[k]; }); + msh.add_physical_group("EnvelopeSurface"); + } + + logger().info("Write {}", file); + msh.save(file, true); +} + +void TriWildMesh::compute_winding_numbers(const std::vector& input_paths) +{ + const auto& faces = get_faces(); + MatrixXd C = MatrixXd::Zero(faces.size(), 2); + for (size_t i = 0; i < faces.size(); i++) { + const auto vs = oriented_tri_vids(faces[i]); + for (size_t v : vs) { + C.row(i) += m_vertex_attribute[v].m_posf; + } + C.row(i) /= 3; + } + + m_tags_count = input_paths.size(); + int64_t input_idx = 0; + for (const std::string& input_path : input_paths) { + MatrixXd V; + MatrixXi E; + io::read_edge_mesh(input_path, V, E); + assert(V.cols() == 3); + assert(E.cols() == 2); + + V = V.block(0, 0, V.rows(), 2).eval(); // only use x,y for winding number + + // compute winding number for V,F + Eigen::VectorXd W; + igl::winding_number(V, E, C, W); + + if (W.maxCoeff() <= 0.5) { + // all removed, let's invert. + logger().info("Correcting winding number"); + for (auto i = 0; i < E.rows(); i++) { + auto temp = E(i, 0); + E(i, 0) = E(i, 1); + E(i, 1) = temp; + } + igl::winding_number(V, E, C, W); + } + + if (W.maxCoeff() <= 0.5) { + logger().warn("No winding number above 0.5 for input_path {}", input_path); + } + + // store winding number in mesh + for (int i = 0; i < faces.size(); ++i) { + const size_t fid = faces[i].fid(*this); + if (W(i) > 0.5) { + m_face_attribute[fid].tags.insert(input_idx); + } + } + ++input_idx; + } +} + +int TriWildMesh::flood_fill() +{ + int current_id = 0; + const auto faces = get_faces(); + std::map visited; + + for (const Tuple& t : faces) { + size_t fid = t.fid(*this); + if (visited.find(fid) != visited.end()) { + continue; + } + + visited[fid] = true; + m_face_attribute[fid].part_id = current_id; + + const Tuple f1 = t; + const Tuple f2 = t.switch_edge(*this); + const Tuple f3 = t.switch_vertex(*this).switch_edge(*this); + + std::queue bfs_queue; + + if (!m_edge_attribute[f1.eid(*this)].m_is_surface_fs) { + auto oppo_t = f1.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + if (!m_edge_attribute[f2.eid(*this)].m_is_surface_fs) { + auto oppo_t = f2.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + if (!m_edge_attribute[f3.eid(*this)].m_is_surface_fs) { + auto oppo_t = f3.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + + while (!bfs_queue.empty()) { + auto tmp = bfs_queue.front(); + bfs_queue.pop(); + size_t tmp_id = tmp.fid(*this); + if (visited.find(tmp_id) != visited.end()) continue; + + visited[tmp_id] = true; + + m_face_attribute[tmp_id].part_id = current_id; + + const Tuple f_tmp1 = tmp; + const Tuple f_tmp2 = tmp.switch_edge(*this); + const Tuple f_tmp3 = tmp.switch_vertex(*this).switch_edge(*this); + + if (!m_edge_attribute[f_tmp1.eid(*this)].m_is_surface_fs) { + auto oppo_t = f_tmp1.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + if (!m_edge_attribute[f_tmp2.eid(*this)].m_is_surface_fs) { + auto oppo_t = f_tmp2.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + if (!m_edge_attribute[f_tmp3.eid(*this)].m_is_surface_fs) { + auto oppo_t = f_tmp3.switch_face(*this); + if (oppo_t.has_value()) { + if (visited.find((*oppo_t).fid(*this)) == visited.end()) { + bfs_queue.push(*oppo_t); + } + } + } + } + + current_id++; + } + return current_id; +} + +void TriWildMesh::partition_mesh() +{ + auto m_vertex_partition_id = partition_TriMesh(*this, NUM_THREADS); + for (size_t i = 0; i < m_vertex_partition_id.size(); i++) { + m_vertex_attribute[i].partition_id = m_vertex_partition_id[i]; + } +} + +void TriWildMesh::partition_mesh_morton() +{ + if (NUM_THREADS == 0) return; + logger().info("Number of parts: {} by morton", NUM_THREADS); + + tbb::task_arena arena(NUM_THREADS); + + arena.execute([&] { + std::vector V_v(vert_capacity()); + + tbb::parallel_for( + tbb::blocked_range(0, V_v.size()), + [&](tbb::blocked_range r) { + for (size_t i = r.begin(); i < r.end(); i++) { + V_v[i] = m_vertex_attribute[i].m_posf; + } + }); + + struct sortstruct + { + size_t order; + Resorting::MortonCode64 morton; + }; + + std::vector list_v; + list_v.resize(V_v.size()); + const int multi = 1000; + // since the morton code requires a correct scale of input vertices, + // we need to scale the vertices if their coordinates are out of range + std::vector V = V_v; // this is for rescaling vertices + Vector2d vmin, vmax; + vmin = V.front(); + vmax = V.front(); + + for (size_t j = 0; j < V.size(); j++) { + for (int i = 0; i < 2; i++) { + vmin(i) = std::min(vmin(i), V[j](i)); + vmax(i) = std::max(vmax(i), V[j](i)); + } + } + + Vector2d center = (vmin + vmax) / 2; + + tbb::parallel_for( + tbb::blocked_range(0, V.size()), + [&](tbb::blocked_range r) { + for (size_t i = r.begin(); i < r.end(); i++) { + V[i] = V[i] - center; + } + }); + + Vector2d scale_point = + vmax - center; // after placing box at origin, vmax and vmin are symetric. + + double xscale, yscale; + xscale = fabs(scale_point[0]); + yscale = fabs(scale_point[1]); + double scale = std::max(xscale, yscale); + if (scale > 300) { + tbb::parallel_for( + tbb::blocked_range(0, V.size()), + [&](tbb::blocked_range r) { + for (size_t i = r.begin(); i < r.end(); i++) { + V[i] = V[i] / scale; + } + }); + } + + tbb::parallel_for( + tbb::blocked_range(0, V.size()), + [&](tbb::blocked_range r) { + for (size_t i = r.begin(); i < r.end(); i++) { + list_v[i].morton = + Resorting::MortonCode64(int(V[i][0] * multi), int(V[i][1] * multi), 0); + list_v[i].order = i; + } + }); + + const auto morton_compare = [](const sortstruct& a, const sortstruct& b) { + return (a.morton < b.morton); + }; + + tbb::parallel_sort(list_v.begin(), list_v.end(), morton_compare); + + size_t interval = list_v.size() / NUM_THREADS + 1; + + tbb::parallel_for( + tbb::blocked_range(0, list_v.size()), + [&](tbb::blocked_range r) { + for (size_t i = r.begin(); i < r.end(); i++) { + m_vertex_attribute[list_v[i].order].partition_id = i / interval; + } + }); + }); +} + +double TriWildMesh::get_length2(const Tuple& l) const +{ + auto& m = *this; + auto& v1 = l; + auto v2 = l.switch_vertex(m); + double length = + (m.m_vertex_attribute[v1.vid(m)].m_posf - m.m_vertex_attribute[v2.vid(m)].m_posf) + .squaredNorm(); + return length; +} + +std::tuple TriWildMesh::get_max_avg_energy() +{ + double max_energy = -1.; + double avg_energy = 0.; + auto cnt = 0; + + for (int i = 0; i < tri_capacity(); i++) { + const Tuple tup = tuple_from_tri(i); + if (!tup.is_valid(*this)) { + continue; + } + const double q = m_face_attribute[tup.fid(*this)].m_quality; + max_energy = std::max(max_energy, q); + avg_energy += q; + cnt++; + } + + avg_energy /= cnt; + + return std::make_tuple(max_energy, avg_energy); +} + + +bool TriWildMesh::is_inverted_f(const Tuple& loc) const +{ + return is_inverted_f(loc.fid(*this)); +} + +bool TriWildMesh::is_inverted_f(const size_t fid) const +{ + auto vs = oriented_tri_vids(fid); + + igl::predicates::exactinit(); + auto res = igl::predicates::orient2d( + m_vertex_attribute[vs[0]].m_posf, + m_vertex_attribute[vs[1]].m_posf, + m_vertex_attribute[vs[2]].m_posf); + if (res == igl::predicates::Orientation::POSITIVE) { + return false; + } + return true; +} + +bool TriWildMesh::is_inverted(const std::array& vs) const +{ + // Return a positive value if the point pd lies below the + // plane passing through pa, pb, and pc; "below" is defined so + // that pa, pb, and pc appear in counterclockwise order when + // viewed from above the plane. + + if (m_vertex_attribute[vs[0]].m_is_rounded && m_vertex_attribute[vs[1]].m_is_rounded && + m_vertex_attribute[vs[2]].m_is_rounded) { + igl::predicates::exactinit(); + auto res = igl::predicates::orient2d( + m_vertex_attribute[vs[0]].m_posf, + m_vertex_attribute[vs[1]].m_posf, + m_vertex_attribute[vs[2]].m_posf); + if (res == igl::predicates::Orientation::POSITIVE) { + return false; + } + return true; + } else { + const Vector2r& v0 = m_vertex_attribute[vs[0]].m_pos; + const Vector2r& v1 = m_vertex_attribute[vs[1]].m_pos; + const Vector2r& v2 = m_vertex_attribute[vs[2]].m_pos; + const Vector2r a = v1 - v0; + const Vector2r b = v2 - v0; + Rational res = a.x() * b.y() - a.y() * b.x(); + if (res > 0) { + return false; + } else { + return true; + } + } +} + +bool TriWildMesh::is_inverted(const Tuple& loc) const +{ + auto vs = oriented_tri_vids(loc); + return is_inverted(vs); +} + +bool TriWildMesh::is_inverted(const size_t fid) const +{ + auto vs = oriented_tri_vids(fid); + return is_inverted(vs); +} + +bool TriWildMesh::round(const Tuple& v) +{ + size_t i = v.vid(*this); + if (m_vertex_attribute[i].m_is_rounded) { + return true; + } + + auto old_pos = m_vertex_attribute[i].m_pos; + m_vertex_attribute[i].m_pos << m_vertex_attribute[i].m_posf[0], m_vertex_attribute[i].m_posf[1]; + auto conn_tets = get_one_ring_tris_for_vertex(v); + m_vertex_attribute[i].m_is_rounded = true; + for (const Tuple& tet : conn_tets) { + if (is_inverted(tet)) { + m_vertex_attribute[i].m_is_rounded = false; + m_vertex_attribute[i].m_pos = old_pos; + return false; + } + } + + return true; +} + +double TriWildMesh::get_quality(const std::array& vs) const +{ + std::array ps; + for (size_t k = 0; k < 3; k++) { + ps[k] = m_vertex_attribute[vs[k]].m_posf; + } + double energy = -1.; + { + std::array T; + for (size_t k = 0; k < 3; k++) + for (size_t j = 0; j < 2; j++) { + T[k * 2 + j] = ps[k][j]; + } + energy = AMIPS2D_energy(T); + } + if (std::isinf(energy) || std::isnan(energy) || energy < 2 - 1e-3) { + return MAX_ENERGY; + } + return energy; +} + +double TriWildMesh::get_quality(const Tuple& loc) const +{ + auto its = oriented_tri_vids(loc); + return get_quality(its); +} + +double TriWildMesh::get_quality(const size_t fid) const +{ + auto its = oriented_tri_vids(fid); + return get_quality(its); +} + +void TriWildMesh::write_vtu(const std::string& path) const +{ + // consolidate_mesh(); + const std::string out_path = path + ".vtu"; + logger().info("Write {}", out_path); + const auto& faces = get_faces(); + const auto edges = get_edges_by_condition([](auto& f) { return f.m_is_surface_fs; }); + + Eigen::MatrixXd V(vert_capacity(), 2); + Eigen::MatrixXi F(tri_capacity(), 3); + Eigen::MatrixXi E(edges.size(), 2); + + V.setZero(); + F.setZero(); + E.setZero(); + + Eigen::VectorXd v_sizing_field(vert_capacity()); + v_sizing_field.setZero(); + + std::vector tags(m_tags_count, MatrixXd(tri_capacity(), 1)); + for (size_t j = 0; j < m_tags_count; ++j) { + tags[j].setZero(); + } + MatrixXd amips(tri_capacity(), 1); + amips.setZero(); + MatrixXd flood_fill(tri_capacity(), 1); + flood_fill.setZero(); + + int index = 0; + for (const Tuple& t : faces) { + size_t fid = t.fid(*this); + amips(index, 0) = m_face_attribute[fid].m_quality; + for (size_t j = 0; j < m_tags_count; ++j) { + tags[j](index, 0) = m_face_attribute[fid].tags.count(j) ? 1 : 0; + } + flood_fill(index, 0) = m_face_attribute[fid].part_id; + + const auto vs = oriented_tri_vids(t); + for (size_t j = 0; j < 3; j++) { + F(index, j) = (int)vs[j]; + } + ++index; + } + + for (size_t i = 0; i < edges.size(); ++i) { + for (size_t j = 0; j < 2; ++j) { + E(i, j) = (int)edges[i][j]; + } + } + + for (const Tuple& v : get_vertices()) { + const size_t vid = v.vid(*this); + V.row(vid) = m_vertex_attribute[vid].m_posf; + v_sizing_field[vid] = m_vertex_attribute[vid].m_sizing_scalar; + } + + paraviewo::VTUWriter writer; + + for (size_t j = 0; j < m_tags_count; ++j) { + writer.add_cell_field(fmt::format("tag_{}", j), tags[j]); + } + writer.add_cell_field("quality", amips); + writer.add_cell_field("flood_fill", flood_fill); + writer.add_field("sizing_field", v_sizing_field); + writer.write_mesh(path + ".vtu", V, F, paraviewo::CellType::Triangle); + + // surface + { + const auto surf_out_path = path + "_surf.vtu"; + paraviewo::VTUWriter surf_writer; + surf_writer.add_field("sizing_field", v_sizing_field); + + logger().info("Write {}", surf_out_path); + surf_writer.write_mesh(surf_out_path, V, E, paraviewo::CellType::Line); + } +} + +std::vector> TriWildMesh::get_edges_by_condition( + std::function cond) const +{ + std::vector> res; + for (const Tuple& e : get_edges()) { + size_t eid = e.eid(*this); + if (cond(m_edge_attribute[eid])) { + res.push_back({{e.vid(*this), e.switch_vertex(*this).vid(*this)}}); + } + } + return res; +} + +bool TriWildMesh::is_edge_on_surface(const Tuple& loc) const +{ + const auto vs = get_edge_vids(loc); + if (!m_vertex_attribute.at(vs[0]).m_is_on_surface || + !m_vertex_attribute.at(vs[1]).m_is_on_surface) { + return false; + } + + const size_t eid = loc.eid(*this); + return m_edge_attribute[eid].m_is_surface_fs; +} +bool TriWildMesh::is_edge_on_surface(const std::array& vids) const +{ + if (!m_vertex_attribute.at(vids[0]).m_is_on_surface || + !m_vertex_attribute.at(vids[1]).m_is_on_surface) { + return false; + } + + const auto [_, eid] = tuple_from_edge(vids); + return m_edge_attribute[eid].m_is_surface_fs; +} +bool TriWildMesh::is_edge_on_bbox(const Tuple& loc) const +{ + const auto vs = get_edge_vids(loc); + if (m_vertex_attribute.at(vs[0]).on_bbox_faces.empty() || + m_vertex_attribute.at(vs[1]).on_bbox_faces.empty()) { + return false; + } + + const size_t eid = loc.eid(*this); + return m_edge_attribute[eid].m_is_bbox_fs >= 0; +} + +bool TriWildMesh::is_edge_on_bbox(const std::array& vids) const +{ + if (m_vertex_attribute.at(vids[0]).on_bbox_faces.empty() || + m_vertex_attribute.at(vids[1]).on_bbox_faces.empty()) { + return false; + } + const auto [_, eid] = tuple_from_edge(vids); + return m_edge_attribute[eid].m_is_bbox_fs >= 0; +} + +} // namespace wmtk::components::triwild diff --git a/components/triwild/wmtk/components/triwild/TriWildMesh.h b/components/triwild/wmtk/components/triwild/TriWildMesh.h new file mode 100644 index 0000000000..3ef3a7a767 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/TriWildMesh.h @@ -0,0 +1,306 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +// clang-format off +#include +#include +#include +#include +#include +#include +#include +#include +// clang-format on + +#include "Parameters.h" + +namespace wmtk::components::triwild { + +// TODO: missing comments on what these attributes are +class VertexAttributes +{ +public: + Vector2d m_posf; + Vector2r m_pos; + bool m_is_rounded = false; + + bool m_is_on_surface = false; + std::vector on_bbox_faces; + + double m_sizing_scalar = 1; + + size_t partition_id = 0; + + VertexAttributes() {} + VertexAttributes(const Vector2r& p); +}; + +class EdgeAttributes +{ +public: + double tag; // TODO: is this used? + + bool m_is_surface_fs = false; // 0; 1 + /** + * Keep track which bbox side the face is on + * -1: none + * 0/1: x min/max + * 2/3: y min/max + * + * This bbox side ID is used to keep the bbox from collapsing. + */ + int m_is_bbox_fs = -1; //-1; 0~3 + + void reset() + { + m_is_surface_fs = false; + m_is_bbox_fs = -1; + } + + void merge(const EdgeAttributes& attr) + { + m_is_surface_fs = m_is_surface_fs || attr.m_is_surface_fs; + if (attr.m_is_bbox_fs >= 0) m_is_bbox_fs = attr.m_is_bbox_fs; + } +}; + +class FaceAttributes +{ +public: + double m_quality; + double m_winding_number = 0; + std::set tags; + int part_id = -1; +}; + +class TriWildMesh : public wmtk::TriMesh +{ +public: + int m_debug_print_counter = 0; + size_t m_tags_count = 0; + std::map m_tag_id_to_name; + std::map m_tag_name_to_id; + + const double MAX_ENERGY = 1e50; + + Parameters& m_params; + std::vector m_V_envelope; + std::vector m_E_envelope; + std::shared_ptr m_envelope; + double m_envelope_eps = -1; + + using VertAttCol = wmtk::AttributeCollection; + using EdgeAttCol = wmtk::AttributeCollection; + using FaceAttCol = wmtk::AttributeCollection; + VertAttCol m_vertex_attribute; + EdgeAttCol m_edge_attribute; + FaceAttCol m_face_attribute; + + tbb::enumerable_thread_specific> m_solver; + + // scaling factors + double m_s_amips = -1; + double m_s_envelope = -1; + + TriWildMesh(Parameters& _m_params, double envelope_eps, int _num_threads = 0) + : m_params(_m_params) + , m_envelope_eps(envelope_eps) + { + NUM_THREADS = _num_threads; + p_vertex_attrs = &m_vertex_attribute; + p_edge_attrs = &m_edge_attribute; + p_face_attrs = &m_face_attribute; + + optimization::deactivate_opt_logger(); + + m_s_amips = 1.; + /** + * eps makes it such that the energy is relative to the envelope thickness. As it's a + * squared energy, we need eps^2. + */ + m_s_envelope = 1. / (m_params.eps * m_params.eps); + + + double& wa = m_params.w_amips; + double& we = m_params.w_envelope; + we = 1 - wa; + logger().info("w_envelope = {}", we); + } + + ~TriWildMesh() {} + + // TODO: this should not be here + void partition_mesh(); + + // TODO: morton should not be here, but inside wmtk + void partition_mesh_morton(); + + size_t get_partition_id(const Tuple& loc) const + { + return m_vertex_attribute[loc.vid(*this)].partition_id; + } + + double get_length2(const Tuple& l) const; + +public: + /** + * @brief Init mesh from IGL-style matrices. + * + * @param V #Vx3 vertices of the tet mesh + * @param F #Fx3 vertex IDs for all faces + * @param E #Ex2 vertex IDs for all constraint edges + */ + void init_mesh(const MatrixXd& V, const MatrixXi& F, const MatrixXi& E); + + void init_surfaces_and_boundaries(); + + void init_envelope(const MatrixXd& V, const MatrixXi& F); + + bool adjust_sizing_field_serial(double max_energy); + + void write_msh_groups(std::string file, const bool write_envelope = true); + + void write_vtu(const std::string& path) const; + + std::vector> get_edges_by_condition( + std::function cond) const; + +public: + void split_all_edges(); + bool split_edge_before(const Tuple& t) override; + bool split_edge_after(const Tuple& loc) override; + + void collapse_all_edges(bool is_limit_length = true); + bool collapse_edge_before(const Tuple& t) override; + bool collapse_edge_after(const Tuple& t) override; + + size_t swap_all_edges(); + /** + * @brief The quality improvement of a swap. + * + * Used to determine the priority and weight of a swap operation. + */ + double swap_weight(const Tuple& t) const; + bool swap_edge_before(const Tuple& t) override; + bool swap_edge_after(const Tuple& t) override; + + void smooth_all_vertices(const size_t n_iters = 1); + bool smooth_before(const Tuple& t) override; + bool smooth_after(const Tuple& t) override; + + /** + * @brief A vector containing the vertex position and all positions of the surface neighbors. + * + * Returns an empty vector if vertex is not on the surface. + */ + std::vector get_surface_assembles(const Tuple& t) const; + std::shared_ptr get_envelope_energy(const Tuple& t) const; + + std::vector> get_amips_assembles(const Tuple& t) const; + std::shared_ptr get_amips_energy(const Tuple& t) const; + + /** + * @brief Inversion check using only floating point numbers. + */ + bool is_inverted_f(const Tuple& loc) const; + bool is_inverted_f(const size_t fid) const; + bool is_inverted(const std::array& vs) const; + bool is_inverted(const Tuple& loc) const; + bool is_inverted(const size_t fid) const; + double get_quality(const std::array& vs) const; + double get_quality(const Tuple& loc) const; + double get_quality(const size_t fid) const; + bool round(const Tuple& loc); + // + bool is_edge_on_surface(const Tuple& loc) const; + bool is_edge_on_surface(const std::array& vids) const; + bool is_edge_on_bbox(const Tuple& loc) const; + bool is_edge_on_bbox(const std::array& vids) const; + + // + void mesh_improvement(int max_its = 80); + + std::tuple local_operations( + const std::array& ops, + bool collapse_limit_length = true); + std::tuple get_max_avg_energy(); + + void compute_winding_numbers(const std::vector& input_paths); + int flood_fill(); + + bool vertex_is_on_surface(const size_t vid) const override + { + return m_vertex_attribute.at(vid).m_is_on_surface || + !m_vertex_attribute.at(vid).on_bbox_faces.empty(); + } + bool edge_is_on_surface(const std::array& vids) const override + { + if (!vertex_is_on_surface(vids[0]) || !vertex_is_on_surface(vids[1])) { + return false; + } + + const auto [_, eid] = tuple_from_edge(vids); + bool on_surface = m_edge_attribute.at(eid).m_is_surface_fs; + bool on_bbox = m_edge_attribute.at(eid).m_is_bbox_fs >= 0; + return on_surface || on_bbox; + } + +private: + ////// Operations + + struct SplitInfoCache + { + // VertexAttributes vertex_info; + size_t v1_id; + size_t v2_id; + std::vector v1_param_type; + std::vector v2_param_type; + + EdgeAttributes old_e_attrs; + + // std::vector>> changed_edges; + std::map changed_edges; + + /** + * All faces incident to the splitted edge, identified by the link vertex (the vertex + * opposite to the splitted edge). + */ + std::map faces; + }; + tbb::enumerable_thread_specific split_cache; + + struct CollapseInfoCache + { + size_t v1_id; + size_t v2_id; + double max_energy; + double edge_length; + bool is_limit_length; + + std::vector>> changed_edges; + // all faces incident to the delete vertex (v1) that are on the tracked surface + std::vector> surface_edges; + std::vector changed_fids; + std::vector changed_energies; + }; + tbb::enumerable_thread_specific collapse_cache; + + + struct SwapInfoCache + { + double max_energy; + std::map changed_edges; + std::set face_tags; + }; + tbb::enumerable_thread_specific swap_cache; +}; + + +} // namespace wmtk::components::triwild diff --git a/components/triwild/wmtk/components/triwild/init_from_delaunay.cpp b/components/triwild/wmtk/components/triwild/init_from_delaunay.cpp new file mode 100644 index 0000000000..ebe0dc78e0 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/init_from_delaunay.cpp @@ -0,0 +1,167 @@ +#include "init_from_delaunay.hpp" + +#include +#include +#include +#include +#include + +namespace wmtk::components::triwild { + +void init_from_delaunay_box_mesh( + const MatrixXd& V, + const MatrixXi& E, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out) +{ + assert(V.cols() == 2); + + // points for delaunay + std::vector points(V.rows()); + // add points from surface + for (int i = 0; i < V.rows(); i++) { + for (int j = 0; j < 2; j++) { + points[i][j] = V(i, j); + } + } + + // bbox + Vector2d box_min = V.colwise().minCoeff(); + Vector2d box_max = V.colwise().maxCoeff(); + + // increase bbox by 30% of diagonal length + const double diagonal_length = (box_max - box_min).norm(); + const double delta = diagonal_length / 15.0; + box_min -= Vector2d(delta, delta); + box_max += Vector2d(delta, delta); + + // add corners of domain + for (int i = 0; i < 4; i++) { + Vector2d p; + std::bitset a(i); + for (int j = 0; j < 2; j++) { + if (a.test(j)) { + p[j] = box_max[j]; + } else { + p[j] = box_min[j]; + } + } + points.push_back({{p[0], p[1]}}); + } + + const double voxel_resolution = diagonal_length / 20.0; + std::array N; // number of grid points per dimension + std::array h; // distance between grid points per dimension + for (int i = 0; i < 2; i++) { + const double D = box_max[i] - box_min[i]; + N[i] = (D / voxel_resolution) + 1; + h[i] = D / N[i]; + } + + std::array, 2> ds; + for (int i = 0; i < 2; i++) { + ds[i].push_back(box_min[i]); + for (int j = 0; j < N[i] - 1; j++) { + ds[i].push_back(box_min[i] + h[i] * (j + 1)); + } + ds[i].push_back(box_max[i]); + } + + wmtk::SampleEnvelope envelope; + // init envelope + { + std::vector V_envelope; + std::vector E_envelope; + for (int i = 0; i < E.rows(); i++) { + E_envelope.push_back(Vector2i(E(i, 0), E(i, 1))); + } + for (int i = 0; i < V.rows(); i++) { + V_envelope.push_back(V.row(i)); + } + envelope.init(V_envelope, E_envelope, 0); + } + + + const double min_dis = voxel_resolution * voxel_resolution / 4; + // double min_dis = state.target_edge_len * state.target_edge_len;//epsilon*2 + for (int i = 0; i < ds[0].size(); i++) { + for (int j = 0; j < ds[1].size(); j++) { + if ((i == 0 || i == ds[0].size() - 1) && (j == 0 || j == ds[1].size() - 1)) { + continue; + } + const Vector2d p(ds[0][i], ds[1][j]); + + Eigen::Vector2d n; + const double sqd = envelope.nearest_point(p, n); + + if (sqd < min_dis) { + continue; + } + points.push_back({{ds[0][i], ds[1][j]}}); + } + } + + // CDT + MatrixXd V_cdt; + V_cdt.resize(points.size(), 2); + for (int i = 0; i < points.size(); i++) { + V_cdt.row(i) = Eigen::Vector2d(points[i][0], points[i][1]); + } + delaunay::constrained_delaunay2D(V_cdt, E, V_out, F_out, E_out); +} + +void init_from_paths( + const std::vector& input_paths, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out) +{ + std::vector Vs; + std::vector Es; + // read input edge meshes + for (const std::string& path : input_paths) { + MatrixXd V; + MatrixXi E; + io::read_edge_mesh(path, V, E); + logger().info("Read edge mesh {}: #V = {}, #E = {}", path, V.rows(), E.rows()); + V = V.block(0, 0, V.rows(), 2).eval(); // only keep x, y + Vs.push_back(V); + Es.push_back(E); + } + if (Vs.size() != Es.size()) { + log_and_throw_error("Vs and Es size mismatch!"); + } + + // generate Delaunay triangulation of the input vertices as the initial mesh + { + MatrixXd V_all; + MatrixXi E_all; + std::vector V_vec; + std::vector E_vec; + for (const auto& V : Vs) { + for (int i = 0; i < V.rows(); i++) { + V_vec.push_back(V.row(i)); + } + } + for (const auto& E : Es) { + for (int i = 0; i < E.rows(); i++) { + E_vec.push_back(E.row(i)); + } + } + V_all.resize(V_vec.size(), 2); + for (int i = 0; i < V_vec.size(); i++) { + V_all.row(i) = V_vec[i]; + } + E_all.resize(E_vec.size(), 2); + for (int i = 0; i < E_vec.size(); i++) { + E_all.row(i) = E_vec[i]; + } + init_from_delaunay_box_mesh(V_all, E_all, V_out, F_out, E_out); + } + + // logger().info("CDT mesh: #V = {}, #F = {}, #E = {}", V_out.rows(), F_out.rows(), + // E_out.rows()); +} + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/init_from_delaunay.hpp b/components/triwild/wmtk/components/triwild/init_from_delaunay.hpp new file mode 100644 index 0000000000..d1464551ba --- /dev/null +++ b/components/triwild/wmtk/components/triwild/init_from_delaunay.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include + +namespace wmtk::components::triwild { + +/** + * @brief Initializes a triangulation from a 2D point set and edge set using constrained Delaunay + * triangulation. The output is a triangulation of the bounding box of the input points that + * contains all edges and vertices from the input. + * + * @param V input vertices (Nx2) + * @param E input edges (Mx2) + * @param V_out output vertices (Kx2) + * @param F_out output faces (Lx3) + * @param E_out output edges (Px2) - the constrained edges from the input + */ +void init_from_delaunay_box_mesh( + const MatrixXd& V, + const MatrixXi& E, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out); + +void init_from_paths( + const std::vector& input_paths, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out); + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/triwild.cpp b/components/triwild/wmtk/components/triwild/triwild.cpp new file mode 100644 index 0000000000..3f26165a30 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/triwild.cpp @@ -0,0 +1,160 @@ +#include "triwild.hpp" + +#include +#include +#include +#include +#include + +#include "Parameters.h" +#include "TriWildMesh.h" +#include "init_from_delaunay.hpp" +#include "triwild_spec.hpp" + +namespace wmtk::components::triwild { + +void triwild(nlohmann::json json_params) +{ + using wmtk::utils::resolve_path; + + // verify input and inject defaults + { + jse::JSE spec_engine; + bool r = spec_engine.verify_json(json_params, triwild_spec); + if (!r) { + log_and_throw_error(spec_engine.log2str()); + } + json_params = spec_engine.inject_defaults(json_params, triwild_spec); + } + const std::filesystem::path root = + json_params.contains("json_input_file") ? json_params["json_input_file"] : ""; + + // logger settings + { + std::string log_file_name = json_params["log_file"]; + if (!log_file_name.empty()) { + log_file_name = resolve_path(root, log_file_name).string(); + wmtk::set_file_logger(log_file_name); + logger().flush_on(spdlog::level::info); + } + } + + std::vector input_paths = json_params["input"]; + for (std::string& p : input_paths) { + p = resolve_path(root, p).string(); + } + + triwild::Parameters params; + + std::string output_path = json_params["output"]; + int NUM_THREADS = json_params["num_threads"]; + int max_its = json_params["max_iterations"]; + + params.epsr = json_params["eps_rel"]; + params.lr = json_params["length_rel"]; + params.stop_energy = json_params["stop_energy"]; + + params.preserve_topology = json_params["preserve_topology"]; + + params.debug_output = json_params["DEBUG_output"]; + params.perform_sanity_checks = json_params["DEBUG_sanity_checks"]; + + // CDT on all input meshes + MatrixXd V; + MatrixXi F; + MatrixXi E; // constraint edges in CDT + init_from_paths(input_paths, V, F, E); + + if (params.debug_output) { + MatrixXd V3(V.rows(), 3); + V3.setZero(); + V3.block(0, 0, V.rows(), 2) = V; + igl::write_triangle_mesh(output_path + "_initial_delaunay.obj", V3, F); + + // write edges + std::ofstream edge_out(output_path + "_initial_edges.obj"); + for (int i = 0; i < E.rows(); i++) { + edge_out << "v " << V(E(i, 0), 0) << " " << V(E(i, 0), 1) << " 0\n"; + edge_out << "v " << V(E(i, 1), 0) << " " << V(E(i, 1), 1) << " 0\n"; + edge_out << "l " << 2 * i + 1 << " " << 2 * i + 2 << "\n"; + } + edge_out.close(); + } + + Vector2d box_min = V.colwise().minCoeff(); + Vector2d box_max = V.colwise().maxCoeff(); + params.init(box_min, box_max); + + TriWildMesh mesh(params, params.eps, NUM_THREADS); + mesh.init_mesh(V, F, E); + + if (params.debug_output) { + mesh.write_vtu(output_path + "_initial"); + } + + /////////mesh improvement + mesh.mesh_improvement(max_its); + mesh.consolidate_mesh(); + + bool all_rounded = true; + for (const auto& v : mesh.get_vertices()) { + if (!mesh.m_vertex_attribute[v.vid(mesh)].m_is_rounded) { + all_rounded = false; + break; + } + } + if (all_rounded) { + logger().info("All vertices are rounded"); + } else { + logger().error("Not all vertices rounded!"); + } + + // apply flood fill + { + int num_parts = mesh.flood_fill(); + logger().info("flood fill parts {}", num_parts); + } + // compute per-input winding number + mesh.compute_winding_numbers(input_paths); + + // double time = timer.getElapsedTime(); + // logger().info("total time {:.4}s", time); + + /////////output + auto [max_energy, avg_energy] = mesh.get_max_avg_energy(); + wmtk::logger().info("final max energy = {} avg = {}", max_energy, avg_energy); + + const std::string report_file = json_params["report"]; + if (!report_file.empty()) { + std::ofstream fout(report_file); + nlohmann::json report; + report["#t"] = mesh.tri_capacity(); + report["#v"] = mesh.vert_capacity(); + report["max_energy"] = max_energy; + report["avg_energy"] = avg_energy; + report["eps"] = params.eps; + report["threads"] = NUM_THREADS; + // report["time"] = time; + report["all_rounded"] = all_rounded; + // report["insertion_and_preprocessing"] = insertion_time; + fout << std::setw(4) << report; + fout.close(); + } + + // check metrics + if (json_params["throw_on_fail"]) { + if (!all_rounded) { + log_and_throw_error("Not all vertices rounded!"); + } + if (max_energy > params.stop_energy) { + log_and_throw_error("Max energy is too large."); + } + } + + mesh.write_vtu(output_path); + mesh.write_msh_groups(output_path + ".msh"); + + logger().info("======= finish ========="); +} + +} // namespace wmtk::components::triwild \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/triwild.hpp b/components/triwild/wmtk/components/triwild/triwild.hpp new file mode 100644 index 0000000000..fa8270c371 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/triwild.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include + +#include "TriWildMesh.h" + + +namespace wmtk::components::triwild { + +void triwild(nlohmann::json json_params); + +} // namespace wmtk::components::triwild diff --git a/components/triwild/wmtk/components/triwild/triwild_spec.hpp b/components/triwild/wmtk/components/triwild/triwild_spec.hpp new file mode 100644 index 0000000000..0c5d3910e2 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/triwild_spec.hpp @@ -0,0 +1,124 @@ +#pragma once +#include +namespace { + +nlohmann::json triwild_spec = R"( +[ + { + "pointer": "/", + "type": "object", + "required": ["application", "input"], + "optional": [ + "output", + "num_threads", + "max_iterations", + "eps_rel", + "length_rel", + "stop_energy", + "preserve_topology", + "throw_on_fail", + "log_file", + "report", + "DEBUG_output", + "DEBUG_sanity_checks" + ] + }, + { + "pointer": "/application", + "type": "string", + "options": ["triwild"], + "doc": "Application name must be triwild." + }, + { + "pointer": "/input", + "type": "list", + "doc": "List of triangular input meshes.", + "min": 1 + }, + { + "pointer": "/input/*", + "type": "string", + "doc": "Triangular input mesh." + }, + { + "pointer": "/output", + "type": "string", + "default": "out", + "doc": "Output file name (without extension)." + }, + { + "pointer": "/num_threads", + "type": "int", + "default": 0, + "doc": "Number of threads used by the application" + }, + { + "pointer": "/max_iterations", + "type": "int", + "default": 80, + "doc": "Maximum iterations before stopping." + }, + { + "pointer": "/eps_rel", + "type": "float", + "default": 2e-3, + "doc": "Envelope thickness relative to the bounding box" + }, + { + "pointer": "/length_rel", + "type": "float", + "default": 5e-2, + "doc": "Target edge length relative to the bounding box" + }, + { + "pointer": "/stop_energy", + "type": "float", + "default": 10, + "doc": "Target energy. If all tets have an energy below this, triwild will stop." + }, + { + "pointer": "/preserve_topology", + "type": "bool", + "default": false, + "doc": "Preserve the topology of the input surface." + }, + { + "pointer": "/throw_on_fail", + "type": "bool", + "default": false, + "doc": "Throw exception if the output does not fulfil the desired criteria. No output will be generated." + }, + { + "pointer": "/log_file", + "type": "string", + "default": "", + "doc": "Logs are not just printed on the terminal but also saved in this file." + }, + { + "pointer": "/report", + "type": "string", + "default": "", + "doc": "A JSON file that stores information about the result and the method execution, e.g., runtime." + }, + { + "pointer": "/DEBUG_output", + "type": "bool", + "default": false, + "doc": "Write the mesh as debug_{}.vtu after every operation." + }, + { + "pointer": "/DEBUG_sanity_checks", + "type": "bool", + "default": false, + "doc": "Perform sanity checks after every operation. This can be very slow and should only be used for debugging." + }, + { + "pointer": "/DEBUG_hausdorff", + "type": "bool", + "default": false, + "doc": "Sanity Check: Compute and report the Hausdorff distance of the output to the input. Should be always smaller than eps." + } +] +)"_json; + +} \ No newline at end of file diff --git a/components/triwild/wmtk/components/triwild/triwild_spec.json b/components/triwild/wmtk/components/triwild/triwild_spec.json new file mode 100644 index 0000000000..8991ea0f87 --- /dev/null +++ b/components/triwild/wmtk/components/triwild/triwild_spec.json @@ -0,0 +1,116 @@ +[ + { + "pointer": "/", + "type": "object", + "required": ["application", "input"], + "optional": [ + "output", + "num_threads", + "max_iterations", + "eps_rel", + "length_rel", + "stop_energy", + "preserve_topology", + "throw_on_fail", + "log_file", + "report", + "DEBUG_output", + "DEBUG_sanity_checks" + ] + }, + { + "pointer": "/application", + "type": "string", + "options": ["triwild"], + "doc": "Application name must be triwild." + }, + { + "pointer": "/input", + "type": "list", + "doc": "List of triangular input meshes.", + "min": 1 + }, + { + "pointer": "/input/*", + "type": "string", + "doc": "Triangular input mesh." + }, + { + "pointer": "/output", + "type": "string", + "default": "out", + "doc": "Output file name (without extension)." + }, + { + "pointer": "/num_threads", + "type": "int", + "default": 0, + "doc": "Number of threads used by the application" + }, + { + "pointer": "/max_iterations", + "type": "int", + "default": 80, + "doc": "Maximum iterations before stopping." + }, + { + "pointer": "/eps_rel", + "type": "float", + "default": 2e-3, + "doc": "Envelope thickness relative to the bounding box" + }, + { + "pointer": "/length_rel", + "type": "float", + "default": 5e-2, + "doc": "Target edge length relative to the bounding box" + }, + { + "pointer": "/stop_energy", + "type": "float", + "default": 10, + "doc": "Target energy. If all tets have an energy below this, triwild will stop." + }, + { + "pointer": "/preserve_topology", + "type": "bool", + "default": false, + "doc": "Preserve the topology of the input surface." + }, + { + "pointer": "/throw_on_fail", + "type": "bool", + "default": false, + "doc": "Throw exception if the output does not fulfil the desired criteria. No output will be generated." + }, + { + "pointer": "/log_file", + "type": "string", + "default": "", + "doc": "Logs are not just printed on the terminal but also saved in this file." + }, + { + "pointer": "/report", + "type": "string", + "default": "", + "doc": "A JSON file that stores information about the result and the method execution, e.g., runtime." + }, + { + "pointer": "/DEBUG_output", + "type": "bool", + "default": false, + "doc": "Write the mesh as debug_{}.vtu after every operation." + }, + { + "pointer": "/DEBUG_sanity_checks", + "type": "bool", + "default": false, + "doc": "Perform sanity checks after every operation. This can be very slow and should only be used for debugging." + }, + { + "pointer": "/DEBUG_hausdorff", + "type": "bool", + "default": false, + "doc": "Sanity Check: Compute and report the Hausdorff distance of the output to the input. Should be always smaller than eps." + } +] diff --git a/src/wmtk/Types.hpp b/src/wmtk/Types.hpp index 180aed124e..be263e7358 100644 --- a/src/wmtk/Types.hpp +++ b/src/wmtk/Types.hpp @@ -48,12 +48,24 @@ using VectorS = Eigen::SparseVector; using VectorSi = VectorS; using VectorSd = VectorS; +inline Vector2r to_rational(const Vector2d& p0) +{ + Vector2r p(p0[0], p0[1]); + return p; +} + inline Vector3r to_rational(const Vector3d& p0) { Vector3r p(p0[0], p0[1], p0[2]); return p; } +inline Vector2d to_double(const Vector2r& p0) +{ + Vector2d p(p0[0].to_double(), p0[1].to_double()); + return p; +} + inline Vector3d to_double(const Vector3r& p0) { Vector3d p(p0[0].to_double(), p0[1].to_double(), p0[2].to_double()); diff --git a/src/wmtk/io/read_edge_mesh.cpp b/src/wmtk/io/read_edge_mesh.cpp new file mode 100644 index 0000000000..7d9331f8ac --- /dev/null +++ b/src/wmtk/io/read_edge_mesh.cpp @@ -0,0 +1,74 @@ +#include "read_edge_mesh.hpp" + +#include +#include +#include +#include +#include +#include + +namespace wmtk::io { + +void read_edge_mesh(const std::string& path, MatrixXd& V, MatrixXi& E) +{ + std::ifstream file(path); + if (!file.is_open()) { + log_and_throw_error("Could not open file: {}", path); + } + + std::vector v_list; + std::vector e_list; + + std::string line; + while (std::getline(file, line)) { + if (line.empty()) { + continue; + } + std::istringstream iss(line); + std::string type; + iss >> type; + + if (type == "v") { + double x = 0, y = 0, z = 0; + iss >> x >> y; + if (iss) { // If it could read 2 + double temp; + if (iss >> temp) { + z = temp; + } + } + v_list.emplace_back(x, y, z); + } else if (type == "l") { + std::string token; + std::vector vs; + while (iss >> token) { + size_t slash_pos = token.find('/'); + std::string v_idx_str = token.substr(0, slash_pos); + if (v_idx_str.empty()) { + continue; + } + + int v_idx = std::stoi(v_idx_str); + if (v_idx < 0) { + v_idx = v_list.size() + v_idx + 1; + } + vs.push_back(v_idx - 1); + } + for (size_t i = 0; i + 1 < vs.size(); i++) { + e_list.emplace_back(vs[i], vs[i + 1]); + } + } + } + + V.resize(v_list.size(), 3); + for (size_t i = 0; i < v_list.size(); i++) { + V.row(i) = v_list[i]; + } + + E.resize(e_list.size(), 2); + for (size_t i = 0; i < e_list.size(); i++) { + E.row(i) = e_list[i]; + } +} + +} // namespace wmtk::io \ No newline at end of file diff --git a/src/wmtk/io/read_edge_mesh.hpp b/src/wmtk/io/read_edge_mesh.hpp new file mode 100644 index 0000000000..10568a0e58 --- /dev/null +++ b/src/wmtk/io/read_edge_mesh.hpp @@ -0,0 +1,21 @@ +#pragma once + +#include + +namespace wmtk::io { + +/** + * @brief Reads an edge mesh from a file. + * The mesh is cleaned by removing duplicated vertices, degenerate faces, and unreferenced vertices. + * + * Duplicated vertices are removed based on the provided tolerances. The tolerance can be specified + * in two ways: absolute and relative. Only one of them can be non-negative. If both are negative, + * duplicated vertices will not be removed. + * + * @param path The file path to read the mesh from. + * @param V Output vertex positions. Size is #V by 3. + * @param E Output edge indices. Size is #E by 2. + */ +void read_edge_mesh(const std::string& path, MatrixXd& V, MatrixXi& E); + +} // namespace wmtk::io \ No newline at end of file diff --git a/src/wmtk/utils/Delaunay.cpp b/src/wmtk/utils/Delaunay.cpp index 8904ebdd93..4c754f3d83 100644 --- a/src/wmtk/utils/Delaunay.cpp +++ b/src/wmtk/utils/Delaunay.cpp @@ -7,6 +7,7 @@ // clang-format on #include +#include #include @@ -181,4 +182,93 @@ auto delaunay2D(const std::vector& points) return {vertices, triangles}; } +void constrained_delaunay2D( + const MatrixXd& V, + const MatrixXi& E, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out) +{ + if (V.cols() != 2) { + log_and_throw_error("Input vertices must be 2D"); + } + if (E.cols() != 2) { + log_and_throw_error("Input edges must be 2D"); + } + + GEO::CDT2d cdt; + // GEO::ExactCDT2d cdt; + std::vector vertex_ids(V.rows()); + + cdt.create_enclosing_rectangle( + V.col(0).minCoeff(), + V.col(1).minCoeff(), + V.col(0).maxCoeff(), + V.col(1).maxCoeff()); + + for (size_t i = 0; i < V.rows(); ++i) { + const GEO::vec2 p(V(i, 0), V(i, 1)); + const GEO::index_t id = cdt.insert(p); + vertex_ids[i] = id; + } + for (size_t i = 0; i < E.rows(); ++i) { + const int v0 = E(i, 0); + const int v1 = E(i, 1); + if (v0 >= V.rows() || v1 >= V.rows()) { + log_and_throw_error("Edge index out of bounds at index {}: {}", i, E.row(i)); + } + cdt.insert_constraint(vertex_ids[v0], vertex_ids[v1]); + } + + GEO::index_t nv = cdt.nv(); + GEO::index_t nt = cdt.nT(); + + // get constrained edges + std::set> constrained_edges; + + for (GEO::index_t t = 0; t < nt; ++t) { + for (GEO::index_t le = 0; le < 3; ++le) { + if (cdt.Tedge_cnstr_first(t, le) == GEO::NO_INDEX) { + continue; + } + + int a = int(cdt.Tv(t, (le + 1) % 3)); + int b = int(cdt.Tv(t, (le + 2) % 3)); + + if (a > b) { + std::swap(a, b); + } + constrained_edges.insert({a, b}); + } + } + + logger().info("CDT: #V = {}, #F = {}, #E_constrained = {}", nv, nt, constrained_edges.size()); + + V_out.resize(nv, 2); + + for (GEO::index_t v = 0; v < nv; ++v) { + const GEO::vec2& p = cdt.point(v); + V_out(v, 0) = p.x; + V_out(v, 1) = p.y; + } + + F_out.resize(nt, 3); + + for (GEO::index_t t = 0; t < nt; ++t) { + F_out(t, 0) = int(cdt.Tv(t, 0)); + F_out(t, 1) = int(cdt.Tv(t, 1)); + F_out(t, 2) = int(cdt.Tv(t, 2)); + } + + E_out.resize(constrained_edges.size(), 2); + { + int idx = 0; + for (const auto& edge : constrained_edges) { + E_out(idx, 0) = edge.first; + E_out(idx, 1) = edge.second; + ++idx; + } + } +} + } // namespace wmtk::delaunay diff --git a/src/wmtk/utils/Delaunay.hpp b/src/wmtk/utils/Delaunay.hpp index 0310771184..f1c020665e 100644 --- a/src/wmtk/utils/Delaunay.hpp +++ b/src/wmtk/utils/Delaunay.hpp @@ -3,6 +3,7 @@ #include #include #include +#include namespace wmtk::delaunay { using Point3D = std::array; @@ -39,4 +40,23 @@ auto delaunay3D(const std::vector& points) auto delaunay2D(const std::vector& points) -> std::pair, std::vector>; +/** + * @brief Compute the constrained Delaunay triangulation of the input 2D points and edges. + * + * The mesh must not have any duplicate vertices or edges. + * + * @param[in] V The input 2D points. + * @param[in] E The constrained input edges. + * @param[out] V_out The output 2D points. + * @param[out] F_out The output triangles. + * @param[out] E_out The constrained input edges in the output mesh. They might differ from the + * input edges if some input edges are split during triangulation. + */ +void constrained_delaunay2D( + const MatrixXd& V, + const MatrixXi& E, + MatrixXd& V_out, + MatrixXi& F_out, + MatrixXi& E_out); + } // namespace wmtk::delaunay diff --git a/tests/integration_tests.cpp b/tests/integration_tests.cpp index 9862fa6754..03550fc18f 100644 --- a/tests/integration_tests.cpp +++ b/tests/integration_tests.cpp @@ -47,7 +47,8 @@ std::vector input_files{ "topological_offset_2d_vertex_input.json", "topological_offset_3d_edge_input.json", "topological_offset_3d.json", - "manifold_extraction_3d.json"}; + "manifold_extraction_3d.json", + "triwild_puzzle.json"}; nlohmann::json load_json(const path& json_input_file)