diff --git a/src/wmtk/simplex/cofaces_single_dimension.cpp b/src/wmtk/simplex/cofaces_single_dimension.cpp index 9e89c63604..921be0565d 100644 --- a/src/wmtk/simplex/cofaces_single_dimension.cpp +++ b/src/wmtk/simplex/cofaces_single_dimension.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -15,34 +16,34 @@ namespace wmtk::simplex { -namespace { - -// wrapper for calling the internal function boundary_with_preserved_face_tuples -std::vector boundary_with_preserved_face_tuples( - const Mesh& mesh, - const std::vector& tuples, - const PrimitiveType pt, - const PrimitiveType coface_pt) -{ - std::vector r; - for (const Tuple& t : tuples) { - const auto tups = - wmtk::simplex::internal::boundary_with_preserved_face_tuples(mesh, t, pt, coface_pt); - r.insert(r.end(), tups.begin(), tups.end()); - } - - const PrimitiveType boundary_pt = get_primitive_type_from_id(get_primitive_type_id(pt) - 1); - std::sort(r.begin(), r.end(), [&](const Tuple& a, const Tuple& b) -> bool { - return utils::SimplexComparisons::less(mesh, boundary_pt, a, b); - }); - auto last = std::unique(r.begin(), r.end(), [&](const Tuple& a, const Tuple& b) -> bool { - return utils::SimplexComparisons::equal(mesh, boundary_pt, a, b); - }); - r.erase(last, r.end()); - - return r; -} -} // namespace +//namespace { +// +//// wrapper for calling the internal function boundary_with_preserved_face_tuples +//std::vector boundary_with_preserved_face_tuples( +// const Mesh& mesh, +// const std::vector& tuples, +// const PrimitiveType pt, +// const PrimitiveType coface_pt) +//{ +// std::vector r; +// for (const Tuple& t : tuples) { +// const auto tups = +// wmtk::simplex::internal::boundary_with_preserved_face_tuples(mesh, t, pt, coface_pt); +// r.insert(r.end(), tups.begin(), tups.end()); +// } +// +// const PrimitiveType boundary_pt = get_primitive_type_from_id(get_primitive_type_id(pt) - 1); +// std::sort(r.begin(), r.end(), [&](const Tuple& a, const Tuple& b) -> bool { +// return utils::SimplexComparisons::less(mesh, boundary_pt, a, b); +// }); +// auto last = std::unique(r.begin(), r.end(), [&](const Tuple& a, const Tuple& b) -> bool { +// return utils::SimplexComparisons::equal(mesh, boundary_pt, a, b); +// }); +// r.erase(last, r.end()); +// +// return r; +//} +//} // namespace std::vector cofaces_single_dimension_tuples( @@ -50,22 +51,9 @@ std::vector cofaces_single_dimension_tuples( const Simplex& my_simplex, PrimitiveType cofaces_type) { - std::vector tuples; - if (my_simplex.primitive_type() == cofaces_type) { - tuples = {my_simplex.tuple()}; - return tuples; - } - - tuples = top_dimension_cofaces_tuples(mesh, my_simplex); - + SimplexCollection cof = cofaces_single_dimension(mesh, my_simplex, cofaces_type); - assert(my_simplex.primitive_type() < cofaces_type); - auto range = wmtk::utils::primitive_range(mesh.top_simplex_type(), cofaces_type); - range.pop_back(); - for (const auto& pt : range) { - tuples = boundary_with_preserved_face_tuples(mesh, tuples, pt, my_simplex.primitive_type()); - } - return tuples; + return cof.simplex_vector_tuples(cofaces_type); } @@ -74,10 +62,9 @@ std::vector cofaces_single_dimension_simplices( const Simplex& simplex, PrimitiveType cofaces_type) { - return utils::tuple_vector_to_homogeneous_simplex_vector( - mesh, - cofaces_single_dimension_tuples(mesh, simplex, cofaces_type), - cofaces_type); + SimplexCollection cof = cofaces_single_dimension(mesh, simplex, cofaces_type); + + return cof.simplex_vector(); } SimplexCollection cofaces_single_dimension( @@ -86,14 +73,166 @@ SimplexCollection cofaces_single_dimension( PrimitiveType cofaces_type, bool sort_and_clean) { - SimplexCollection collection( - mesh, - cofaces_single_dimension_simplices(mesh, my_simplex, cofaces_type)); + simplex::SimplexCollection cofaces(mesh); + + cofaces_single_dimension(cofaces, my_simplex, cofaces_type, sort_and_clean); + + return cofaces; +} + +void cofaces_single_dimension_general( + SimplexCollection& collection, + const Simplex& my_simplex, + PrimitiveType cofaces_type, + bool sort_and_clean) +{ + if (my_simplex.primitive_type() == cofaces_type) { + collection.add(my_simplex); + return; + } + + std::vector tuples = top_dimension_cofaces_tuples(collection.mesh(), my_simplex); + + assert(my_simplex.primitive_type() < cofaces_type); + + for (const Tuple& t : tuples) { + simplices_preserving_primitive_types( + collection, + t, + collection.mesh().top_simplex_type(), + my_simplex.primitive_type(), + cofaces_type); + } + + if (sort_and_clean) { + collection.sort_and_clean(); + } +} + +namespace { +/** + * @brief Special case of a TriMesh where we want to get the edges adjacent to a vertex. + */ +void cofaces_single_dimension_tri_vertex_edges( + SimplexCollection& collection, + const Simplex& my_simplex, + bool sort_and_clean) +{ + assert(my_simplex.primitive_type() == PrimitiveType::Vertex); + + const TriMesh& mesh = static_cast(collection.mesh()); + + const Tuple& t_in = my_simplex.tuple(); + assert(mesh.is_valid_slow(t_in)); + + Tuple t = t_in; + do { + collection.add(PrimitiveType::Edge, t); + + if (mesh.is_boundary_edge(t)) { + break; + } + t = mesh.switch_tuples(t, {PrimitiveType::Triangle, PrimitiveType::Edge}); + } while (t != t_in); + + if (t == t_in && !mesh.is_boundary_edge(t)) { + return; + } + + t = mesh.switch_edge(t_in); + collection.add(PrimitiveType::Edge, t); + + if (mesh.is_boundary_edge(t)) { + return; + } + t = mesh.switch_tuples(t, {PrimitiveType::Triangle, PrimitiveType::Edge}); + + do { + collection.add(PrimitiveType::Edge, t); + + if (mesh.is_boundary_edge(t)) { + break; + } + t = mesh.switch_tuples(t, {PrimitiveType::Triangle, PrimitiveType::Edge}); + } while (true); + if (sort_and_clean) { collection.sort_and_clean(); } +} + +void cofaces_single_dimension_tet_edge_triangles( + SimplexCollection& collection, + const Simplex& my_simplex, + bool sort_and_clean) +{ + assert(my_simplex.primitive_type() == PrimitiveType::Vertex); + + const TetMesh& mesh = static_cast(collection.mesh()); + + const Tuple& t_in = my_simplex.tuple(); + assert(mesh.is_valid_slow(t_in)); + Tuple t = t_in; + do { + collection.add(PrimitiveType::Triangle, t); + + if (mesh.is_boundary_face(t)) { + break; + } + t = mesh.switch_tuples(t, {PrimitiveType::Tetrahedron, PrimitiveType::Triangle}); + } while (t != t_in); + + if (t == t_in && !mesh.is_boundary_face(t)) { + return; + } + + t = mesh.switch_face(t_in); + collection.add(PrimitiveType::Triangle, t); + + if (mesh.is_boundary_face(t)) { + return; + } + t = mesh.switch_tuples(t, {PrimitiveType::Tetrahedron, PrimitiveType::Triangle}); + + do { + collection.add(PrimitiveType::Triangle, t); + + if (mesh.is_boundary_face(t)) { + break; + } + t = mesh.switch_tuples(t, {PrimitiveType::Tetrahedron, PrimitiveType::Triangle}); + } while (true); + + if (sort_and_clean) { + collection.sort_and_clean(); + } +} +} // namespace + +void cofaces_single_dimension( + SimplexCollection& collection, + const Simplex& my_simplex, + PrimitiveType cofaces_type, + bool sort_and_clean) +{ + assert(my_simplex.primitive_type() < cofaces_type); + + const Mesh& m = collection.mesh(); + if (m.top_simplex_type() == PrimitiveType::Triangle && + my_simplex.primitive_type() == PrimitiveType::Vertex && + cofaces_type == PrimitiveType::Edge) { + cofaces_single_dimension_tri_vertex_edges(collection, my_simplex, sort_and_clean); + return; + } + + if (m.top_simplex_type() == PrimitiveType::Tetrahedron && + my_simplex.primitive_type() == PrimitiveType::Edge && + cofaces_type == PrimitiveType::Triangle) { + cofaces_single_dimension_tet_edge_triangles(collection, my_simplex, sort_and_clean); + return; + } - return collection; + cofaces_single_dimension_general(collection, my_simplex, cofaces_type, sort_and_clean); } } // namespace wmtk::simplex diff --git a/src/wmtk/simplex/cofaces_single_dimension.hpp b/src/wmtk/simplex/cofaces_single_dimension.hpp index 012162f440..58506af8c5 100644 --- a/src/wmtk/simplex/cofaces_single_dimension.hpp +++ b/src/wmtk/simplex/cofaces_single_dimension.hpp @@ -26,6 +26,12 @@ SimplexCollection cofaces_single_dimension( PrimitiveType cofaces_type, bool sort_and_clean = true); +void cofaces_single_dimension( + SimplexCollection& collection, + const Simplex& my_simplex, + PrimitiveType cofaces_type, + bool sort_and_clean = true); + std::vector cofaces_single_dimension_tuples( const Mesh& mesh, const Simplex& my_simplex, diff --git a/src/wmtk/simplex/tuples_preserving_primitive_types.cpp b/src/wmtk/simplex/tuples_preserving_primitive_types.cpp index 92e2d8fdb6..f8ab75697c 100644 --- a/src/wmtk/simplex/tuples_preserving_primitive_types.cpp +++ b/src/wmtk/simplex/tuples_preserving_primitive_types.cpp @@ -12,7 +12,9 @@ std::vector tuples_preserving_primitive_types( const PrimitiveType simplex_ptype, const PrimitiveType face_ptype) { - std::vector switch_tuple_types = + // old and slow implementation + + const std::vector switch_tuple_types = wmtk::utils::primitive_range(simplex_ptype, face_ptype); Tuple t_iter = t; @@ -28,4 +30,33 @@ std::vector tuples_preserving_primitive_types( return intersection_tuples; } +void simplices_preserving_primitive_types( + SimplexCollection& collection, + const Tuple& t, + const PrimitiveType simplex_ptype, + const PrimitiveType face_ptype, + const PrimitiveType pt_return) +{ + if (simplex_ptype == pt_return) { + collection.add(pt_return, t); + return; + } + + const Mesh& mesh = collection.mesh(); + + const int8_t pt0_id = get_primitive_type_id(face_ptype); + const int8_t pt1_id = get_primitive_type_id(simplex_ptype); + assert(pt0_id <= pt1_id); + + Tuple t_iter = t; + + do { + collection.add(pt_return, t_iter); + for (int8_t i = pt0_id + 1; i < pt1_id; ++i) { + const PrimitiveType pt = get_primitive_type_from_id(i); + t_iter = mesh.switch_tuple(t_iter, pt); + } + } while (t != t_iter); +} + } // namespace wmtk::simplex \ No newline at end of file diff --git a/src/wmtk/simplex/tuples_preserving_primitive_types.hpp b/src/wmtk/simplex/tuples_preserving_primitive_types.hpp index 3f48a20d96..204dd0f91e 100644 --- a/src/wmtk/simplex/tuples_preserving_primitive_types.hpp +++ b/src/wmtk/simplex/tuples_preserving_primitive_types.hpp @@ -2,6 +2,7 @@ #include #include +#include namespace wmtk::simplex { @@ -19,4 +20,19 @@ std::vector tuples_preserving_primitive_types( const Tuple& t, const PrimitiveType ptype1, const PrimitiveType ptype2); + +/** + * @brief Compute all simplices that contain simplex(ptype1, t) and that are contained by + * simplex(ptype2, t). Simplices can appear multiple times. + * + * If ptype1 and ptype2 are the same, only the input tuple is returned. + * + * The return tuples are guaranteed to contain both input simplices. + */ +void simplices_preserving_primitive_types( + SimplexCollection& collection, + const Tuple& t, + const PrimitiveType ptype1, + const PrimitiveType ptype2, + const PrimitiveType pt_return); } // namespace wmtk::simplex \ No newline at end of file diff --git a/tests/multimesh/test_multi_mesh.cpp b/tests/multimesh/test_multi_mesh.cpp index 43d3c054ec..e721a5bea6 100644 --- a/tests/multimesh/test_multi_mesh.cpp +++ b/tests/multimesh/test_multi_mesh.cpp @@ -1052,10 +1052,10 @@ TEST_CASE("test_split_multi_mesh", "[multimesh][2D]") CHECK(parent.fv_from_fid(4) == Vector3l(3, 5, 0)); CHECK(child0.fv_from_fid(1) == Vector3l(3, 1, 2)); CHECK(child0.fv_from_fid(2) == Vector3l(0, 3, 2)); - CHECK(child1.fv_from_fid(4) == Vector3l(4, 1, 2)); - CHECK(child1.fv_from_fid(5) == Vector3l(0, 4, 2)); - CHECK(child1.fv_from_fid(2) == Vector3l(3, 1, 4)); - CHECK(child1.fv_from_fid(3) == Vector3l(3, 4, 0)); + CHECK(child1.fv_from_fid(2) == Vector3l(4, 1, 2)); + CHECK(child1.fv_from_fid(3) == Vector3l(0, 4, 2)); + CHECK(child1.fv_from_fid(4) == Vector3l(3, 1, 4)); + CHECK(child1.fv_from_fid(5) == Vector3l(3, 4, 0)); CHECK(child2.fv_from_fid(2) == Vector3l(0, 2, 4)); CHECK(child2.fv_from_fid(3) == Vector3l(7, 1, 2)); CHECK(child2.fv_from_fid(4) == Vector3l(0, 7, 2)); @@ -1088,8 +1088,8 @@ TEST_CASE("test_split_multi_mesh", "[multimesh][2D]") CHECK(child0.fv_from_fid(3) == Vector3l(0, 4, 2)); CHECK(child0.fv_from_fid(4) == Vector3l(4, 3, 2)); - CHECK(child1.fv_from_fid(4) == Vector3l(4, 1, 2)); - CHECK(child1.fv_from_fid(2) == Vector3l(3, 1, 4)); + CHECK(child1.fv_from_fid(2) == Vector3l(4, 1, 2)); + CHECK(child1.fv_from_fid(4) == Vector3l(3, 1, 4)); CHECK(child1.fv_from_fid(8) == Vector3l(0, 5, 2)); CHECK(child1.fv_from_fid(9) == Vector3l(5, 4, 2)); CHECK(child1.fv_from_fid(6) == Vector3l(3, 5, 0)); @@ -1337,4 +1337,4 @@ TEST_CASE("test_deregister_child_mesh", "[multimesh]") CHECK(c1_mul_manager.is_root()); CHECK_FALSE(c0_mul_manager.is_root()); } -} \ No newline at end of file +} diff --git a/tests/test_performance.cpp b/tests/test_performance.cpp index 800a56946e..1c0b998113 100644 --- a/tests/test_performance.cpp +++ b/tests/test_performance.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -528,6 +529,38 @@ TEST_CASE("navigation_performance_tri", "[simplex][performance][.]") } logger().info("sum = {}", counter); } + + // cofaces tests + counter = 0; + if (false) { + POLYSOLVE_SCOPED_STOPWATCH("coface edges (sort)", logger()); + for (size_t i = 0; i < n_repetitions; ++i) { + for (const Tuple& t : vertices) { + const auto neighs = simplex::cofaces_single_dimension( + m, + simplex::Simplex::vertex(m, t), + PrimitiveType::Edge); + counter += neighs.size(); + } + } + logger().info("sum = {}", counter); + } + + counter = 0; + if (false) { + POLYSOLVE_SCOPED_STOPWATCH("coface edges (no sort)", logger()); + for (size_t i = 0; i < n_repetitions; ++i) { + for (const Tuple& t : vertices) { + const auto neighs = simplex::cofaces_single_dimension( + m, + simplex::Simplex::vertex(m, t), + PrimitiveType::Edge, + false); + counter += neighs.size(); + } + } + logger().info("sum = {}", counter); + } } TEST_CASE("navigation_performance_tet", "[simplex][performance][.]") @@ -603,4 +636,24 @@ TEST_CASE("navigation_performance_tet", "[simplex][performance][.]") } logger().info("sum = {}", counter); } + + // cofaces + { + const auto edges = m.get_all(PrimitiveType::Edge); + + counter = 0; + if (false) { + POLYSOLVE_SCOPED_STOPWATCH("coface triangles of edges (sort)", logger()); + for (size_t i = 0; i < n_repetitions; ++i) { + for (const Tuple& t : edges) { + const auto neighs = simplex::cofaces_single_dimension( + m, + simplex::Simplex::edge(m, t), + PrimitiveType::Triangle); + counter += neighs.size(); + } + } + logger().info("sum = {}", counter); + } + } } \ No newline at end of file diff --git a/tests/test_simplex_collection.cpp b/tests/test_simplex_collection.cpp index 1d6a6ea81e..da360580c5 100644 --- a/tests/test_simplex_collection.cpp +++ b/tests/test_simplex_collection.cpp @@ -1548,11 +1548,8 @@ TEST_CASE("simplex_cofaces_single_dimension", "[simplex_collection][2D]") REQUIRE(tc.size() == 3); SimplexCollection sc( m, - simplex::utils::tuple_vector_to_homogeneous_simplex_vector( - m, - tc, - PrimitiveType::Triangle)); - sc.sort(); + simplex::utils::tuple_vector_to_homogeneous_simplex_vector(m, tc, PrimitiveType::Edge)); + sc.sort_and_clean(); const auto& cells = sc.simplex_vector();