From 8454a60fcc459e09e68b4bbaa2b2d075fae2e447 Mon Sep 17 00:00:00 2001 From: Zhicheng Wang <1627343141@qq.com> Date: Mon, 17 Nov 2025 00:41:45 +0800 Subject: [PATCH] fix known bug; improve these steps --- .../src/connect_by_topo/pair_faces.cpp | 244 ++++++++++-------- .../build_tetrahedron_and_adjacency.cpp | 36 +-- network_process/src/prim_gen/extract_mesh.cpp | 239 +++++++++-------- .../algorithm/find_mismatch_permutation.hpp | 26 ++ 4 files changed, 309 insertions(+), 236 deletions(-) create mode 100644 shared_module/algorithm/find_mismatch_permutation.hpp diff --git a/network_process/src/connect_by_topo/pair_faces.cpp b/network_process/src/connect_by_topo/pair_faces.cpp index bb0236c..be11ad7 100644 --- a/network_process/src/connect_by_topo/pair_faces.cpp +++ b/network_process/src/connect_by_topo/pair_faces.cpp @@ -223,8 +223,7 @@ void pair_patches_in_tets(const stl_vector_mp &containi // map: (tet_id, tet_face_id) -> iso_face_id flat_hash_map_mp iso_face_id_of_face{}; for (const auto &iso_face_id : character_edge.face_indices) { - for (const auto &t : iso_faces[iso_face_id].headers) - iso_face_id_of_face.lazy_emplace(t, [&](const auto &ctor) { ctor(t, iso_face_id); }); + for (const auto &t : iso_faces[iso_face_id].headers) iso_face_id_of_face[t] = iso_face_id; } // find identical tet boundary planes incident to iso-edge // map: (tet_Id, tet_plane_Id) -> (oppo_tet_Id, oppo_tet_plane_Id) @@ -275,120 +274,149 @@ void pair_patches_in_tets(const stl_vector_mp &containi } } } + // find identical faces (tet_id, tet_face_id) on tet boundary planes incident to iso-edge // map: (tet_Id, tet_face_Id) -> (oppo_tet_Id, oppo_tet_face_Id) flat_hash_map_mp opposite_face{}; - // keys of vertices on the boundary of tetrahedron - stl_vector_mp tet_vert_keys{}; - stl_vector_mp> tet_edge_keys{}; - stl_vector_mp> tet_face_keys{}; - // hash table for faces on the boundary of tetrahedron - // map: (i,j,k) -> (tet_Id, tet_face_Id) - flat_hash_map_mp, face_header_t> face_on_tetFace{}; - // auxiliary data - flat_hash_set_mp boundary_vertices{}; - flat_hash_set_mp boundary_faces{}; - // map: local vert index --> boundary vert index - stl_vector_mp boundary_vId_of_vert{}; - stl_vector_mp face_verts{}; - pod_key_t<3> key3; - pod_key_t<5> key5; - std::array implicit_pIds; - std::array boundary_pIds; - std::array not_boundary_pIds; - for (const auto tet_index : containing_tetIds) { - // non-empty tet i - const auto &arrangement = tetrahedron_arrangements[tet_index]; - const auto &vertices = arrangement.vertices; - const auto &faces = arrangement.faces; - auto start_index = tetrahedron_active_subface_start_index[tet_index]; - auto num_func = tetrahedron_active_subface_start_index[tet_index + 1] - start_index; - // find vertices and faces on tet boundary incident to iso-edge - boundary_vertices.clear(); - boundary_faces.clear(); - for (size_t i = 0; i < faces.size(); ++i) { - const auto &face = faces[i]; - if (face.supporting_plane < 4 - && identical_tet_planes.find({tet_index, face.supporting_plane}) != identical_tet_planes.end()) { - boundary_faces.emplace(i); - for (const auto &vertex_index : face.vertices) boundary_vertices.emplace(vertex_index); + { + flat_hash_map_mp face_header_cache{}; + for (const auto tet_index : containing_tetIds) { + const auto &arrangement = tetrahedron_arrangements[tet_index]; + const auto &vertices = arrangement.vertices; + const auto &faces = arrangement.faces; + auto start_index = tetrahedron_active_subface_start_index[tet_index]; + auto num_func = tetrahedron_active_subface_start_index[tet_index + 1] - start_index; + // find faces on tet boundary incident to iso-edge + for (uint32_t face_index = 0; face_index < faces.size(); ++face_index) { + const auto &face = faces[face_index]; + if (face.supporting_plane < 4 + && identical_tet_planes.find({tet_index, face.supporting_plane}) != identical_tet_planes.end()) { + face_header_t face_header{tet_index, face_index}; + auto iso_face_index = iso_face_id_of_face.at(face_header); + auto [iterator, is_new] = face_header_cache.try_emplace(iso_face_index, face_header); + const auto &[_, other_face_header] = *iterator; + if (!is_new) { + // face inserted before, pair the two faces + opposite_face[face_header] = other_face_header; + opposite_face[other_face_header] = face_header; + } + } } } - boundary_vId_of_vert.clear(); - boundary_vId_of_vert.reserve(vertices.size()); - // create boundary vertices - for (const auto &vertex_index : boundary_vertices) { - const auto &vertex = vertices[vertex_index]; - auto &local_vert_to_boundary_vert = boundary_vId_of_vert.emplace_back(invalid_index); + } + // flat_hash_map_mp opposite_face{}; + // // keys of vertices on the boundary of tetrahedron + // stl_vector_mp tet_vert_keys{}; + // stl_vector_mp> tet_edge_keys{}; + // stl_vector_mp> tet_face_keys{}; + // // hash table for faces on the boundary of tetrahedron + // // map: (i,j,k) -> (tet_Id, tet_face_Id) + // flat_hash_map_mp, face_header_t> face_on_tetFace{}; + // // auxiliary data + // flat_hash_set_mp boundary_vertices{}; + // flat_hash_set_mp boundary_faces{}; + // // map: local vert index --> boundary vert index + // stl_vector_mp boundary_vId_of_vert{}; + // stl_vector_mp face_verts{}; + // pod_key_t<3> key3; + // pod_key_t<5> key5; + // std::array implicit_pIds; + // std::array boundary_pIds; + // std::array not_boundary_pIds; + // for (const auto tet_index : containing_tetIds) { + // // non-empty tet i + // const auto &arrangement = tetrahedron_arrangements[tet_index]; + // const auto &vertices = arrangement.vertices; + // const auto &faces = arrangement.faces; + // auto start_index = tetrahedron_active_subface_start_index[tet_index]; + // auto num_func = tetrahedron_active_subface_start_index[tet_index + 1] - start_index; + // // find vertices and faces on tet boundary incident to iso-edge + // boundary_vertices.clear(); + // boundary_faces.clear(); + // for (size_t i = 0; i < faces.size(); ++i) { + // const auto &face = faces[i]; + // if (face.supporting_plane < 4 + // && identical_tet_planes.find({tet_index, face.supporting_plane}) != identical_tet_planes.end()) { + // boundary_faces.emplace(i); + // for (const auto &vertex_index : face.vertices) boundary_vertices.emplace(vertex_index); + // } + // } - uint32_t num_boundary_planes{}; - uint32_t num_impl_planes{}; - // vertex.size() == 3 - for (const auto &vert_plane : vertex) { - if (vert_plane >= 4) { // plane 0,1,2,3 are tet boundaries - implicit_pIds[num_impl_planes++] = active_subface_indices[vert_plane - 4 + start_index]; - } else { - boundary_pIds[num_boundary_planes++] = vert_plane; - } - } - std::sort(boundary_pIds.begin(), boundary_pIds.begin() + num_boundary_planes); - std::set_difference(tet_boundary_faces.begin(), - tet_boundary_faces.end(), - boundary_pIds.begin(), - boundary_pIds.begin() + num_impl_planes, - not_boundary_pIds.begin()); + // boundary_vId_of_vert.clear(); + // boundary_vId_of_vert.resize(vertices.size()); + // // create boundary vertices + // for (const auto &vertex_index : boundary_vertices) { + // const auto &vertex = vertices[vertex_index]; + // auto &local_vert_to_boundary_vert = boundary_vId_of_vert[vertex_index]; - const auto &tet_vertices = tetrahedrons[tet_index]; - if (num_boundary_planes == 1) { // on tet face - key5 = {tet_vertices[not_boundary_pIds[0]], - tet_vertices[not_boundary_pIds[1]], - tet_vertices[not_boundary_pIds[2]], - implicit_pIds[0], - implicit_pIds[1]}; - auto iter = std::find(tet_face_keys.begin(), tet_face_keys.end(), key5); - if (iter != tet_face_keys.end()) - local_vert_to_boundary_vert = std::distance(tet_face_keys.begin(), iter); - else { - local_vert_to_boundary_vert = tet_face_keys.size(); - tet_face_keys.emplace_back(key5); - } - } else if (num_boundary_planes == 2) { // on tet edge - key3 = {tet_vertices[not_boundary_pIds[0]], tet_vertices[not_boundary_pIds[1]], implicit_pIds[0]}; - auto iter = std::find(tet_edge_keys.begin(), tet_edge_keys.end(), key3); - if (iter != tet_edge_keys.end()) - local_vert_to_boundary_vert = std::distance(tet_edge_keys.begin(), iter); - else { - local_vert_to_boundary_vert = tet_edge_keys.size(); - tet_edge_keys.emplace_back(key3); - } - } else if (num_boundary_planes == 3) { // on tet vertex - auto key = tet_vertices[not_boundary_pIds[0]]; - auto iter = std::find(tet_vert_keys.begin(), tet_vert_keys.end(), key); - if (iter != tet_vert_keys.end()) - local_vert_to_boundary_vert = std::distance(tet_vert_keys.begin(), iter); - else { - local_vert_to_boundary_vert = tet_vert_keys.size(); - tet_vert_keys.emplace_back(key); - } - } - } - // pair boundary faces - for (const auto &face_index : boundary_faces) { - const auto &face = faces[face_index]; + // uint32_t num_boundary_planes{}; + // uint32_t num_impl_planes{}; + // // vertex.size() == 3 + // for (const auto &vert_plane : vertex) { + // if (vert_plane >= 4) { // plane 0,1,2,3 are tet boundaries + // implicit_pIds[num_impl_planes++] = active_subface_indices[vert_plane - 4 + start_index]; + // } else { + // boundary_pIds[num_boundary_planes++] = vert_plane; + // } + // } + // std::sort(boundary_pIds.begin(), boundary_pIds.begin() + num_boundary_planes); + // std::set_difference(tet_boundary_faces.begin(), + // tet_boundary_faces.end(), + // boundary_pIds.begin(), + // boundary_pIds.begin() + num_boundary_planes, + // not_boundary_pIds.begin()); - face_verts.clear(); - for (auto vId : face.vertices) face_verts.emplace_back(boundary_vId_of_vert[vId]); - compute_iso_face_key(face_verts, key3); - face_header_t face_header{tet_index, face_index}; - auto iter_inserted = face_on_tetFace.try_emplace(key3, face_header); - if (!iter_inserted.second) { - // face inserted before, pair the two faces - const auto &other_face_header = iter_inserted.first->second; - opposite_face[face_header] = other_face_header; - opposite_face[other_face_header] = face_header; - } - } - } + // const auto &tet_vertices = tetrahedrons[tet_index]; + // if (num_boundary_planes == 1) { // on tet face + // key5 = {tet_vertices[not_boundary_pIds[0]], + // tet_vertices[not_boundary_pIds[1]], + // tet_vertices[not_boundary_pIds[2]], + // implicit_pIds[0], + // implicit_pIds[1]}; + // auto iter = std::find(tet_face_keys.begin(), tet_face_keys.end(), key5); + // if (iter != tet_face_keys.end()) + // local_vert_to_boundary_vert = std::distance(tet_face_keys.begin(), iter); + // else { + // local_vert_to_boundary_vert = tet_face_keys.size(); + // tet_face_keys.emplace_back(key5); + // } + // } else if (num_boundary_planes == 2) { // on tet edge + // key3 = {tet_vertices[not_boundary_pIds[0]], tet_vertices[not_boundary_pIds[1]], implicit_pIds[0]}; + // auto iter = std::find(tet_edge_keys.begin(), tet_edge_keys.end(), key3); + // if (iter != tet_edge_keys.end()) + // local_vert_to_boundary_vert = std::distance(tet_edge_keys.begin(), iter); + // else { + // local_vert_to_boundary_vert = tet_edge_keys.size(); + // tet_edge_keys.emplace_back(key3); + // } + // } else if (num_boundary_planes == 3) { // on tet vertex + // auto key = tet_vertices[not_boundary_pIds[0]]; + // auto iter = std::find(tet_vert_keys.begin(), tet_vert_keys.end(), key); + // if (iter != tet_vert_keys.end()) + // local_vert_to_boundary_vert = std::distance(tet_vert_keys.begin(), iter); + // else { + // local_vert_to_boundary_vert = tet_vert_keys.size(); + // tet_vert_keys.emplace_back(key); + // } + // } + // } + // // pair boundary faces + // for (const auto &face_index : boundary_faces) { + // const auto &face = faces[face_index]; + + // face_verts.clear(); + // for (auto vId : face.vertices) face_verts.emplace_back(boundary_vId_of_vert[vId]); + // compute_iso_face_key(face_verts, key3); + // face_header_t face_header{tet_index, face_index}; + // auto iter_inserted = face_on_tetFace.try_emplace(key3, face_header); + // if (!iter_inserted.second) { + // // face inserted before, pair the two faces + // const auto &other_face_header = iter_inserted.first->second; + // opposite_face[face_header] = other_face_header; + // opposite_face[other_face_header] = face_header; + // } + // } + // } // find the half-iso-face of the local face in tet with given orientation // the orientation of an iso-face is defined by the smallest-index implicit function passing the iso-face diff --git a/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp index 9079aeb..2f28386 100644 --- a/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp +++ b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp @@ -35,19 +35,20 @@ void build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t& for (size_t i = 0; i < pos_hash_max(culled_scene_bg_mesh_info); i += 2) { auto [x, y, z] = get_grid_pos(i, culled_scene_bg_mesh_info); auto v0 = hash_vert_pos(x, y, z, scene_bg_mesh_info); - auto v1 = hash_increment_x(v0, scene_bg_mesh_info); - auto v2 = hash_increment_y(v1, scene_bg_mesh_info); - auto v3 = hash_increment_y(v0, scene_bg_mesh_info); - auto v4 = hash_increment_z(v0, scene_bg_mesh_info); - auto v5 = hash_increment_x(v4, scene_bg_mesh_info); - auto v6 = hash_increment_y(v5, scene_bg_mesh_info); - auto v7 = hash_increment_y(v4, scene_bg_mesh_info); + auto v1 = hash_increment_x(v0, scene_bg_mesh_info); // +yz + auto v2 = hash_increment_y(v1, scene_bg_mesh_info); // +yz+z + auto v3 = hash_increment_y(v0, scene_bg_mesh_info); // +z + auto v4 = hash_increment_z(v0, scene_bg_mesh_info); // +1 + auto v5 = hash_increment_x(v4, scene_bg_mesh_info); // +yz+1 + auto v6 = hash_increment_y(v5, scene_bg_mesh_info); // +yz+z+1 + auto v7 = hash_increment_y(v4, scene_bg_mesh_info); // +z+1 + // vertex lexigraphical ordering: v0, v4, v3, v1, v7, v5, v2, v6 - tetrahedrons.emplace_back(std::array{v1, v3, v4, v6}); - tetrahedrons.emplace_back(std::array{v3, v4, v6, v7}); - tetrahedrons.emplace_back(std::array{v0, v1, v3, v4}); - tetrahedrons.emplace_back(std::array{v1, v2, v3, v6}); - tetrahedrons.emplace_back(std::array{v1, v4, v5, v6}); + tetrahedrons.emplace_back(std::array{v4, v3, v1, v6}); + tetrahedrons.emplace_back(std::array{v4, v3, v7, v6}); + tetrahedrons.emplace_back(std::array{v0, v4, v3, v1}); + tetrahedrons.emplace_back(std::array{v3, v1, v2, v6}); + tetrahedrons.emplace_back(std::array{v4, v1, v5, v6}); insert_or_compare_vertex_adjacency(v1, v0); insert_or_compare_vertex_adjacency(v2, v3); @@ -68,12 +69,13 @@ void build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t& auto v5 = hash_increment_x(v4, scene_bg_mesh_info); auto v6 = hash_increment_y(v5, scene_bg_mesh_info); auto v7 = hash_increment_y(v4, scene_bg_mesh_info); + // vertex lexigraphical ordering: v0, v4, v3, v1, v7, v5, v2, v6 - tetrahedrons.emplace_back(std::array{v0, v2, v5, v7}); - tetrahedrons.emplace_back(std::array{v0, v2, v3, v7}); - tetrahedrons.emplace_back(std::array{v0, v4, v5, v7}); - tetrahedrons.emplace_back(std::array{v2, v5, v6, v7}); - tetrahedrons.emplace_back(std::array{v0, v1, v2, v5}); + tetrahedrons.emplace_back(std::array{v0, v7, v5, v2}); + tetrahedrons.emplace_back(std::array{v0, v3, v7, v2}); + tetrahedrons.emplace_back(std::array{v0, v4, v7, v5}); + tetrahedrons.emplace_back(std::array{v7, v5, v2, v6}); + tetrahedrons.emplace_back(std::array{v0, v1, v5, v2}); insert_or_compare_vertex_adjacency(v1, v0); insert_or_compare_vertex_adjacency(v2, v0); diff --git a/network_process/src/prim_gen/extract_mesh.cpp b/network_process/src/prim_gen/extract_mesh.cpp index 37bf55f..b9ff0ca 100644 --- a/network_process/src/prim_gen/extract_mesh.cpp +++ b/network_process/src/prim_gen/extract_mesh.cpp @@ -1,8 +1,19 @@ #include #include +#include const std::array tet_boundary_faces = {0u, 1u, 2u, 3u}; +void contiguous_hash(size_t& seed, uint32_t value) +{ + seed = seed ^ (std::hash()(value) + 0x9e3779b9 + (seed << 6) + (seed >> 2)); +} + +void contiguous_hash(size_t& seed, size_t value) +{ + seed = seed ^ (std::hash()(value) + 0x9e3779b9 + (seed << 6) + (seed >> 2)); +} + void extract_iso_mesh(const std::array& tet_active_subface_counts, const stl_vector_mp& tetrahedron_arrangements, const scene_bg_mesh_info_t& scene_bg_mesh_info, @@ -24,10 +35,10 @@ void extract_iso_mesh(const std::array& tet_active_s iso_faces.reserve(max_num_face); // hash table for vertices on the boundary of tetrahedron - flat_hash_map_mp vert_on_tetVert{}; - flat_hash_map_mp, uint32_t> vert_on_tetEdge{}; + flat_hash_map_mp vert_on_tetVert{}; + flat_hash_map_mp vert_on_tetEdge{}; vert_on_tetEdge.reserve(num_1_func + 3 * num_2_func + num_more_func); - flat_hash_map_mp, uint32_t> vert_on_tetFace{}; + flat_hash_map_mp vert_on_tetFace{}; vert_on_tetFace.reserve(num_2_func + 6 * num_more_func); // hash table for faces on the boundary of tetrahedron @@ -42,9 +53,12 @@ void extract_iso_mesh(const std::array& tet_active_s face_verts.reserve(4); pod_key_t<3> key3; pod_key_t<5> key5; - std::array implicit_pIds; + std::array used_pId; + std::array vIds2; + std::array vIds3; + std::array implicit_pIds; std::array boundary_pIds; - std::array not_boundary_pIds; + std::array not_boundary_vIds{}; // for (uint32_t tet_index = 0; tet_index < tetrahedrons.size(); tet_index++) { @@ -63,131 +77,132 @@ void extract_iso_mesh(const std::array& tet_active_s const auto& face = faces[i]; if (face.supporting_plane > 3) { // plane 0,1,2,3 are tet boundaries tet_iso_faces.emplace(i); - for (const auto& vId : face.vertices) tet_iso_verts.emplace(vId); + tet_iso_verts.insert(face.vertices.begin(), face.vertices.end()); } } } else { for (size_t i = 0; i < faces.size(); ++i) { - const auto& face = faces[i]; - const auto pId = face.supporting_plane; - const auto uId = arrangement.unique_plane_indices[pId]; - for (const auto& plane_index : arrangement.unique_planes[uId]) { - if (plane_index > 3) { // plane 0,1,2,3 are tet boundaries - tet_iso_faces.emplace(i); - for (const auto& vId : face.vertices) tet_iso_verts.emplace(vId); - } + const auto& face = faces[i]; + const auto pId = face.supporting_plane; + const auto uId = arrangement.unique_plane_indices[pId]; + // plane 0,1,2,3 are tet boundaries + auto implicit_plane_iter = std::find_if(arrangement.unique_planes[uId].begin(), + arrangement.unique_planes[uId].end(), + [](uint32_t plane_index) { return plane_index > 3; }); + if (implicit_plane_iter != arrangement.unique_planes[uId].end()) { + tet_iso_faces.emplace(*implicit_plane_iter); + tet_iso_verts.insert(face.vertices.begin(), face.vertices.end()); } } } iso_vId_of_vert.clear(); iso_vId_of_vert.resize(vertices.size()); + const auto& tet_vertices = tetrahedrons[tet_index]; // create iso-vertices for (const auto& vertex_index : tet_iso_verts) { const auto& vertex = vertices[vertex_index]; auto& local_vert_to_iso_vert = iso_vId_of_vert[vertex_index]; - uint32_t num_boundary_planes{}; - uint32_t num_impl_planes{}; + auto boundary_plane_end = boundary_pIds.begin(); + auto impl_plane_end = implicit_pIds.begin(); for (const auto& vertex_plane : vertex) { if (vertex_plane >= 4) { // plane 0,1,2,3 are tet boundaries - implicit_pIds[num_impl_planes++] = - static_cast(active_subface_indices[vertex_plane - 4 + start_index]); + const auto impl_plane = active_subface_indices[vertex_plane - 4 + start_index]; + *impl_plane_end = impl_plane; + ++impl_plane_end; } else { - boundary_pIds[num_boundary_planes++] = vertex_plane; + *boundary_plane_end = vertex_plane; + ++boundary_plane_end; } } - std::sort(boundary_pIds.begin(), boundary_pIds.begin() + num_boundary_planes); - std::set_difference(tet_boundary_faces.begin(), - tet_boundary_faces.end(), - boundary_pIds.begin(), - boundary_pIds.begin() + num_boundary_planes, - not_boundary_pIds.begin()); - - const auto& tet_vertices = tetrahedrons[tet_index]; - vertex_header_t iso_vert_header{tet_index, vertex_index, 4 - num_boundary_planes}; + not_boundary_vIds.fill(invalid_index); + auto end_iter = algorithm::find_mismatch_permutation(boundary_pIds.begin(), + boundary_plane_end, + size_t{4}, + not_boundary_vIds.begin()); + std::transform(not_boundary_vIds.begin(), end_iter, not_boundary_vIds.begin(), [&](uint32_t local_index) { + return tet_vertices[local_index]; + }); + size_t hash_key = + XXH3_64bits(implicit_pIds.data(), sizeof(uint32_t) * std::distance(implicit_pIds.begin(), impl_plane_end)); + for (auto iter = not_boundary_vIds.begin(); iter != end_iter - 1; ++iter) { + assert(*iter < *(iter + 1)); + contiguous_hash(hash_key, *iter); + } + contiguous_hash(hash_key, *(end_iter - 1)); + + uint32_t num_boundary_planes = std::distance(boundary_pIds.begin(), boundary_plane_end); + const vertex_header_t vertex_header{tet_index, vertex_index, 4 - num_boundary_planes}; + auto mapped_not_boundary_vIds = not_boundary_vIds; + std::transform(mapped_not_boundary_vIds.begin(), + mapped_not_boundary_vIds.begin() + (end_iter - not_boundary_vIds.begin()), + mapped_not_boundary_vIds.begin(), + [&](uint32_t vId) { return vertex_indices_mapping.at(vId); }); + + auto generate_new_iso_vert = [&](auto&& vert_hash, auto&& f) { + const auto& [_, iso_vert_index] = *vert_hash.lazy_emplace(hash_key, [&](auto&& ctor) { + ctor(hash_key, static_cast(iso_verts.size())); + + auto& iso_vert = iso_verts.emplace_back(); + iso_vert.header = vertex_header; + iso_vert.simplex_vertex_indices = not_boundary_vIds; + + f(); + }); + local_vert_to_iso_vert = iso_vert_index; + }; + if (num_boundary_planes == 0) { // in tet cell local_vert_to_iso_vert = static_cast(iso_verts.size()); auto& iso_vert = iso_verts.emplace_back(); - iso_vert.simplex_vertex_indices = tet_vertices; - iso_vert.subface_indices = implicit_pIds; - iso_vert.header = iso_vert_header; - - const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[0])), - vert_face1 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[1])), - vert_face2 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[2])), - vert_face3 = vertex_infos.row(vertex_indices_mapping.at(tet_vertices[3])); - const auto f1 = std::array{vert_face0(implicit_pIds[0]), - vert_face1(implicit_pIds[0]), - vert_face2(implicit_pIds[0]), - vert_face3(implicit_pIds[0])}; - const auto f2 = std::array{vert_face0(implicit_pIds[1]), - vert_face1(implicit_pIds[1]), - vert_face2(implicit_pIds[1]), - vert_face3(implicit_pIds[1])}; - const auto f3 = std::array{vert_face0(implicit_pIds[2]), - vert_face1(implicit_pIds[2]), - vert_face2(implicit_pIds[2]), - vert_face3(implicit_pIds[2])}; - const auto coord = compute_barycentric_coords(f1, f2, f3); - iso_pts.emplace_back(coord[0] * get_vert_pos(tet_vertices[0], scene_bg_mesh_info) - + coord[1] * get_vert_pos(tet_vertices[1], scene_bg_mesh_info) - + coord[2] * get_vert_pos(tet_vertices[2], scene_bg_mesh_info) - + coord[3] * get_vert_pos(tet_vertices[3], scene_bg_mesh_info)); + iso_vert.header = vertex_header; + iso_vert.simplex_vertex_indices = not_boundary_vIds; + + const auto &impl0 = vertex_infos.col(implicit_pIds[0]), &impl1 = vertex_infos.col(implicit_pIds[1]), + &impl2 = vertex_infos.col(implicit_pIds[2]); + const std::array f1 = {impl0[mapped_not_boundary_vIds[0]], + impl0[mapped_not_boundary_vIds[1]], + impl0[mapped_not_boundary_vIds[2]], + impl0[mapped_not_boundary_vIds[3]]}, + f2 = {impl1[mapped_not_boundary_vIds[0]], + impl1[mapped_not_boundary_vIds[1]], + impl1[mapped_not_boundary_vIds[2]], + impl1[mapped_not_boundary_vIds[3]]}, + f3 = {impl2[mapped_not_boundary_vIds[0]], + impl2[mapped_not_boundary_vIds[1]], + impl2[mapped_not_boundary_vIds[2]], + impl2[mapped_not_boundary_vIds[3]]}; + const auto coord = compute_barycentric_coords(f1, f2, f3); + iso_pts.emplace_back(coord[0] * get_vert_pos(not_boundary_vIds[0], scene_bg_mesh_info) + + coord[1] * get_vert_pos(not_boundary_vIds[1], scene_bg_mesh_info) + + coord[2] * get_vert_pos(not_boundary_vIds[2], scene_bg_mesh_info) + + coord[3] * get_vert_pos(not_boundary_vIds[3], scene_bg_mesh_info)); } else if (num_boundary_planes == 1) { // on tet face - key5 = {tet_vertices[not_boundary_pIds[0]], - tet_vertices[not_boundary_pIds[1]], - tet_vertices[not_boundary_pIds[2]], - implicit_pIds[0], - implicit_pIds[1]}; - auto iter_inserted = vert_on_tetFace.try_emplace(key5, static_cast(iso_verts.size())); - if (iter_inserted.second) { - auto& iso_vert = iso_verts.emplace_back(); - iso_vert.simplex_vertex_indices = {key5[0], key5[1], key5[2]}; - iso_vert.subface_indices = {implicit_pIds[0], implicit_pIds[1]}; - iso_vert.header = iso_vert_header; - - // - const auto &vert_face0 = vertex_infos.row(vertex_indices_mapping.at(key5[0])), - vert_face1 = vertex_infos.row(vertex_indices_mapping.at(key5[1])), - vert_face2 = vertex_infos.row(vertex_indices_mapping.at(key5[2])); - const auto f1 = - std::array{vert_face0(implicit_pIds[0]), vert_face1(implicit_pIds[0]), vert_face2(implicit_pIds[0])}; - const auto f2 = - std::array{vert_face0(implicit_pIds[1]), vert_face1(implicit_pIds[1]), vert_face2(implicit_pIds[1])}; - const auto coord = compute_barycentric_coords(f1, f2); - iso_pts.emplace_back(coord[0] * get_vert_pos(key5[0], scene_bg_mesh_info) - + coord[1] * get_vert_pos(key5[1], scene_bg_mesh_info) - + coord[2] * get_vert_pos(key5[2], scene_bg_mesh_info)); - } - local_vert_to_iso_vert = iter_inserted.first->second; + generate_new_iso_vert(vert_on_tetFace, [&] { + const auto & impl0 = vertex_infos.col(implicit_pIds[0]), &impl1 = vertex_infos.col(implicit_pIds[1]); + const std::array f1 = {impl0[mapped_not_boundary_vIds[0]], + impl0[mapped_not_boundary_vIds[1]], + impl0[mapped_not_boundary_vIds[2]]}, + f2 = {impl1[mapped_not_boundary_vIds[0]], + impl1[mapped_not_boundary_vIds[1]], + impl1[mapped_not_boundary_vIds[2]]}; + const auto coord = compute_barycentric_coords(f1, f2); + iso_pts.emplace_back(coord[0] * get_vert_pos(not_boundary_vIds[0], scene_bg_mesh_info) + + coord[1] * get_vert_pos(not_boundary_vIds[1], scene_bg_mesh_info) + + coord[2] * get_vert_pos(not_boundary_vIds[2], scene_bg_mesh_info)); + }); } else if (num_boundary_planes == 2) { // on tet edge - key3 = {tet_vertices[not_boundary_pIds[0]], tet_vertices[not_boundary_pIds[1]], implicit_pIds[0]}; - auto iter_inserted = vert_on_tetEdge.try_emplace(key3, static_cast(iso_verts.size() - 1)); - if (iter_inserted.second) { - auto& iso_vert = iso_verts.emplace_back(); - iso_vert.simplex_vertex_indices = {key3[0], key3[1]}; - iso_vert.subface_indices = {implicit_pIds[0]}; - iso_vert.header = iso_vert_header; - - const auto f1 = vertex_infos(vertex_indices_mapping.at(key3[0]), implicit_pIds[0]); - const auto f2 = vertex_infos(vertex_indices_mapping.at(key3[1]), implicit_pIds[0]); + generate_new_iso_vert(vert_on_tetEdge, [&] { + const auto f1 = vertex_infos(vertex_indices_mapping.at(not_boundary_vIds[0]), implicit_pIds[0]); + const auto f2 = vertex_infos(vertex_indices_mapping.at(not_boundary_vIds[1]), implicit_pIds[0]); const auto coord = compute_barycentric_coords(f1, f2); - iso_pts.emplace_back(coord[0] * get_vert_pos(key3[0], scene_bg_mesh_info) - + coord[1] * get_vert_pos(key3[1], scene_bg_mesh_info)); - } - local_vert_to_iso_vert = iter_inserted.first->second; + iso_pts.emplace_back(coord[0] * get_vert_pos(not_boundary_vIds[0], scene_bg_mesh_info) + + coord[1] * get_vert_pos(not_boundary_vIds[1], scene_bg_mesh_info)); + }); } else { // on tet vertex - auto key = tet_vertices[not_boundary_pIds[0]]; - auto iter_inserted = vert_on_tetVert.try_emplace(key, static_cast(iso_verts.size() - 1)); - if (iter_inserted.second) { - auto& iso_vert = iso_verts.emplace_back(); - iso_vert.simplex_vertex_indices = {static_cast(key)}; - iso_vert.header = iso_vert_header; - - iso_pts.emplace_back(get_vert_pos(iso_vert.simplex_vertex_indices[0], scene_bg_mesh_info)); - } - local_vert_to_iso_vert = iter_inserted.first->second; + generate_new_iso_vert(vert_on_tetVert, + [&] { iso_pts.emplace_back(get_vert_pos(not_boundary_vIds[0], scene_bg_mesh_info)); }); } } // create iso-faces @@ -202,18 +217,20 @@ void extract_iso_mesh(const std::array& tet_active_s bool is_new_iso_face{true}; if (face.negative_cell == invalid_index) { compute_iso_face_key(face_verts, key3); - auto [iter, is_new] = face_on_tetFace.try_emplace(key3, static_cast(iso_faces.size())); - if (!is_new) { - iso_face_index = iter->second; - is_new_iso_face = false; - } - } - if (is_new_iso_face) { + const auto& [_, iso_face_index] = *face_on_tetFace.lazy_emplace(key3, [&](auto&& ctor) { + ctor(key3, static_cast(iso_faces.size())); + + auto& iso_face = iso_faces.emplace_back(); + iso_face.vertex_indices = face_verts; + iso_face.subface_index = active_subface_indices[face.supporting_plane - 4 + start_index]; + }); + iso_faces[iso_face_index].headers.emplace_back(face_header); + } else { // face not on tet boundary auto& iso_face = iso_faces.emplace_back(); iso_face.vertex_indices = face_verts; iso_face.subface_index = active_subface_indices[face.supporting_plane - 4 + start_index]; + iso_face.headers.emplace_back(face_header); } - iso_faces[iso_face_index].headers.emplace_back(face_header); } } diff --git a/shared_module/algorithm/find_mismatch_permutation.hpp b/shared_module/algorithm/find_mismatch_permutation.hpp new file mode 100644 index 0000000..72cd74e --- /dev/null +++ b/shared_module/algorithm/find_mismatch_permutation.hpp @@ -0,0 +1,26 @@ +#include + +namespace algorithm +{ +template +OutputIt find_mismatch_permutation(InputIt first, InputIt last, size_t range, OutputIt dest) +{ + using input_value_type = typename std::iterator_traits::value_type; + using output_value_type = typename std::iterator_traits::value_type; + static_assert(std::is_integral_v, "InputIt must point to an integral type"); + static_assert(std::is_integral_v, "OutputIt must point to an integral type"); + static_assert(std::is_convertible_v, + "InputIt value type must be convertible to OutputIt value type"); + using value_type = typename std::iterator_traits::value_type; + auto dst = dest; + for (value_type i = 0; i < range; ++i) { + auto iter = std::find(first, last, i); + if (iter == last) { + *dst = i; + dst++; + } + } + + return dst; +} +} // namespace algorithm \ No newline at end of file