|
|
|
@ -1,8 +1,19 @@ |
|
|
|
#include <prim_gen.hpp> |
|
|
|
#include <helper.hpp> |
|
|
|
#include <algorithm/find_mismatch_permutation.hpp> |
|
|
|
|
|
|
|
const std::array tet_boundary_faces = {0u, 1u, 2u, 3u}; |
|
|
|
|
|
|
|
void contiguous_hash(size_t& seed, uint32_t value) |
|
|
|
{ |
|
|
|
seed = seed ^ (std::hash<uint32_t>()(value) + 0x9e3779b9 + (seed << 6) + (seed >> 2)); |
|
|
|
} |
|
|
|
|
|
|
|
void contiguous_hash(size_t& seed, size_t value) |
|
|
|
{ |
|
|
|
seed = seed ^ (std::hash<size_t>()(value) + 0x9e3779b9 + (seed << 6) + (seed >> 2)); |
|
|
|
} |
|
|
|
|
|
|
|
void extract_iso_mesh(const std::array<uint32_t, 3>& tet_active_subface_counts, |
|
|
|
const stl_vector_mp<arrangement_t>& tetrahedron_arrangements, |
|
|
|
const scene_bg_mesh_info_t& scene_bg_mesh_info, |
|
|
|
@ -24,10 +35,10 @@ void extract_iso_mesh(const std::array<uint32_t, 3>& tet_active_s |
|
|
|
iso_faces.reserve(max_num_face); |
|
|
|
|
|
|
|
// hash table for vertices on the boundary of tetrahedron
|
|
|
|
flat_hash_map_mp<uint32_t, uint32_t> vert_on_tetVert{}; |
|
|
|
flat_hash_map_mp<pod_key_t<3>, uint32_t> vert_on_tetEdge{}; |
|
|
|
flat_hash_map_mp<uint32_t, uint32_t> vert_on_tetVert{}; |
|
|
|
flat_hash_map_mp<size_t, uint32_t> vert_on_tetEdge{}; |
|
|
|
vert_on_tetEdge.reserve(num_1_func + 3 * num_2_func + num_more_func); |
|
|
|
flat_hash_map_mp<pod_key_t<5>, uint32_t> vert_on_tetFace{}; |
|
|
|
flat_hash_map_mp<size_t, uint32_t> 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<uint32_t, 3>& tet_active_s |
|
|
|
face_verts.reserve(4); |
|
|
|
pod_key_t<3> key3; |
|
|
|
pod_key_t<5> key5; |
|
|
|
std::array<uint16_t, 3> implicit_pIds; |
|
|
|
std::array<bool, 4> used_pId; |
|
|
|
std::array<uint32_t, 2> vIds2; |
|
|
|
std::array<uint32_t, 3> vIds3; |
|
|
|
std::array<uint32_t, 3> implicit_pIds; |
|
|
|
std::array<uint32_t, 3> boundary_pIds; |
|
|
|
std::array<uint32_t, 3> not_boundary_pIds; |
|
|
|
std::array<uint32_t, 4> 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<uint32_t, 3>& 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<uint16_t>(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<uint32_t>(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<uint32_t>(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<uint32_t>(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<uint32_t>(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<uint32_t>(iso_verts.size() - 1)); |
|
|
|
if (iter_inserted.second) { |
|
|
|
auto& iso_vert = iso_verts.emplace_back(); |
|
|
|
iso_vert.simplex_vertex_indices = {static_cast<uint16_t>(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<uint32_t, 3>& 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<uint32_t>(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<uint32_t>(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); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
|