diff --git a/components/wmtk_components/CMakeLists.txt b/components/wmtk_components/CMakeLists.txt index 014a04d88e..0992568b0f 100644 --- a/components/wmtk_components/CMakeLists.txt +++ b/components/wmtk_components/CMakeLists.txt @@ -9,4 +9,5 @@ add_subdirectory(input) add_subdirectory(isotropic_remeshing) add_subdirectory(mesh_info) add_subdirectory(output) +add_subdirectory(extract_subset) add_subdirectory(regular_space) \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/CMakeLists.txt b/components/wmtk_components/extract_subset/CMakeLists.txt new file mode 100644 index 0000000000..b67c632b25 --- /dev/null +++ b/components/wmtk_components/extract_subset/CMakeLists.txt @@ -0,0 +1,19 @@ +set(SRC_FILES + #internal/extract_subset_2d.hpp + #internal/extract_subset_2d.cpp + #internal/topology_separate_2d.cpp + #internal/topology_separate_2d.hpp + #internal/extract_subset_3d.hpp + #internal/extract_subset_3d.cpp + #internal/topology_separate_3d.cpp + #internal/topology_separate_3d.hpp + internal/utils.hpp + internal/utils.cpp + internal/generate_submesh.cpp + internal/generate_submesh.hpp + internal/new_topology_separate.cpp + internal/new_topology_separate.hpp + extract_subset.hpp + extract_subset.cpp + ) +target_sources(wildmeshing_components PRIVATE ${SRC_FILES}) \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/extract_subset.cpp b/components/wmtk_components/extract_subset/extract_subset.cpp new file mode 100644 index 0000000000..e411cba0fc --- /dev/null +++ b/components/wmtk_components/extract_subset/extract_subset.cpp @@ -0,0 +1,44 @@ +#include "extract_subset.hpp" +namespace wmtk { +namespace components { + +Eigen::VectorX& vector2tag(Eigen::VectorX& ret, std::vector vector) +{ + ret.resize(vector.size()); + for (int i = 0; i < vector.size(); ++i) ret.row(i) << vector[i]; + return ret; +} + +std::unique_ptr extract_subset(wmtk::Mesh& m, const std::vector& tag_vec, bool pos) +{ + wmtk::PrimitiveType topType = m.top_simplex_type(); + // tag vector must have the same size as the number of simplices in the mesh + assert(tag_vec.size() == m.get_all(topType).size()); + if (pos) { // if user asks to preserve geometry, then geometry must be provided + try { + m.get_attribute_handle("position", wmtk::PrimitiveType::Vertex); + } catch (const std::exception& e) { + throw std::runtime_error("input mesh doesn't have position attributes!"); + } + } + + Eigen::VectorX tag; + wmtk::MeshAttributeHandle tag_handle = + wmtk::mesh_utils::set_matrix_attribute(vector2tag(tag, tag_vec), "tag", topType, m); + + std::unique_ptr submesh; // Declare the submesh variable here + switch (m.top_cell_dimension()) { + case 2: + case 3: + submesh = internal::generate_submesh( + m, + tag_handle, + pos); // Assign the value inside the switch statement + // return submesh; + return internal::topology_separate(*(submesh.get()), pos); + default: throw std::runtime_error("Invalid mesh dimension in extracting subset!"); + } +} + +} // namespace components +} // namespace wmtk \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/extract_subset.hpp b/components/wmtk_components/extract_subset/extract_subset.hpp new file mode 100644 index 0000000000..a11915329b --- /dev/null +++ b/components/wmtk_components/extract_subset/extract_subset.hpp @@ -0,0 +1,17 @@ +#pragma once + +#include "internal/generate_submesh.hpp" +#include "internal/new_topology_separate.hpp" +namespace wmtk { + +namespace components { + +/* +This function provides a unified interface for extracting a subset of a mesh, 2d or 3d, with or +without preserving geometry. +*/ +std::unique_ptr +extract_subset(wmtk::Mesh& m, const std::vector& tag_vec, bool pos); +// wmtk::Mesh& extract_subset(wmtk::Mesh& m, const std::vector& tag_vec, bool pos); +} // namespace components +} // namespace wmtk \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.cpp b/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.cpp new file mode 100644 index 0000000000..468bd07873 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.cpp @@ -0,0 +1,105 @@ +#include "extract_subset_2d.hpp" +#include +namespace wmtk::components::internal { + +wmtk::TriMesh +extract_subset_2d(wmtk::TriMesh m, wmtk::MeshAttributeHandle tag_handle, bool pos) +{ + wmtk::Accessor tag_acc = m.create_accessor(tag_handle); + std::vector faces = m.get_all(wmtk::PrimitiveType::Face); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + int nb_vertex = vertices.size(); + int nb_tri = faces.size(); + + // a tag on each "real" vertex, true if tagged inside + std::map vertices_in_bool; + for (long t = 0; t < nb_vertex; ++t) vertices_in_bool.insert({t, false}); + + // both init to 0, increment by count later + long nb_vertex_in = 0, nb_tri_in = 0; + + // store the temporary "id" of the tagged triangles + std::vector tag_tri_index; + for (size_t i = 0; i < nb_tri; ++i) { + long tri_tag = tag_acc.const_scalar_attribute(faces.at(i)); + switch (tri_tag) { + // inside: store the temp id of this tri + case 1: + tag_tri_index.push_back(i); + break; + // outside: do nothing + case 0: break; + // neither: runtime error + default: throw std::runtime_error("illegal tag!"); + } + } + nb_tri_in = tag_tri_index.size(); + assert(nb_tri_in <= nb_tri); + // std::cout << "# of tri inside = " << nb_tri_in << std::endl; + + // for the tagged tri, mark their "real" vertices as inside (duplicates handled by boolean) + // current algo for bug fixing: O(N^2), go over all vertices and look for match, + // only assign tag to inside ones + // TODO: improve the algorithm to achieve O(N) + for (size_t i = 0; i < nb_tri_in; ++i) { + Simplex s = Simplex::face(faces[tag_tri_index[i]]); + std::vector tuple_list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + for (wmtk::Tuple t : tuple_list) vertices_in_bool[find_vertex_index(m, t)] = true; + } + + // std::cout << "# of vertex inside = " << vertices_in_bool.size() << std::endl; + // construct a map from old tuple to temp new "id" of a "real" vertex + std::map old2new; + for (long i = 0; i < nb_vertex; ++i) { + if (vertices_in_bool[i]) { + // std::cout << "inside! nb_vertex_in = " << nb_vertex_in << std::endl; + // old vertex tuple t mapped to new vertex id j, where j increases by count + old2new.insert({i, nb_vertex_in}); + nb_vertex_in++; + } + } + // std::cout << "# of pairs in map = " << old2new.size() << std::endl; + // for (auto [_, i] : old2new) std::cout << i << std::endl; + + wmtk::TriMesh mesh; + wmtk::RowVectors3l tris; + tris.resize(nb_tri_in, 3); + // only put in the extracted ones + for (size_t i = 0; i < nb_tri_in; ++i) { + Simplex s = Simplex::face(faces[tag_tri_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(3, -1); + for (int index = 0; index < 3; ++index) + data[index] = old2new[find_vertex_index(m, list[index])]; + tris.row(i) << data[0], data[1], data[2]; + } + // for (size_t i = 0; i < nb_tri_in; ++i) { + // std::cout << tris.row(i)[0] << tris.row(i)[1] << tris.row(i)[2] << std::endl; + // } + mesh.initialize(tris); // init the topology + + // if told to extract and preserve the coordinates + if (pos) { + Eigen::MatrixXd points_in; + points_in.resize(nb_vertex_in, 2); + wmtk::MeshAttributeHandle pos_handle = + m.get_attribute_handle("position", PrimitiveType::Vertex); + wmtk::ConstAccessor pos_acc = m.create_const_accessor(pos_handle); + for (const Tuple& t : vertices) { + // ignore the outside vertices + long old_index = find_vertex_index(m, t); + if (vertices_in_bool[old_index]) { + points_in.row(old2new[old_index]) = pos_acc.const_vector_attribute(t); + } + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + mesh); + } + } + return mesh; +} +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.hpp b/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.hpp new file mode 100644 index 0000000000..392bc3eb23 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/extract_subset_2d.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { + +wmtk::TriMesh +extract_subset_2d(wmtk::TriMesh m, wmtk::MeshAttributeHandle taghandle, bool pos); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.cpp b/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.cpp new file mode 100644 index 0000000000..f07a0b15c1 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.cpp @@ -0,0 +1,97 @@ +#include "extract_subset_3d.hpp" + +namespace wmtk::components::internal { +wmtk::TetMesh +extract_subset_3d(wmtk::TetMesh m, wmtk::MeshAttributeHandle taghandle, bool pos) +{ + wmtk::Accessor tag_acc = m.create_accessor(taghandle); + std::vector tets = m.get_all(wmtk::PrimitiveType::Tetrahedron); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + int nb_vertex = vertices.size(); + int nb_tet = tets.size(); + + std::map vertices_in_bool; + for (long t = 0; t < nb_vertex; ++t) vertices_in_bool.insert({t, false}); + + long nb_vertex_in = 0, nb_tet_in = 0; + + // store the temporary "id" of the tagged tets + std::vector tag_tet_index; + for (size_t i = 0; i < nb_tet; ++i) { + long tri_tag = tag_acc.const_scalar_attribute(tets.at(i)); + std::cout << tri_tag << " "; + switch (tri_tag) { + // inside: store the temp id of this tri + case 1: + nb_tet_in++; + tag_tet_index.push_back(i); + break; + // outside: do nothing + case 0: break; + // neither: runtime error + default: throw std::runtime_error("illegal tag!"); + } + } + assert(nb_tet_in <= nb_tet); + + // for the tagged tri, mark their "real" vertices as inside (duplicates handled by boolean) + // current algo for bug fixing: O(N^2), go over all vertices and look for match, + // only assign tag to inside ones + for (size_t i = 0; i < nb_tet_in; ++i) { + Simplex s = Simplex::tetrahedron(tets[tag_tet_index[i]]); + std::vector tuple_list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + for (wmtk::Tuple t : tuple_list) + vertices_in_bool[find_vertex_index(m, t)] = true; + } + + // construct a map from old tuple to temp new "id" of a "real" vertex + std::map old2new; + for (long i = 0; i < nb_vertex; ++i) { + if (vertices_in_bool[i]) { + // std::cout << "inside! nb_vertex_in = " << nb_vertex_in << std::endl; + // old vertex tuple t mapped to new vertex id j, where j increases by count + old2new.insert({i, nb_vertex_in}); + nb_vertex_in++; + } + } + + wmtk::TetMesh mesh; + wmtk::RowVectors4l tris; + tris.resize(nb_tet_in, 4); + // only put in the extracted ones + for (size_t i = 0; i < nb_tet_in; ++i) { + Simplex s = Simplex::tetrahedron(tets[tag_tet_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(4, -1); + for (int index = 0; index < 4; ++index) + data[index] = old2new[find_vertex_index(m, list[index])]; + tris.row(i) << data[0], data[1], data[2], data[3]; + } + mesh.initialize(tris); // init the topology + + // if told to extract and preserve the coordinates + if (pos) { + Eigen::MatrixXd points_in; + points_in.resize(nb_vertex_in, 3); + wmtk::MeshAttributeHandle pos_handle = + m.get_attribute_handle("position", PrimitiveType::Vertex); + wmtk::ConstAccessor pos_acc = m.create_const_accessor(pos_handle); + for (const Tuple& t : vertices) { + // ignore the outside vertices + long old_index = find_vertex_index(m, t); + if (vertices_in_bool[old_index]) { + points_in.row(old2new[old_index]) = pos_acc.const_vector_attribute(t); + } + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + mesh); + } + } + return mesh; +} + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.hpp b/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.hpp new file mode 100644 index 0000000000..bf15a15b91 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/extract_subset_3d.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { + +wmtk::TetMesh +extract_subset_3d(wmtk::TetMesh m, wmtk::MeshAttributeHandle taghandle, bool pos); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.cpp b/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.cpp new file mode 100644 index 0000000000..10a2e2e648 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.cpp @@ -0,0 +1,139 @@ +#include "topology_separate_2d.hpp" +#include +namespace wmtk::components::internal { + +wmtk::TriMesh topology_separate_2d(wmtk::TriMesh m) +{ + /* + Algorithm idea: + first partition the whole mesh to be m edge-connected components + for each vertex as the connecting vertex of multiple(n) edge-connected components, + for edge-cc i in n, if this vertex is shared by a_i groups of + vertex-connected but not edge-connected triangles, then create a_i copies of this vertex + assign one to each group. So a total of \Sigma_{i = 1}^{n} a_i copies of the same vertex will be + created + + Algorithm steps: (for 2d) + 1. find all the ecc relationships in this trimesh + 2. for each "real" vertex, + init the default nb_cp = 1 to this vertex + if found 2 different triangles vertex-connected at this vertex but not edge-connected, + then assign nb_cp = \Sigma_{i = 1}^{k} a_i , where k is the number of ecc adj to this + vertex, a_i is the number of groups of vertex-connected but not edge-connected triangles in the + i-th ecc + 3. init a counter starting from 0. For vertex starting from 0, for each "real" vertex: + if only 1 copy, assign vertex with a single current counter value, then increment + else assign vertex with {curr, ..., ctr + nb_cp - 1}, then increment by num of copies + 4. use the new ids to reconstruct the mesh + 5. return the new mesh + */ + std::vector faces = m.get_all(wmtk::PrimitiveType::Face); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + int nb_vertex = vertices.size(); + int nb_tri = faces.size(); + std::vector vertex_cp(nb_vertex, 1); // how many copies should we make on this vertex + + // std::cout << "# of tris = " << nb_tri << std::endl; + // std::cout << "# of vertices = " << nb_vertex << std::endl; + + // Prior work: build an face-adjacency list of triangle faces + std::vector> adj_list_faces(nb_tri, std::vector()); + for (long i = 0; i < nb_tri; ++i) { + for (long j = i; j < nb_tri; ++j) { + long edge_con = + edge_connected(m, wmtk::Simplex::face(faces[i]), wmtk::Simplex::face(faces[j])); + if (edge_con != -1) { + if (i != j) { + adj_list_faces[i].push_back(j); + adj_list_faces[j].push_back(i); + } else { + adj_list_faces[i].push_back(j); + } + } + } + } + std::cout << "adj_list_tets = " << std::endl; + for (int i = 0; i < nb_tri; ++i) { + std::cout << i << ": "; + for (auto j : adj_list_faces[i]) std::cout << j << " "; + std::cout << std::endl; + } + + // Step 1: constuct a list of edge-connected components + std::vector> face_cc_list; + std::vector visited_faces(nb_tri, false); + auto condition = [](long face, std::vector&) noexcept { return true; }; + std::vector nullvector = std::vector(nb_tri); + std::iota(nullvector.begin(), nullvector.end(), 1); + for (long i = 0; i < nb_tri; ++i) { + if (visited_faces[i] == false) { + std::vector cc; + dfs(i, visited_faces, cc, adj_list_faces, condition, nullvector); + face_cc_list.push_back(cc); + } + } + long nb_cc = face_cc_list.size(); + // std::cout << "# of cc = " << nb_cc << std::endl; + // print_vv(face_cc_list); + std::vector nb_cc_vec(nb_cc); + for (long i = 0; i < nb_cc; ++i) nb_cc_vec[i] = face_cc_list[i].size(); + + // Step 2: for each vertex on the boundary, count number of group tris around it + // start a version of the algo where we loop over vertices instead of triangles + std::map>> ccav_vector; + for (long i = 0; i < nb_vertex; ++i) { + if (m.is_boundary(vertices[i], PrimitiveType::Vertex)) { + // std::cout << "vertex " << i << " is on boundary. The adjacent faces are: "; + std::vector adj_faces = adj_faces_of_vertex(m, i); + // for (auto j : adj_faces) std::cout << j << " "; + // std::cout << std::endl; + std::vector> ccav = cc_around_vertex(m, adj_faces, adj_list_faces); + vertex_cp[i] = ccav.size(); + ccav_vector[i] = ccav; + } + } + // for (auto j : vertex_cp) std::cout << j << " "; + // std::cout << std::endl; + + // Step 3: assign queues to each vertex + int counter = 0; + std::vector> new_id_of_vertex(nb_vertex, std::vector()); + for (int i = 0; i < nb_vertex; ++i) { + for (int j = 0; j < vertex_cp[i]; ++j) { + new_id_of_vertex[i].push_back(counter); + counter++; + } + } + + // Step 4: reconstruct the mesh + wmtk::TriMesh mesh; + wmtk::RowVectors3l tris; + tris.resize(nb_tri, 3); + for (long i = 0; i < nb_tri; ++i) { + std::vector list = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::face(faces[i]), + PrimitiveType::Vertex); + std::vector data(3, -1); + for (int index = 0; index < 3; ++index) { + long id_v = find_vertex_index(m, list[index]); + if (vertex_cp[id_v] == 1) + data[index] = new_id_of_vertex[id_v][0]; + else { + for (long j = 0; j < ccav_vector[id_v].size(); ++j) { + if (std::find(ccav_vector[id_v][j].begin(), ccav_vector[id_v][j].end(), i) != + ccav_vector[id_v][j].end()) { + data[index] = new_id_of_vertex[id_v][j]; + break; + } + } + } + } + // std::cout << "data = " << data[0] << ", " << data[1] << ", " << data[2] << std::endl; + tris.row(i) << data[0], data[1], data[2]; + } + mesh.initialize(tris); // init the topology + + return mesh; +} +} // namespace wmtk::components::internal diff --git a/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.hpp b/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.hpp new file mode 100644 index 0000000000..22ef7cd5c7 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/topology_separate_2d.hpp @@ -0,0 +1,13 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { + +wmtk::TriMesh topology_separate_2d(wmtk::TriMesh m); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.cpp b/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.cpp new file mode 100644 index 0000000000..3368f8506f --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.cpp @@ -0,0 +1,234 @@ +#include "topology_separate_3d.hpp" + +namespace wmtk::components::internal { + +wmtk::TetMesh topology_separate_3d_old(wmtk::TetMesh m) +{ + std::vector edges = m.get_all(wmtk::PrimitiveType::Edge); + std::vector tets = m.get_all(wmtk::PrimitiveType::Tetrahedron); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + long nb_vertex = vertices.size(); + long nb_tet = tets.size(); + long nb_edge = edges.size(); + + std::vector edge_cp(nb_edge, 1); // how many copies should we make on this edge + std::vector vertex_cp(nb_vertex, 0); // how many copies should we make on this vertex + + // std::cout << "\n# of tets = " << nb_tet << std::endl; + // std::cout << "# of edges = " << nb_edge << std::endl; + // std::cout << "# of vertices = " << nb_vertex << std::endl; + + // for (int i = 0; i < nb_tet; ++i) { + // std::cout << "tet " << i << ": \n"; + // auto ee = wmtk::simplex::faces_single_dimension( + // m, + // wmtk::Simplex::tetrahedron(tets[i]), + // PrimitiveType::Edge); + // for (int j = 0; j < 6; ++j) { + // auto vers = wmtk::simplex::faces_single_dimension(m, wmtk::Simplex::edge(ee[j]), + // PrimitiveType::Vertex); std::cout << "edge " << find_vertex_index(m, vers[0]) << + // find_vertex_index(m, vers[1]) << "with index " << find_edge_index(m, ee[j]) << + // std::endl; + // } + // } + + // Prior work: build an face-adjacency list of tets + std::vector> adj_list_tets(nb_tet, std::vector()); + for (long i = 0; i < nb_tet; ++i) { + for (long j = i; j < nb_tet; ++j) { + long face_con = face_connected( + m, + wmtk::Simplex::tetrahedron(tets[i]), + wmtk::Simplex::tetrahedron(tets[j])); + if (face_con != -1) { + if (i != j) { + adj_list_tets[i].push_back(j); + adj_list_tets[j].push_back(i); + } else { + adj_list_tets[i].push_back(j); + } + } + } + } + // std::cout << "adj_list_tets = " << std::endl; + // for (int i = 0; i < nb_tet; ++i) { + // std::cout << i << ": "; + // for (auto j : adj_list_tets[i]) std::cout << j << " "; + // std::cout << std::endl; + // } + + // Step 1: constuct a list of face-connected components + std::vector> face_cc_list; + std::vector visited_faces(nb_tet, false); + auto condition = [](long face, std::vector&) noexcept { return true; }; + std::vector nullvector = std::vector(nb_tet); + std::iota(nullvector.begin(), nullvector.end(), 1); + for (long i = 0; i < nb_tet; ++i) { + if (visited_faces[i] == false) { + std::vector cc; + dfs(i, visited_faces, cc, adj_list_tets, condition, nullvector); + face_cc_list.push_back(cc); + } + } + long nb_cc = face_cc_list.size(); + // std::cout << "# of cc = " << nb_cc << std::endl; + std::vector nb_cc_vec(nb_cc); + for (long i = 0; i < nb_cc; ++i) nb_cc_vec[i] = face_cc_list[i].size(); + + // Step 1.5: construct the inverse, tet -> cc + std::vector tet_cc(nb_tet); + for (long i = 0; i < face_cc_list.size(); ++i) { + for (long j : face_cc_list[i]) tet_cc[j] = i; + } + + // Step 2: for each edge on the boundary, count number of group tets around it + // std::map>> ccav_vector; + std::map> ccav_vector; + for (long i = 0; i < nb_edge; ++i) { + if (m.is_boundary(edges[i], PrimitiveType::Edge)) { + // std::cout << "edge " << i << " is on boundary. "; + // auto vers = wmtk::simplex::faces_single_dimension(m, wmtk::Simplex::edge(edges[i]), + // PrimitiveType::Vertex); std::cout << " edge " << find_vertex_index(m, vers[0]) << + // find_vertex_index(m, vers[1]) << std::endl; + // std::cout << "The adjacent faces are: "; + std::vector adj_tets = adj_tets_of_edge(m, i); + // for (auto j : adj_tets) std::cout << j << " "; + // std::cout << std::endl; + std::vector> ccav = tet_cc_around_tuple(m, adj_tets, adj_list_tets); + edge_cp[i] = ccav.size(); + std::vector cc; + for (std::vector j : ccav) cc.push_back(tet_cc[j[0]]); + ccav_vector[i] = cc; + // std::cout << "ccav.size() = " << ccav.size() << std::endl; + // ccav_vector[i] = ccav; + } + } + // for (long j : edge_cp) std::cout << j << " "; + // std::cout << std::endl; + + // Step 2.5(with question): load the number of copies of edges to vertices + // CAREFUL: copies of a vertex should be the number of face-connected components in tets around + // all non-manifold edges incident at this vertex + std::vector> vertex_cc_set(nb_vertex, std::set()); + for (int i = 0; i < nb_edge; ++i) { + if (edge_cp[i] != 1) { + std::vector edge_vertices = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::edge(edges[i]), + PrimitiveType::Vertex); + vertex_cc_set[find_vertex_index(m, edge_vertices[0])].insert( + ccav_vector[i].begin(), + ccav_vector[i].end()); + vertex_cc_set[find_vertex_index(m, edge_vertices[1])].insert( + ccav_vector[i].begin(), + ccav_vector[i].end()); + // vertex_cp[find_vertex_index(m, edge_vertices[0])] += edge_cp[i]; + // vertex_cp[find_vertex_index(m, edge_vertices[1])] += edge_cp[i]; + } + } + for (int i = 0; i < nb_vertex; ++i) vertex_cp[i] += vertex_cc_set[i].size(); + // for (long j : vertex_cp) std::cout << j << " "; + // std::cout << std::endl; + + // Step 2.9: deal with vertices on the boundary, count number of group tets around it. Make sure + // the already edge-connected ones are not in count + std::map>> ccav_vector_vertex; + for (long i = 0; i < nb_vertex; ++i) { + if (m.is_boundary(vertices[i], PrimitiveType::Vertex)) { + // std::cout << "vertex " << i << " is on boundary. The adjacent faces are: "; + std::vector adj_tets = adj_tets_of_vertex(m, i); + // for (long j : adj_tets) std::cout << j << " "; + // std::cout << std::endl; + std::vector> ccav = tet_cc_around_tuple(m, adj_tets, adj_list_tets); + // std::cout << "ccav.size() = " << ccav.size() << std::endl; + // for more than 1 cc, we need to check if 2 cc are connected by edges, which we have + // already counted + if (ccav.size() > 1) { + std::vector cc_flag(ccav.size(), true); + // std::cout << "before processing, cc_flag = "; + // for (bool f : cc_flag) std::cout << f << " "; + // std::cout << std::endl; + for (int j = 0; j < ccav.size(); ++j) { + for (int k = j + 1; k < ccav.size(); ++k) { + if (cc_flag[j] == false && cc_flag[k] == false) continue; + if (edge_connected( + m, + wmtk::Simplex::tetrahedron(tets[ccav[j][0]]), + wmtk::Simplex::tetrahedron(tets[ccav[k][0]]))) { + cc_flag[j] = false; + cc_flag[k] = false; + break; + } + } + } + // std::cout << "after processing, cc_flag = "; + // for (bool f : cc_flag) std::cout << f << " "; + // std::cout << std::endl; + long cc_to_count = 0; + for (int j = 0; j < cc_flag.size(); ++j) { + if (cc_flag[j] == true) cc_to_count++; + } + vertex_cp[i] += cc_to_count; + } + ccav_vector_vertex[i] = ccav; + } + // else { + // std::cout << "vertex " << i << " is not on boundary. The adjacent faces are: "; + // std::vector adj_tets = adj_tets_of_vertex(m, i); + // for (auto j : adj_tets) std::cout << j << " "; + // std::cout << std::endl; + // } + } + // for (long j : vertex_cp) std::cout << j << " "; + // std::cout << std::endl; + + // Step 3: assign new id to each vertex + int counter = 0; + std::vector> new_id_of_vertex(nb_vertex, std::vector()); + for (int i = 0; i < nb_vertex; ++i) { + if (vertex_cp[i] == 0) { + new_id_of_vertex[i].push_back(counter); + counter++; + continue; + } + for (int j = 0; j < vertex_cp[i]; ++j) { + new_id_of_vertex[i].push_back(counter); + counter++; + } + } + + // Step 4: reconstruct the mesh + wmtk::TetMesh mesh; + wmtk::RowVectors4l tris; + tris.resize(nb_tet, 4); + for (long i = 0; i < nb_tet; ++i) { + std::vector list = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::tetrahedron(tets[i]), + PrimitiveType::Vertex); + std::vector data(4, -1); + for (int index = 0; index < 4; ++index) { + long id_v = find_vertex_index(m, list[index]); + if (vertex_cp[id_v] == 1) + data[index] = new_id_of_vertex[id_v][0]; + else { + for (long j = 0; j < ccav_vector_vertex[id_v].size(); ++j) { + if (std::find( + ccav_vector_vertex[id_v][j].begin(), + ccav_vector_vertex[id_v][j].end(), + i) != ccav_vector_vertex[id_v][j].end()) { + data[index] = new_id_of_vertex[id_v][j]; + break; + } + } + } + } + // std::cout << "data = " << data[0] << ", " << data[1] << ", " << data[2] << ", " << data[3] + // << std::endl; + tris.row(i) << data[0], data[1], data[2], data[3]; + } + mesh.initialize(tris); // init the topology + return mesh; + // return m; +} +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.hpp b/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.hpp new file mode 100644 index 0000000000..ca06b59e55 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/attic/topology_separate_3d.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { +wmtk::TetMesh topology_separate_3d_old(wmtk::TetMesh m); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/generate_submesh.cpp b/components/wmtk_components/extract_subset/internal/generate_submesh.cpp new file mode 100644 index 0000000000..a1c0949844 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/generate_submesh.cpp @@ -0,0 +1,281 @@ +#include "generate_submesh.hpp" + +namespace wmtk::components::internal { + +// Getting submesh and operate on it is essential because extracting subset will change connectivity +std::unique_ptr +generate_submesh(wmtk::Mesh& m, wmtk::MeshAttributeHandle tag_handle, bool pos) +{ + /* + new Algo for this first part of plainly getting a subset from the original mesh: + 1. register a new attribute to each vertex and init as -1 + 2. store all the indices of tagged top simplices + 3. for each tagged top dim simplex, for each vertex of the simplex, + if the vertex has attribute -1, then assign it a new index, increment the counter (leave it + there if already having an index) + + Reconstruction: For each top dim simplex, get the tag if tag is 1 + get attribute of all vertices of the simplex, and store them in a new mesh + if pos is true, get the position of the vertices and store them in the new mesh + */ + + int top_simplex_dim = m.top_cell_dimension(); + if (top_simplex_dim != 2 && top_simplex_dim != 3) + throw std::runtime_error("Invalid top dimension in separating topology!"); + wmtk::Accessor tag_acc = m.create_accessor(tag_handle); + wmtk::PrimitiveType topType = m.top_simplex_type(); + std::vector top_simplices = m.get_all(topType); + long top_simplex_count = top_simplices.size(); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + int nb_vertex = vertices.size(); + std::map vertices_in_bool; + for (long t = 0; t < nb_vertex; ++t) vertices_in_bool.insert({t, false}); + /* + long nb_vertex_in = 0, nb_cell_in = 0; + std::vector tag_simplex_index; + for (size_t i = 0; i < top_simplex_count; ++i) { + long cell_tag = tag_acc.const_scalar_attribute(top_simplices.at(i)); + switch (cell_tag) { + // inside: store the temp id of this cell + case 1: tag_simplex_index.push_back(i); break; + // outside: do nothing + case 0: break; + // neither: runtime error + default: throw std::runtime_error("illegal tag!"); + } + } + nb_cell_in = tag_simplex_index.size(); + assert(nb_cell_in <= top_simplex_count); + std::cout << "# of cell inside = " << nb_cell_in << std::endl; + + // TODO: improve the algorithm to achieve O(N) + for (size_t i = 0; i < nb_cell_in; ++i) { + Simplex s = (top_simplex_dim == 2) + ? Simplex::face(top_simplices[tag_simplex_index[i]]) + : Simplex::tetrahedron(top_simplices[tag_simplex_index[i]]); + std::vector tuple_list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + for (wmtk::Tuple t : tuple_list) vertices_in_bool[find_vertex_index(m, t)] = true; + } + + std::map old2new; + for (long i = 0; i < nb_vertex; ++i) { + if (vertices_in_bool[i]) { + // std::cout << "inside! vertex_index = " << i << std::endl; + // old vertex tuple t mapped to new vertex id j, where j increases by count + old2new.insert({i, nb_vertex_in}); + nb_vertex_in++; + } + } + std::cout << "nb_vertex_in = " << nb_vertex_in << std::endl; + + wmtk::TriMesh tri_ext_mesh; + wmtk::RowVectors3l tri_exts; + wmtk::TetMesh tet_ext_mesh; + wmtk::RowVectors4l tet_exts; + + if (top_simplex_dim == 2) { + tri_exts.resize(nb_cell_in, 3); + for (int i = 0; i < nb_cell_in; ++i) { + Simplex s = Simplex::face(top_simplices[tag_simplex_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(3, -1); + for (int index = 0; index < 3; ++index) + data[index] = old2new[find_vertex_index(m, list[index])]; + tri_exts.row(i) << data[0], data[1], data[2]; + } + // std::cout << "tri_exts = " << tri_exts << std::endl; + tri_ext_mesh.initialize(tri_exts); + assert(tri_ext_mesh.get_all(topType).size() == nb_cell_in); + if (!pos) return std::make_unique(tri_ext_mesh); + } else if (top_simplex_dim == 3) { + tet_exts.resize(nb_cell_in, m.top_cell_dimension() + 1); + for (size_t i = 0; i < nb_cell_in; ++i) { + Simplex s = Simplex::tetrahedron(top_simplices[tag_simplex_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(m.top_cell_dimension() + 1, -1); + for (int index = 0; index < 4; ++index) + data[index] = old2new[find_vertex_index(m, list[index])]; + tet_exts.row(i) << data[0], data[1], data[2], data[3]; + } + tet_ext_mesh.initialize(tet_exts); + assert(tet_ext_mesh.get_all(topType).size() == nb_cell_in); + if (!pos) return std::make_unique(tet_ext_mesh); + } + + // if told to, extract and preserve the coordinates + if (pos) { + Eigen::MatrixXd points_in; + points_in.resize(nb_vertex_in, m.top_cell_dimension()); + wmtk::MeshAttributeHandle pos_handle = + m.get_attribute_handle("position", PrimitiveType::Vertex); + wmtk::ConstAccessor pos_acc = m.create_const_accessor(pos_handle); + for (const Tuple& t : vertices) { + // ignore the outside vertices + long old_index = find_vertex_index(m, t); + if (vertices_in_bool[old_index]) { + points_in.row(old2new[old_index]) = pos_acc.const_vector_attribute(t); + } + } + // call the set_matrix_attribute function according to the top dimension + switch (m.top_cell_dimension()) { + case 2: + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + tri_ext_mesh); + top_simplices = tri_ext_mesh.get_all(topType); + return std::make_unique(tri_ext_mesh); + case 3: + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + tet_ext_mesh); + top_simplices = tet_ext_mesh.get_all(topType); + return std::make_unique(tet_ext_mesh); + default: throw std::runtime_error("Invalid top dimension in separating topology!"); + } + } + throw std::runtime_error("Should not reach here!"); + + */ + // The following is an implementation of the algo listed in the comment above + ///* + // Step 1: register a new attribute to each vertex and init as -1 + wmtk::VectorXl submesh_index_vector; + submesh_index_vector.resize(nb_vertex, 1); + for (long i = 0; i < nb_vertex; ++i) submesh_index_vector.row(i) << -1; + wmtk::MeshAttributeHandle new_ver_index_handle = wmtk::mesh_utils::set_matrix_attribute( + submesh_index_vector, + "submesh_index", + wmtk::PrimitiveType::Vertex, + m); + wmtk::Accessor new_ver_index_acc = m.create_accessor(new_ver_index_handle); + + long nb_vertex_in = 0, nb_cell_in = 0; + // Step 2: store all the indices of tagged top simplices + std::vector tag_simplex_index; + for (size_t i = 0; i < top_simplex_count; ++i) { + long tri_tag = tag_acc.const_scalar_attribute(top_simplices.at(i)); + switch (tri_tag) { + // inside: store the temp id of this tri + case 1: tag_simplex_index.push_back(i); break; + // outside: do nothing + case 0: break; + // neither: runtime error + default: throw std::runtime_error("illegal tag!"); + } + } + nb_cell_in = tag_simplex_index.size(); + // std::cout << "# of cell inside = " << nb_cell_in << std::endl; + // std::cout << "index of cells inside = " << tag_simplex_index << std::endl; + assert(nb_cell_in <= top_simplex_count); + + // Step 3.1: for each tagged top dim simplex + for (size_t index : tag_simplex_index) { + // std::cout << "index = " << index; + // std::cout << "tag = " << tag_acc.const_scalar_attribute(top_simplices.at(index)) + // << std::endl; + Simplex s = (top_simplex_dim == 2) ? Simplex::face(top_simplices[index]) + : Simplex::tetrahedron(top_simplices[index]); + std::vector vertices_in_simplex = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + // for each vertex of the simplex + for (wmtk::Tuple t : vertices_in_simplex) { + // std::cout << "current vertex id = " << find_vertex_index(m, t) << std::endl; + // std::cout << "before: attribute = " << new_ver_index_acc.scalar_attribute(t) + // << std::endl; + // Step 3.2: if the vertex has attribute -1, then assign it a new index, + if (new_ver_index_acc.scalar_attribute(t) == -1) { + new_ver_index_acc.scalar_attribute(t) = nb_vertex_in; + nb_vertex_in++; + } + // std::cout << "after: attribute = " << new_ver_index_acc.scalar_attribute(t) + // << std::endl; + } + } + // std::cout << "nb_vertex_in = " << nb_vertex_in << std::endl; + + wmtk::TriMesh tri_ext_mesh; + wmtk::RowVectors3l tri_exts; + wmtk::TetMesh tet_ext_mesh; + wmtk::RowVectors4l tet_exts; + if (top_simplex_dim == 2) { + tri_exts.resize(nb_cell_in, 3); + for (int i = 0; i < nb_cell_in; ++i) { + if (tag_acc.const_scalar_attribute(top_simplices.at(tag_simplex_index[i])) == 0) + continue; + Simplex s = Simplex::face(top_simplices[tag_simplex_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(3, -1); + for (int index = 0; index < 3; ++index) { + data[index] = new_ver_index_acc.scalar_attribute(list[index]); + // std::cout << "data[" << index << "] = " << data[index] << std::endl; + } + tri_exts.row(i) << data[0], data[1], data[2]; + } + // std::cout << "tri_exts = " << tri_exts << std::endl; + tri_ext_mesh.initialize(tri_exts); + assert(tri_ext_mesh.get_all(topType).size() == nb_cell_in); + if (!pos) return std::make_unique(tri_ext_mesh); + } else if (top_simplex_dim == 3) { + tet_exts.resize(nb_cell_in, 4); + for (size_t i = 0; i < nb_cell_in; ++i) { + if (tag_acc.const_scalar_attribute(top_simplices.at(tag_simplex_index[i])) == 0) + continue; + Simplex s = Simplex::tetrahedron(top_simplices[tag_simplex_index[i]]); + std::vector list = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + std::vector data(4, -1); + for (int index = 0; index < 4; ++index) + data[index] = new_ver_index_acc.scalar_attribute(list[index]); + tet_exts.row(i) << data[0], data[1], data[2], data[3]; + } + tet_ext_mesh.initialize(tet_exts); + assert(tet_ext_mesh.get_all(topType).size() == nb_cell_in); + if (!pos) return std::make_unique(tet_ext_mesh); + } + + // if told to, extract and preserve the coordinates + if (pos) { + Eigen::MatrixXd points_in; + points_in.resize(nb_vertex_in, m.top_cell_dimension()); + wmtk::MeshAttributeHandle pos_handle = + m.get_attribute_handle("position", PrimitiveType::Vertex); + wmtk::ConstAccessor pos_acc = m.create_const_accessor(pos_handle); + for (int i = 0; i < nb_vertex; ++i) { + if (new_ver_index_acc.scalar_attribute(vertices[i]) != -1) { + points_in.row(new_ver_index_acc.scalar_attribute(vertices[i])) = + pos_acc.const_vector_attribute(vertices[i]); + } + } + // call the set_matrix_attribute function according to the top dimension + switch (m.top_cell_dimension()) { + case 2: + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + tri_ext_mesh); + top_simplices = tri_ext_mesh.get_all(topType); + return std::make_unique(tri_ext_mesh); + case 3: + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + tet_ext_mesh); + top_simplices = tet_ext_mesh.get_all(topType); + return std::make_unique(tet_ext_mesh); + default: throw std::runtime_error("Invalid top dimension in separating topology!"); + } + } + throw std::runtime_error("Should not reach here!"); + //*/ +} +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/generate_submesh.hpp b/components/wmtk_components/extract_subset/internal/generate_submesh.hpp new file mode 100644 index 0000000000..5bd1d188b3 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/generate_submesh.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { + +std::unique_ptr +generate_submesh(wmtk::Mesh& m, wmtk::MeshAttributeHandle taghandle, bool pos); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/new_topology_separate.cpp b/components/wmtk_components/extract_subset/internal/new_topology_separate.cpp new file mode 100644 index 0000000000..ceb2762db0 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/new_topology_separate.cpp @@ -0,0 +1,158 @@ +#include "new_topology_separate.hpp" +#include +namespace wmtk::components::internal { + +// general function to separate topology, regardless of dimension +std::unique_ptr topology_separate(wmtk::Mesh& m, bool pos) +{ + // First extract the non-manifold submesh, then execute the following algo + /*Algo outline: + 1. get subset of top simplices with tag 1, reconstruct a submesh that is non-manifold + 2. for each top simplex, get all corners and assign them a unique index + 3. for each corner, get all top simplices sharing the corner and update their duplicate index + 4. create a new mesh and copy the topology + 5. if pos is true, copy the geometry as well + */ + int top_simplex_dim = m.top_cell_dimension(); + if (top_simplex_dim != 2 && top_simplex_dim != 3) + throw std::runtime_error("Invalid top dimension in separating topology!"); + wmtk::PrimitiveType topType = m.top_simplex_type(); + std::vector top_simplices = m.get_all(topType); + long top_simplex_count = top_simplices.size(); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + int nb_vertex = vertices.size(); + + // Now begins the second part of the algorithm: make it topologically manifold + long counter = 0; + + // first, register a vector attribute to store the corner ids for each top dimension simplex + wmtk::RowVectors4l dup; + dup.resize(top_simplex_count, top_simplex_dim + 1); + for (long i = 0; i < top_simplex_count; ++i) { + switch (top_simplex_dim) { + case 2: dup.row(i) << -1, -1, -1; break; + case 3: dup.row(i) << -1, -1, -1, -1; break; + default: break; + } + } + wmtk::MeshAttributeHandle duplicate_handle = + wmtk::mesh_utils::set_matrix_attribute(dup, "duplicate_index", topType, m); + wmtk::Accessor dup_acc = m.create_accessor(duplicate_handle); + + // second, go over all top dimension simplices and adjust the duplicate index + for (long i = 0; i < top_simplex_count; ++i) { + auto v = dup_acc.vector_attribute(top_simplices[i]); + // Question: Why can't I use the constructor here to build a Simplex? + Simplex s = (top_simplex_dim == 2) ? Simplex::face(top_simplices[i]) + : Simplex::tetrahedron(top_simplices[i]); + // get all corners for current top simplex + std::vector corners = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + // std::cout << "Hello1, # of corners = " << corners.size() << std::endl; + for (long j = 0; j < corners.size(); ++j) { + // check whether it has been visited + // std::cout << "local_vid = " << corners[j].get_local_vid() + // << ", gid = " << find_vertex_index(m, corners[j]) + // << ", v[local_vid] = " << v[corners[j].get_local_vid()] << std::endl; + // BE CAREFUL: the vertex id is not the same as the index of the corner in the simplex + if (v[corners[j].get_local_vid()] != -1) continue; + // if the corner has not been assigned a duplicate index, assign it + v[corners[j].get_local_vid()] = counter; + // std::cout << "verify: v[j] = " + // << dup_acc.vector_attribute(top_simplices[i])[corners[j].get_local_vid()] + // << std::endl; + + // find all top simplices sharing the same corner vertex, + // update duplicate index of theie corner accordingly + + // get all top dimension simplices sharing the vertex and are face-connected + wmtk::simplex::SimplexCollection sc = + wmtk::simplex::top_dimension_cofaces(m, Simplex::vertex(corners[j])); + // std::cout << "# of adj corners simplices = " << sc.simplex_vector().size() << + // std::endl; + for (wmtk::Simplex adj_simplex : sc) { + // tuple for a top dimension simplex would be the same as tuple for the corner + wmtk::Tuple adj_corner_tuple = adj_simplex.tuple(); + // std::cout << "vertex id = " << find_vertex_index(m, adj_corner_tuple) << + // std::endl; + auto adj_vector = dup_acc.vector_attribute(adj_corner_tuple); + long k = adj_corner_tuple.get_local_vid(); + // std::cout << "before adjusting = " << adj_vector[k] << std::endl; + if (adj_vector[k] == counter) continue; + if (adj_vector[k] != -1 && adj_vector[k] != counter) + throw std::runtime_error("Duplicate index conflict!"); + adj_vector[k] = counter; + // std::cout << "after adjusting = " << + // dup_acc.vector_attribute(adj_corner_tuple)[k] + // << std::endl; + } + // finally, increment the counter + counter++; + } + } + + // third, create a new mesh and copy the topology + Eigen::MatrixXd points_in; + if (pos) { + std::map vertices_pos_visited; + for (long t = 0; t < counter + 1; ++t) vertices_pos_visited.insert({t, false}); + points_in.resize(counter + 1, top_simplex_dim); + wmtk::MeshAttributeHandle pos_handle = + m.get_attribute_handle("position", PrimitiveType::Vertex); + wmtk::ConstAccessor pos_acc = m.create_const_accessor(pos_handle); + for (int i = 0; i < top_simplex_count; ++i) { + Simplex s = (top_simplex_dim == 2) ? Simplex::face(top_simplices[i]) + : Simplex::tetrahedron(top_simplices[i]); + std::vector corners = + wmtk::simplex::faces_single_dimension(m, s, PrimitiveType::Vertex); + for (int j = 0; j < top_simplex_dim + 1; ++j) { + long index = dup_acc.vector_attribute(top_simplices[i])[j]; + if (!vertices_pos_visited[index]) { + points_in.row(index) = pos_acc.const_vector_attribute(corners[j]); + vertices_pos_visited[index] = true; + } + } + } + } + + switch (m.top_cell_dimension()) { + case 2: { + wmtk::TriMesh mesh; + wmtk::RowVectors3l tris; + tris.resize(top_simplex_count, 3); + for (long i = 0; i < top_simplex_count; ++i) { + auto v = dup_acc.vector_attribute(top_simplices[i]); + tris.row(i) << v[0], v[1], v[2]; + } + mesh.initialize(tris); + if (pos) { + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + mesh); + } + return std::make_unique(mesh); + } + case 3: { + wmtk::TetMesh mesh; + wmtk::RowVectors4l tets; + tets.resize(top_simplex_count, 4); + for (long i = 0; i < top_simplex_count; ++i) { + auto v = dup_acc.vector_attribute(top_simplices[i]); + tets.row(i) << v[0], v[1], v[2], v[3]; + } + mesh.initialize(tets); + if (pos) { + wmtk::mesh_utils::set_matrix_attribute( + points_in, + "position", + wmtk::PrimitiveType::Vertex, + mesh); + } + return std::make_unique(mesh); + } + default: throw std::runtime_error("Invalid mesh dimension in separating topology!"); + } +} +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/new_topology_separate.hpp b/components/wmtk_components/extract_subset/internal/new_topology_separate.hpp new file mode 100644 index 0000000000..d18f42e741 --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/new_topology_separate.hpp @@ -0,0 +1,15 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include "utils.hpp" + +namespace wmtk::components::internal { + +std::unique_ptr topology_separate(wmtk::Mesh& m, bool pos); +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/utils.cpp b/components/wmtk_components/extract_subset/internal/utils.cpp new file mode 100644 index 0000000000..410983f2da --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/utils.cpp @@ -0,0 +1,225 @@ +#include "utils.hpp" + +namespace wmtk::components::internal { +// template +long connected( + const wmtk::Mesh& m, + wmtk::Simplex i, + wmtk::Simplex j, + std::function extractor, + wmtk::PrimitiveType type) +{ + std::vector primitives = m.get_all(type); + std::vector i_tuple_list = wmtk::simplex::faces_single_dimension(m, i, type); + std::vector j_tuple_list = wmtk::simplex::faces_single_dimension(m, j, type); + for (int a = 0; a < i_tuple_list.size(); ++a) { + for (int b = 0; b < j_tuple_list.size(); ++b) { + if (m.simplices_are_equal(extractor(i_tuple_list[a]), extractor(j_tuple_list[b]))) { + return find_index(m, i_tuple_list[a], extractor, type); + } + } + } + return -1; +} + +long face_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j) +{ + return connected( + m, + i, + j, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::face(tuple); }, + wmtk::PrimitiveType::Face); +} + +long edge_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j) +{ + return connected( + m, + i, + j, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::edge(tuple); }, + wmtk::PrimitiveType::Edge); +} + +long vertex_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j) +{ + return connected( + m, + i, + j, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::vertex(tuple); }, + wmtk::PrimitiveType::Vertex); +} + +long find_index( + const wmtk::Mesh& m, + wmtk::Tuple t, + std::function extractFunction, + wmtk::PrimitiveType type) +{ + std::vector primitives = m.get_all(type); + for (int i = 0; i < primitives.size(); ++i) { + if (m.simplices_are_equal(extractFunction(primitives[i]), extractFunction(t))) { + return i; + } + } + return -1; +} + +long find_edge_index(const wmtk::Mesh& m, wmtk::Tuple t) +{ + return find_index( + m, + t, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::edge(tuple); }, + wmtk::PrimitiveType::Edge); +} + +long find_vertex_index(const wmtk::Mesh& m, wmtk::Tuple t) +{ + return find_index( + m, + t, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::vertex(tuple); }, + wmtk::PrimitiveType::Vertex); +} + +long find_face_index(const wmtk::Mesh& m, wmtk::Tuple t) +{ + return find_index( + m, + t, + [](const wmtk::Tuple& tuple) { return wmtk::Simplex::face(tuple); }, + wmtk::PrimitiveType::Face); +} + +std::vector adj_faces_of_vertex(const wmtk::TriMesh& m, long i) +{ + // Algo: given a vertex, traverse all faces in the mesh, for each face, find all the vertices + // if the vertex we are checking is in the list, then add the face index to a list to return + std::vector faces = m.get_all(wmtk::PrimitiveType::Face); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + std::vector adj_faces; + for (wmtk::Tuple face : faces) { + std::vector face_vertices = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::face(face), + PrimitiveType::Vertex); + for (wmtk::Tuple vertex : face_vertices) { + if (m.simplices_are_equal( + wmtk::Simplex::vertex(vertex), + wmtk::Simplex::vertex(vertices[i]))) { + adj_faces.push_back(find_face_index(m, face)); + break; + } + } + } + return adj_faces; +} + +std::vector adj_tets_of_edge(const wmtk::TetMesh& m, long i) +{ + std::vector tets = m.get_all(wmtk::PrimitiveType::Tetrahedron); + std::vector edges = m.get_all(wmtk::PrimitiveType::Edge); + std::vector adj_tets; + for (long j = 0; j < tets.size(); ++j) { + wmtk::Tuple tet = tets[j]; + std::vector tet_edges = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::tetrahedron(tet), + PrimitiveType::Edge); + for (wmtk::Tuple edge : tet_edges) { + if (m.simplices_are_equal(wmtk::Simplex::edge(edge), wmtk::Simplex::edge(edges[i]))) { + adj_tets.push_back(j); + break; + } + } + } + return adj_tets; +} + +std::vector adj_tets_of_vertex(const wmtk::TetMesh& m, long i) +{ + std::vector tets = m.get_all(wmtk::PrimitiveType::Tetrahedron); + std::vector vertices = m.get_all(wmtk::PrimitiveType::Vertex); + std::vector adj_tets; + for (long j = 0; j < tets.size(); ++j) { + wmtk::Tuple tet = tets[j]; + std::vector tet_vertices = wmtk::simplex::faces_single_dimension( + m, + wmtk::Simplex::tetrahedron(tet), + PrimitiveType::Vertex); + for (wmtk::Tuple vertex : tet_vertices) { + if (m.simplices_are_equal( + wmtk::Simplex::vertex(vertex), + wmtk::Simplex::vertex(vertices[i]))) { + adj_tets.push_back(j); + break; + } + } + } + return adj_tets; +} + +void dfs( + long start, + std::vector& visited, + std::vector& cc, + const std::vector>& adj, + const std::function&)>& condition, + std::vector& candidates) +{ + visited[start] = true; + cc.push_back(start); + for (long j : adj[start]) { + if (!visited[j] && condition(j, candidates)) { + dfs(j, visited, cc, adj, condition, candidates); + } + } +} + +std::vector> cc_around_vertex( + const wmtk::TriMesh& m, + std::vector& adj_faces, + std::vector>& adj_list_faces) +{ + // use exactly the same also as finding connected components in the whole mesh to find cc here + std::vector> face_cc_list; + std::vector visited_faces(m.get_all(wmtk::PrimitiveType::Face).size(), false); + auto condition = [](long face, std::vector& candidates) { + return std::find(candidates.begin(), candidates.end(), face) != candidates.end(); + }; + for (long i = 0; i < adj_faces.size(); ++i) { + // std::cout << "adj_faces[i] = " << adj_faces[i] << std::endl; + if (visited_faces[adj_faces[i]] == false) { + std::vector cc; + dfs(adj_faces[i], visited_faces, cc, adj_list_faces, condition, adj_faces); + face_cc_list.push_back(cc); + } + } + return face_cc_list; +} + +std::vector> tet_cc_around_tuple( + const wmtk::TetMesh& m, + std::vector& adj_tets, + std::vector>& adj_list_tets) +{ + std::vector> tet_cc_list; + std::vector visited_tets(m.get_all(wmtk::PrimitiveType::Tetrahedron).size(), false); + auto condition = [](long face, std::vector& candidates) { + return std::find(candidates.begin(), candidates.end(), face) != candidates.end(); + }; + for (long i = 0; i < adj_tets.size(); ++i) { + // std::cout << "adj_faces[i] = " << adj_faces[i] << std::endl; + if (visited_tets[adj_tets[i]] == false) { + std::vector cc; + dfs(adj_tets[i], visited_tets, cc, adj_list_tets, condition, adj_tets); + tet_cc_list.push_back(cc); + } + } + return tet_cc_list; +} + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wmtk_components/extract_subset/internal/utils.hpp b/components/wmtk_components/extract_subset/internal/utils.hpp new file mode 100644 index 0000000000..f4647d8e6e --- /dev/null +++ b/components/wmtk_components/extract_subset/internal/utils.hpp @@ -0,0 +1,61 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace wmtk::components::internal { + +// template +long connected( + const wmtk::Mesh& m, + wmtk::Simplex i, + wmtk::Simplex j, + std::function extractor, + wmtk::PrimitiveType type); + +long face_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j); + +long edge_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j); + +long vertex_connected(const wmtk::Mesh& m, wmtk::Simplex i, wmtk::Simplex j); + +long find_index( + const wmtk::Mesh& m, + wmtk::Tuple t, + std::function extractFunction, + wmtk::PrimitiveType type); + +long find_edge_index(const wmtk::Mesh& m, wmtk::Tuple t); + +long find_vertex_index(const wmtk::Mesh& m, wmtk::Tuple t); + +long find_face_index(const wmtk::Mesh& m, wmtk::Tuple t); + +std::vector adj_faces_of_vertex(const wmtk::TriMesh& m, long i); + +std::vector adj_tets_of_edge(const wmtk::TetMesh& m, long i); + +std::vector adj_tets_of_vertex(const wmtk::TetMesh& m, long i); + +void dfs( + long start, + std::vector& visited, + std::vector& cc, + const std::vector>& adj, + const std::function&)>& condition, + std::vector& candidates); + +std::vector> cc_around_vertex( + const wmtk::TriMesh& m, + std::vector& adj_faces, + std::vector>& adj_list_faces); + +std::vector> tet_cc_around_tuple( + const wmtk::TetMesh& m, + std::vector& adj_tets, + std::vector>& adj_list_tets); + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/src/wmtk/TetMesh.cpp b/src/wmtk/TetMesh.cpp index 65546b142e..1c92780eef 100644 --- a/src/wmtk/TetMesh.cpp +++ b/src/wmtk/TetMesh.cpp @@ -388,10 +388,15 @@ bool TetMesh::is_boundary_face(const Tuple& tuple) const return tt_accessor.vector_attribute(tuple)(tuple.m_local_fid) < 0; } -bool TetMesh::is_boundary_edge(const Tuple& vertex) const +bool TetMesh::is_boundary_edge(const Tuple& edge) const { - assert(false); - throw std::runtime_error("NotImplemented"); + const SimplicialComplex neigh = SimplicialComplex::open_star(*this, Simplex::edge(edge)); + for (const Simplex& s : neigh.get_faces()) { + if (is_boundary(s.tuple())) { + return true; + } + } + return false; } bool TetMesh::is_boundary_vertex(const Tuple& vertex) const { diff --git a/src/wmtk/Tuple.cpp b/src/wmtk/Tuple.cpp index 7df9066071..d13c89772e 100644 --- a/src/wmtk/Tuple.cpp +++ b/src/wmtk/Tuple.cpp @@ -59,4 +59,11 @@ Tuple Tuple::with_updated_hash(long new_hash) const return Tuple(m_local_vid, m_local_eid, m_local_fid, m_global_cid, new_hash); } +std::string Tuple::to_string() const +{ + return std::to_string(m_local_vid) + " " + std::to_string(m_local_eid) + " " + + std::to_string(m_local_fid) + " " + std::to_string(m_global_cid) + " " + + std::to_string(m_hash); +} + } // namespace wmtk diff --git a/src/wmtk/Tuple.hpp b/src/wmtk/Tuple.hpp index bcee43587a..e0c7444e66 100644 --- a/src/wmtk/Tuple.hpp +++ b/src/wmtk/Tuple.hpp @@ -69,5 +69,10 @@ class Tuple bool is_null() const; Tuple with_updated_hash(long new_hash) const; + std::string to_string() const; + long get_local_vid() const { return m_local_vid; } + long get_local_eid() const { return m_local_eid; } + long get_local_fid() const { return m_local_fid; } + long get_global_cid() const { return m_global_cid; } }; } // namespace wmtk diff --git a/tests/components/CMakeLists.txt b/tests/components/CMakeLists.txt index 98057181e7..6b6ba89ed2 100644 --- a/tests/components/CMakeLists.txt +++ b/tests/components/CMakeLists.txt @@ -5,6 +5,7 @@ set(TEST_SOURCES test_component_mesh_info.cpp test_component_output.cpp test_component_isotropic_remeshing.cpp + test_component_extract_subset.cpp #test_smoothing.cpp test_component_regular_space.cpp ) diff --git a/tests/components/test_component_extract_subset.cpp b/tests/components/test_component_extract_subset.cpp new file mode 100644 index 0000000000..99fb22367c --- /dev/null +++ b/tests/components/test_component_extract_subset.cpp @@ -0,0 +1,648 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "../tools/DEBUG_TriMesh.hpp" +#include "../tools/TetMesh_examples.hpp" +#include "../tools/TriMesh_examples.hpp" +#include "wmtk_components/delaunay/internal/delaunay_2d.hpp" +#include "wmtk_components/delaunay/internal/delaunay_3d.hpp" +// #include +// #include + +long test_size_calculation(long n) +{ + // For a mesh with n faces/tets, there are 2^n different configurations of tag assignments, so 2^n possible subsets + // we want to test all of them, according to Coupon Collector's Problem, https://en.wikipedia.org/wiki/Coupon_collector%27s_problem + // the expected number of trials to collect them all is \Theta(2 ^n log(2^n)) = 2^n * n * log(2) + return long(ceil(pow(2, n) * n * log(2)) + 1); +} + +bool is_valid_mesh(const wmtk::TriMesh& tm) +{ + return true; +} + +bool is_connected(std::map>& connections) +{ + std::set visited_vertices; + std::stack stack; + stack.push(connections.begin()->first); + while (!stack.empty()) { + long current_vertex = stack.top(); + stack.pop(); + if (visited_vertices.count(current_vertex) == 0) { + visited_vertices.insert(current_vertex); + for (long neighbor : connections[current_vertex]) { + stack.push(neighbor); + } + } + } + return visited_vertices.size() == connections.size(); +} + +std::map> get_connection(const wmtk::TriMesh& tm, std::set& index_set) +{ + std::vector edges = tm.get_all(wmtk::PrimitiveType::Edge); + std::map> connections; + for (long edgeindex : index_set) { + wmtk::Tuple edgeTuple = edges[edgeindex]; + std::vector edgeVertexList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::edge(edgeTuple), + wmtk::PrimitiveType::Vertex); + std::vector vertex_index; + vertex_index.reserve(2); + for (wmtk::Tuple t : edgeVertexList) { + vertex_index.push_back(wmtk::components::internal::find_vertex_index(tm, t)); + } + + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) { + if (i != j) { + connections[vertex_index[i]].insert(vertex_index[j]); + } + } + } + // long v1 = wmtk::components::internal::find_vertex_index(tm, edgeVertexList[0]); + // long v2 = wmtk::components::internal::find_vertex_index(tm, edgeVertexList[1]); + // if (!connections.count(v1)) { + // std::vector nodes; + // nodes.push_back(v2); + // connections[v1] = nodes; + // } else + // connections[v1].push_back(v2); + // if (!connections.count(v2)) { + // std::vector nodes; + // nodes.push_back(v1); + // connections[v2] = nodes; + // } else + // connections[v2].push_back(v1); + } + return connections; +} + +std::map> get_connection_3d(const wmtk::TetMesh& tm, std::set& index_set) +{ + std::vector faces = tm.get_all(wmtk::PrimitiveType::Face); + std::map> connections; + for (long faceindex : index_set) { + wmtk::Tuple faceTuple = faces[faceindex]; + std::vector faceVertexList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::face(faceTuple), + wmtk::PrimitiveType::Vertex); + std::vector vertex_index; + vertex_index.reserve(3); + for (wmtk::Tuple t : faceVertexList) { + vertex_index.push_back(wmtk::components::internal::find_vertex_index(tm, t)); + } + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + if (i != j) connections[vertex_index[i]].insert(vertex_index[j]); + } + } + } + return connections; +} + +// Reference: https://www.geeksforgeeks.org/determining-topology-formed-in-a-graph/ +bool is_circle(const wmtk::TriMesh& tm, std::set index_set) +{ + std::map> connections = get_connection(tm, index_set); + if (index_set.size() != connections.size()) return false; + bool isRing = all_of(connections.begin(), connections.end(), [](auto& nodes) { + return nodes.second.size() == 2; + }); + bool connected = is_connected(connections); + return isRing && connected; +} + +// Reference: https://www.geeksforgeeks.org/determining-topology-formed-in-a-graph/ +bool is_line(const wmtk::TriMesh& tm, std::set index_set) +{ + if (index_set.size() == 1) return true; + std::map> connections = get_connection(tm, index_set); + if (index_set.size() != connections.size() - 1) return false; + long deg1 = 0, deg2 = 0; + for (auto& nodes : connections) { + if (nodes.second.size() == 1) + deg1++; + else if (nodes.second.size() == 2) + deg2++; + else + return false; + } + return deg1 == 2 && deg2 == connections.size() - 2; +} + +bool is_manifold_2d(const wmtk::TriMesh& tm) +{ + std::map> vertexLinkEdges; + std::vector faces = tm.get_all(wmtk::PrimitiveType::Face); + std::vector vertices = tm.get_all(wmtk::PrimitiveType::Vertex); + for (long vid = 0; vid < vertices.size(); ++vid) { + // wmtk::simplex::SimplexCollection sc = wmtk::simplex::top_dimension_cofaces(tm, + // wmtk::Simplex::vertex(vertices[vid])); + + // find all faces that are connected to the vertex + // TODO: replace adj_faces_of_vertex() with better implementation + std::vector adj_faces = wmtk::components::internal::adj_faces_of_vertex(tm, vid); + for (long fid : adj_faces) { + wmtk::Tuple faceTuple = faces[fid]; + std::vector edgeList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::face(faceTuple), + wmtk::PrimitiveType::Edge); + for (wmtk::Tuple edgeTuple : edgeList) { + std::vector edgeVertexList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::edge(edgeTuple), + wmtk::PrimitiveType::Vertex); + if (!tm.simplices_are_equal( + wmtk::Simplex::vertex(edgeVertexList[0]), + wmtk::Simplex::vertex(vertices[vid])) && + !tm.simplices_are_equal( + wmtk::Simplex::vertex(edgeVertexList[1]), + wmtk::Simplex::vertex(vertices[vid]))) { + // TODO: replace find_edge_index() with better implementation + vertexLinkEdges[vid].insert( + wmtk::components::internal::find_edge_index(tm, edgeTuple)); + } + } + } + } + + for (auto& [vid, edgeSet] : vertexLinkEdges) { + // for vertices on the boundary, the link needs to be a 1-ball, which is a line + if (tm.is_boundary(vertices[vid], wmtk::PrimitiveType::Vertex)) { + // std::cout << "Vertex " << vid << " is on the boundary." << std::endl; + // std::all_of(edgeSet.begin(), edgeSet.end(), [](long e) { + // std::cout << e << " "; + // return true; + // }); + // std::cout << std::endl; + if (!is_line(tm, edgeSet)) { + // std::cout << "Vertex " << vid << " doesn't have a line link." << std::endl; + return false; + } + } + // for vertices inside the mesh, the link needs to be a 1-sphere, which is a circle + else { + // std::cout << "Vertex " << vid << " is not on the boundary." << std::endl; + // std::all_of(edgeSet.begin(), edgeSet.end(), [](long e) { + // std::cout << e << " "; + // return true; + // }); + // std::cout << std::endl; + if (!is_circle(tm, edgeSet)) { + // std::cout << "Vertex " << vid << " doesn't have a circle link." << std::endl; + return false; + } + } + } + return true; +} + +bool is_disk(const wmtk::TetMesh& tm, std::set index_set) +{ + std::map> connections = get_connection_3d(tm, index_set); + + // display all items in connections + // std::cout << "Items in connections:" << std::endl; + // for (const auto& pair : connections) { + // std::cout << "Key: " << pair.first << ", Values: "; + // for (const auto& value : pair.second) { + // std::cout << value << " "; + // } + // std::cout << std::endl; + // } + std::vector all_faces = tm.get_all(wmtk::PrimitiveType::Face); + std::set edge_indexs_in_link; + std::set vertex_indexs_in_link; + for (long index : index_set) { + wmtk::Tuple face = all_faces[index]; + std::vector edgeList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::face(face), + wmtk::PrimitiveType::Edge); + std::vector vertexList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::face(face), + wmtk::PrimitiveType::Vertex); + for (wmtk::Tuple vertex : vertexList) { + vertex_indexs_in_link.insert(wmtk::components::internal::find_vertex_index(tm, vertex)); + } + for (wmtk::Tuple edge : edgeList) { + edge_indexs_in_link.insert(wmtk::components::internal::find_edge_index(tm, edge)); + } + } + // std::cout << "vertex_indexs_in_link.size() = " << vertex_indexs_in_link.size() << std::endl; + // std::cout << "edge_indexs_in_link.size() = " << edge_indexs_in_link.size() << std::endl; + // std::cout << "index_set.size() = " << index_set.size() << std::endl; + long euler_char = + vertex_indexs_in_link.size() - edge_indexs_in_link.size() + index_set.size() + 1; + // std::cout << "euler_char = " << euler_char << std::endl; + if (euler_char != 2) return false; + bool isRing = all_of(connections.begin(), connections.end(), [](auto& nodes) { + return nodes.second.size() >= 2; + }); + bool connected = is_connected(connections); + // std::cout << "isRIng = " << isRing << ", connected = " << connected << std::endl; + return isRing && connected; +} + +bool is_sphere(const wmtk::TetMesh& tm, std::set index_set) +{ + // long euler_char = tm.capacity(wmtk::PrimitiveType::Vertex) - + // tm.capacity(wmtk::PrimitiveType::Edge) + + // tm.capacity(wmtk::PrimitiveType::Face); + // if (euler_char != 2) return false; + std::map> connections = get_connection_3d(tm, index_set); + bool isRing = all_of(connections.begin(), connections.end(), [](auto& nodes) { + return nodes.second.size() >= 3; + }); + bool connected = is_connected(connections); + return isRing && connected; +} + +bool is_manifold_3d(const wmtk::TetMesh& tm) +{ + std::map> vertexLinkFaces; + std::vector tets = tm.get_all(wmtk::PrimitiveType::Tetrahedron); + std::vector faces = tm.get_all(wmtk::PrimitiveType::Face); + std::vector vertices = tm.get_all(wmtk::PrimitiveType::Vertex); + for (long vid = 0; vid < vertices.size(); ++vid) { + std::vector adj_tets = wmtk::components::internal::adj_tets_of_vertex(tm, vid); + for (long fid : adj_tets) { + wmtk::Tuple tetTuple = tets[fid]; + std::vector faceList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::tetrahedron(tetTuple), + wmtk::PrimitiveType::Face); + bool find_link_face = false; + for (wmtk::Tuple faceTuple : faceList) { + std::vector faceVertexList = wmtk::simplex::faces_single_dimension( + tm, + wmtk::Simplex::face(faceTuple), + wmtk::PrimitiveType::Vertex); + if (std::none_of(faceVertexList.begin(), faceVertexList.end(), [&](wmtk::Tuple t) { + return tm.simplices_are_equal( + wmtk::Simplex::vertex(t), + wmtk::Simplex::vertex(vertices[vid])); + })) { + vertexLinkFaces[vid].insert( + wmtk::components::internal::find_face_index(tm, faceTuple)); + find_link_face = true; + } + } + if (!find_link_face) throw std::runtime_error("Link face not found!"); + } + } + + for (auto& [vid, faceSet] : vertexLinkFaces) { + // for vertices on the boundary, the link needs to be a 2-ball, which is a disk + // if (tm.is_boundary(vertices[vid], wmtk::PrimitiveType::Vertex)) { + if (tm.is_boundary(vertices[vid], wmtk::PrimitiveType::Vertex)) { + // std::cout << "Vertex " << vid << " is on the boundary." << std::endl; + // std::all_of(faceSet.begin(), faceSet.end(), [](long e) { + // std::cout << e << " "; + // return true; + // }); + // std::cout << std::endl; + if (!is_disk(tm, faceSet)) { + // std::cout << "Vertex " << vid << " on the boundary doesn't have a disk link." + // << std::endl; + return false; + } + } + // for vertices inside the mesh, the link needs to be a 1-sphere, which is a circle + else { + // std::cout << "Vertex " << vid << " is not on the boundary." << std::endl; + // std::all_of(faceSet.begin(), faceSet.end(), [](long e) { + // std::cout << e << " "; + // return true; + // }); + // std::cout << std::endl; + if (!is_sphere(tm, faceSet)) { + // std::cout << "Vertex " << vid << " doesn't have a sphere link." << std::endl; + return false; + } + } + } + return true; +} + +void random_trimesh_test_executor(const wmtk::TriMesh& m, const unsigned long test_size) +{ + wmtk::tests::DEBUG_TriMesh tm = m; + long top_dimen_count = m.get_all(wmtk::PrimitiveType::Face).size(); + std::vector tag_vector(top_dimen_count, 0); + for (size_t i = 0; i < test_size; ++i) { + std::random_device rd{}; + std::mt19937 mt{rd()}; + std::uniform_int_distribution tag{0, 1}; + for (int j = 0; j < tag_vector.size(); ++j) { + tag_vector[j] = tag(mt); + } + if (std::reduce(tag_vector.begin(), tag_vector.end()) == 0) { + std::fill(tag_vector.begin(), tag_vector.end(), 0); + continue; + } + // std::cout << "Tag: "; + // std::all_of(tag_vector.begin(), tag_vector.end(), [](int i) { + // std::cout << i; + // return true; + // }); + std::unique_ptr new_tm = + wmtk::components::extract_subset(tm, tag_vector, false); + bool after = is_manifold_2d(*dynamic_cast(new_tm.get())); + // std::cout << "After: manifold = " << after << std::endl; + CHECK(after); + std::fill(tag_vector.begin(), tag_vector.end(), 0); + } +} + +long random_tetmesh_test_executor(const wmtk::TetMesh& m, const unsigned long test_size) +{ + wmtk::TetMesh tm = m; + long top_dimen_count = m.get_all(m.top_simplex_type()).size(); + std::vector tag_vector(top_dimen_count, 0); + for (size_t i = 0; i < test_size; ++i) { + std::random_device rd{}; + std::mt19937 mt{rd()}; + std::uniform_int_distribution tag{0, 1}; + for (int j = 0; j < tag_vector.size(); ++j) { + tag_vector[j] = tag(mt); + } + if (std::reduce(tag_vector.begin(), tag_vector.end()) == 0) { + std::fill(tag_vector.begin(), tag_vector.end(), 0); + continue; + } + // std::cout << "Tag: "; + // std::all_of(tag_vector.begin(), tag_vector.end(), [](int i) { + // std::cout << i; + // return true; + // }); + std::unique_ptr new_tm = wmtk::components::extract_subset(tm, tag_vector, true); + bool after = is_manifold_3d(*dynamic_cast(new_tm.get())); + // std::cout << "After: manifold = " << after << std::endl; + CHECK(after); + if (!after) { + paraviewo::VTUWriter writer; + // writer.write_mesh("man_ext_3d_random.vtu", vertices, faces); + return -1; + } + } + std::fill(tag_vector.begin(), tag_vector.end(), 0); + return 0; +} + +TEST_CASE("2d_9tri_with_a_hole_test_case", "[components][extract_subset][2D]") +{ + wmtk::tests::DEBUG_TriMesh tm = wmtk::tests::nine_triangles_with_a_hole(); + const unsigned long test_size = test_size_calculation(tm.capacity(wmtk::PrimitiveType::Face)); + std::cout << "Test size: " << test_size << std::endl; + random_trimesh_test_executor(tm, test_size); +} + +TEST_CASE("component_3+4_test_case", "[components][extract_subset][2D][manual]") +{ + wmtk::tests::DEBUG_TriMesh tm; + wmtk::RowVectors3l tris; + tris.resize(46, 3); + tris.row(0) << 0, 1, 2; + tris.row(1) << 0, 2, 3; + tris.row(2) << 1, 2, 4; + tris.row(3) << 2, 3, 4; + tris.row(4) << 0, 3, 5; + tris.row(5) << 3, 4, 5; + tris.row(6) << 0, 5, 6; + tris.row(7) << 5, 6, 7; + tris.row(8) << 4, 5, 8; + tris.row(9) << 0, 6, 7; + tris.row(10) << 5, 7, 8; + tris.row(11) << 0, 7, 9; + tris.row(12) << 7, 8, 9; + tris.row(13) << 0, 9, 10; + tris.row(14) << 9, 10, 12; + tris.row(15) << 9, 12, 11; + tris.row(16) << 9, 8, 11; + tris.row(17) << 0, 10, 13; + tris.row(18) << 10, 12, 13; + tris.row(19) << 13, 12, 15; + tris.row(20) << 12, 15, 16; + tris.row(21) << 11, 12, 16; + tris.row(22) << 11, 16, 17; + tris.row(23) << 15, 16, 17; + tris.row(24) << 13, 15, 17; + tris.row(25) << 0, 13, 14; + tris.row(26) << 13, 14, 17; + tris.row(27) << 0, 14, 20; + tris.row(28) << 14, 18, 20; + tris.row(29) << 14, 18, 17; + tris.row(30) << 17, 18, 19; + tris.row(31) << 0, 20, 24; + tris.row(32) << 20, 18, 24; + tris.row(33) << 18, 19, 24; + tris.row(34) << 19, 21, 24; + tris.row(35) << 0, 24, 23; + tris.row(36) << 24, 21, 23; + tris.row(37) << 21, 22, 23; + tris.row(38) << 0, 23, 27; + tris.row(39) << 27, 23, 26; + tris.row(40) << 23, 22, 26; + tris.row(41) << 0, 27, 29; + tris.row(42) << 29, 27, 28; + tris.row(43) << 28, 27, 26; + tris.row(44) << 0, 25, 29; + tris.row(45) << 25, 28, 29; + tm.initialize(tris); + + std::vector tag_vector(tm.get_all(wmtk::PrimitiveType::Face).size(), 0); + std::vector id = {0, 1, 2, 3, 5, 6, 7, 8, 10, 11, 12, 25, 26, 29, + 30, 31, 32, 33, 34, 36, 37, 38, 39, 40, 42, 43, 44, 45}; + for (int i : id) tag_vector[i] = 1; + std::cout << "Before: nb_vertex = " << tm.capacity(wmtk::PrimitiveType::Vertex) + << std::endl; // 30 -> 25 + std::cout << "Before: nb_face = " << tm.capacity(wmtk::PrimitiveType::Face) + << std::endl; // 46 -> 28 + std::unique_ptr new_tm = wmtk::components::extract_subset(tm, tag_vector, false); + + CHECK(is_valid_mesh(*dynamic_cast(new_tm.get()))); + CHECK(new_tm->capacity(wmtk::PrimitiveType::Vertex) == 31); + CHECK(new_tm->capacity(wmtk::PrimitiveType::Face) == 28); + bool after = is_manifold_2d(*dynamic_cast(new_tm.get())); + // std::cout << "After: manifold = " << after << std::endl; + CHECK(after); + // new_tm.print_vf(); +} + +TEST_CASE("random_2D_test_from_manext_branch", "[components][extract_subset][2D][random]") +{ + unsigned int nb_points = 20; // 20 + double range = 10.0; + const size_t tagass_loop = 100; // 100 + const size_t pntgen_loop = 6; // 10 + const double prob = 0.2; + + // test for 10 iterations, each with 10 more vertices, so 10~100 + for (size_t i = 0; i < pntgen_loop; ++i) { + wmtk::TriMesh tm; + wmtk::RowVectors3l tris; + wmtk::RowVectors2d points(nb_points, 2); + std::random_device rd{}; + std::mt19937 gen(rd()); + std::uniform_real_distribution dis(0, range); + for (size_t j = 0; j < nb_points; ++j) { + // generate 2 random doubles between 0 and the given range + points.row(j) << dis(gen), dis(gen); + } + + Eigen::MatrixXd vertices; + Eigen::MatrixXi faces; + std::tie(vertices, faces) = wmtk::components::internal::delaunay_2d(points); + unsigned int nb_triangles = faces.rows(); + unsigned int nb_vertices = vertices.rows(); + std::cout << "Man-ext 2D test: total tri num=" << nb_triangles << "\n"; + tris.resize(nb_triangles, 3); + for (unsigned int j = 0; j < nb_triangles; ++j) { + tris.row(j) << faces(j, 0), faces(j, 1), faces(j, 2); + } + tm.initialize(tris); + wmtk::mesh_utils::set_matrix_attribute( + vertices, + "position", + wmtk::PrimitiveType::Vertex, + tm); + random_trimesh_test_executor(tm, tagass_loop); + nb_points += 10; + range += 10.0; + } +} + +TEST_CASE("2_non_manifold_vertices", "[components][extract_subset][3D][manual][3]") +{ + wmtk::TetMesh tm; + wmtk::RowVectors tets; + tets.resize(3, 4); + tets.row(0) << 0, 1, 2, 3; + tets.row(1) << 0, 2, 3, 4; + tets.row(2) << 1, 4, 5, 6; + tm.initialize(tets); + std::cout << "Before: # of vertices = " << tm.get_all(wmtk::PrimitiveType::Vertex).size() + << std::endl; + bool before = is_manifold_3d(tm); + std::cout << "Before: manifold = " << before << std::endl; + std::vector tag_vector = {1, 1, 1}; + std::unique_ptr topo_tm = wmtk::components::extract_subset(tm, tag_vector, false); + std::cout << "After, # of vertices = " << topo_tm->get_all(wmtk::PrimitiveType::Vertex).size() + << std::endl; + bool after = is_manifold_3d(*dynamic_cast(topo_tm.get())); + std::cout << "After: manifold = " << after << std::endl; + CHECK(after); +} + +TEST_CASE("2_non_manifold_edges", "[components][extract_subset][3D][manual][3]") +{ + wmtk::TetMesh tm; + wmtk::RowVectors tets; + tets.resize(3, 4); + tets.row(0) << 0, 1, 2, 3; + tets.row(1) << 0, 2, 3, 4; + tets.row(2) << 1, 3, 4, 5; + tm.initialize(tets); + bool before = is_manifold_3d(tm); + std::cout << "Before: manifold = " << before << std::endl; + std::vector tag_vector = {1, 1, 1}; + std::unique_ptr topo_tm = wmtk::components::extract_subset(tm, tag_vector, false); + std::cout << "After: # of vertices = " << topo_tm->get_all(wmtk::PrimitiveType::Vertex).size() + << std::endl; + bool after = is_manifold_3d(*dynamic_cast(topo_tm.get())); + std::cout << "; After: manifold = " << after << std::endl; + CHECK(after); +} + +TEST_CASE("six_cycle_tets", "[components][extract_subset][3D][manual][6]") +{ + wmtk::TetMesh tm = wmtk::tests_3d::six_cycle_tets(); + const unsigned long test_size = 10; // total cases + std::vector tag_vector(tm.capacity(wmtk::PrimitiveType::Tetrahedron), 0); + for (size_t i = 0; i < test_size; ++i) { + std::mt19937 mt{i}; + std::uniform_int_distribution tag{0, 1}; + for (int j = 0; j < tag_vector.size(); ++j) { + tag_vector[j] = tag(mt); + } + if (std::reduce(tag_vector.begin(), tag_vector.end()) == 0) { + std::fill(tag_vector.begin(), tag_vector.end(), 0); + continue; + } + // std::all_of(tag_vector.begin(), tag_vector.end(), [](int i) { + // std::cout << i << " "; + // return true; + // }); + std::unique_ptr topo_tm = + wmtk::components::extract_subset(tm, tag_vector, false); + bool after = is_manifold_3d(*(dynamic_cast(topo_tm.get()))); + std::cout << "After: manifold = " << after << std::endl; + CHECK(after); + std::fill(tag_vector.begin(), tag_vector.end(), 0); + } +} + +TEST_CASE("random_3D_test", "[components][extract_subset][2D][random]") +{ + std::cout << "Random 3D test" << std::endl; + unsigned int nb_points = 7; // 20 + double range = 10.0; + long tagass_loop = 40; // 100 + const size_t pntgen_loop = 5; // 10 + const double prob = 0.2; + + // test for 10 iterations, each with 10 more vertices, so 10~100 + for (size_t i = 0; i < pntgen_loop; ++i) { + wmtk::TetMesh tm; + wmtk::RowVectors4l tets; + wmtk::RowVectors3d points(nb_points, 3); + std::random_device rd{}; + std::mt19937 gen(rd()); + std::uniform_real_distribution dis(0, range); + for (size_t j = 0; j < nb_points; ++j) { + // generate 3 random doubles between 0 and the given range + points.row(j) << dis(gen), dis(gen), dis(gen); + } + + Eigen::MatrixXd vertices; + Eigen::MatrixXi cells; + std::tie(vertices, cells) = wmtk::components::internal::delaunay_3d(points); + unsigned int nb_tets = cells.rows(); + unsigned int nb_vertices = vertices.rows(); + std::cout << "Man-ext 3D test: total tet num=" << nb_tets << "\n"; + tets.resize(nb_tets, 4); + for (unsigned int j = 0; j < nb_tets; ++j) { + tets.row(j) << cells(j, 0), cells(j, 1), cells(j, 2), cells(j, 3); + } + tm.initialize(tets); + wmtk::mesh_utils::set_matrix_attribute( + vertices, + "position", + wmtk::PrimitiveType::Vertex, + tm); + tagass_loop = test_size_calculation(nb_tets); + long ret = random_tetmesh_test_executor(tm, tagass_loop); + if (ret == -1) { + paraviewo::VTUWriter writer; + writer.write_mesh("man_ext_3d_random.vtu", vertices, cells); + throw std::runtime_error("Manifold test failed!"); + } + // nb_points += 10; + range += 10.0; + } +} \ No newline at end of file