diff --git a/components/CMakeLists.txt b/components/CMakeLists.txt index e02b58fd63..22d0f6757c 100644 --- a/components/CMakeLists.txt +++ b/components/CMakeLists.txt @@ -39,6 +39,8 @@ add_subdirectory("CDT/wmtk/components/CDT") add_subdirectory("edge_insertion/wmtk/components/edge_insertion") +add_subdirectory("simplicial_embedding/wmtk/components/simplicial_embedding") + # add_component("export_cache") # add_component("import_cache") # add_component("mesh_info") diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/CMakeLists.txt b/components/simplicial_embedding/wmtk/components/simplicial_embedding/CMakeLists.txt new file mode 100644 index 0000000000..e631d58e24 --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/CMakeLists.txt @@ -0,0 +1,20 @@ +set(COMPONENT_NAME simplicial_embedding) +add_component(${COMPONENT_NAME}) + +if(NOT ${WMTK_ENABLE_COMPONENT_${COMPONENT_NAME}}) + return() +endif() + +set(SRC_FILES + SimplicialEmbeddingOptions.hpp + internal/SimplicialEmbedding.hpp + internal/SimplicialEmbedding.cpp + simplicial_embedding.hpp + simplicial_embedding.cpp + ) + + +target_sources(wmtk_${COMPONENT_NAME} PRIVATE ${SRC_FILES}) + +add_component_test(wmtk::${COMPONENT_NAME} + ${PROJECT_SOURCE_DIR}/components/tests/simplicial_embedding/test_simplicial_embedding.cpp) \ No newline at end of file diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/SimplicialEmbeddingOptions.hpp b/components/simplicial_embedding/wmtk/components/simplicial_embedding/SimplicialEmbeddingOptions.hpp new file mode 100644 index 0000000000..8f2ecbc701 --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/SimplicialEmbeddingOptions.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include + +namespace wmtk::components { + +struct SimplicialEmbeddingOptions +{ + /** + * All simplex dimensions must have an int64_t scalar attribute where the tags are stored. + */ + std::map tag_attributes; + /** + * The value that should be simplicially embedded. + */ + int64_t value; + /** + * Other attributes that should be processed with the default behavior. + */ + std::vector pass_through_attributes; + /** + * If false, this component forms a simplicial complex out of the tags, i.e., all faces of + * tagged simplices are also tagged. + */ + bool generate_simplicial_embedding = true; +}; + +} // namespace wmtk::components \ No newline at end of file diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.cpp b/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.cpp new file mode 100644 index 0000000000..59a8a2c167 --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.cpp @@ -0,0 +1,311 @@ +#include "SimplicialEmbedding.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace wmtk::components::internal { + +namespace { + +class TagAttribute +{ +public: + wmtk::attribute::Accessor m_accessor; + PrimitiveType m_ptype; + int64_t m_val; + + TagAttribute( + Mesh& m, + const attribute::MeshAttributeHandle& attribute, + PrimitiveType ptype, + int64_t val) + : m_accessor(m.create_accessor(attribute)) + , m_ptype(ptype) + , m_val(val) + {} + + TagAttribute(TagAttribute&) = delete; +}; +} // namespace + +SimplicialEmbedding::SimplicialEmbedding( + Mesh& mesh, + const std::vector& label_attributes, + const int64_t& value, + const std::vector& pass_through_attributes) + : m_mesh(mesh) + , m_label_attributes(label_attributes) + , m_value(value) + , m_pass_through_attributes(pass_through_attributes) +{} + +void SimplicialEmbedding::regularize_tags(bool generate_simplicial_embedding) +{ + using namespace operations; + + std::vector todo_handles; + for (size_t i = 0; i < m_label_attributes.size(); ++i) { + attribute::MeshAttributeHandle todo_handle = + m_mesh.register_attribute("todo", m_label_attributes[i].primitive_type(), 1); + todo_handles.emplace_back(todo_handle); + } + + std::deque tag_attributes; + for (size_t i = 0; i < m_label_attributes.size(); ++i) { + tag_attributes.emplace_back( + m_mesh, + m_label_attributes[i], + m_label_attributes[i].primitive_type(), + m_value); + } + + // make sure tag vector is complete and sorted in descending order + for (size_t i = 0; i < tag_attributes.size(); ++i) { + TagAttribute& a = tag_attributes[i]; + if (get_primitive_type_id(a.m_ptype) != m_mesh.top_cell_dimension() - i) { + log_and_throw_error("Tag array must be sorted in descending order starting with " + "the top simplex type up to vertex."); + } + } + + + // tag all faces of attributes + for (size_t attr_it = 0; attr_it < tag_attributes.size() - 1; ++attr_it) { + const TagAttribute& ta = tag_attributes[attr_it]; + for (const Tuple& t : m_mesh.get_all(ta.m_ptype)) { + if (ta.m_accessor.const_scalar_attribute(t) != ta.m_val) { + continue; // t is not tagged + } + + const PrimitiveType face_ptype = + get_primitive_type_from_id(get_primitive_type_id(ta.m_ptype) - 1); + const auto faces = simplex::faces_single_dimension_tuples( + m_mesh, + simplex::Simplex(m_mesh, ta.m_ptype, t), + face_ptype); + + TagAttribute& face_ta = tag_attributes[attr_it + 1]; + for (const Tuple& f : faces) { + face_ta.m_accessor.scalar_attribute(f) = face_ta.m_val; + } + } + } + + if (!generate_simplicial_embedding) { + return; + } + + // check for inversions + if (m_mesh.has_attribute("vertices", PrimitiveType::Vertex)) { + const auto pos_handle = + m_mesh.get_attribute_handle("vertices", PrimitiveType::Vertex); + if (pos_handle.dimension() == m_mesh.top_cell_dimension() - 1) { + SimplexInversionInvariant inv(pos_handle.mesh(), pos_handle.as()); + + if (!inv.after({}, m_mesh.get_all(m_mesh.top_simplex_type()))) { + logger().error("Mesh is not inversion free BEFORE simplicial embedding!"); + } + } + } + + // split untagged simplices that have only tagged faces + for (size_t attr_it = 0; attr_it < tag_attributes.size() - 1;) { + const TagAttribute& ta = tag_attributes[attr_it]; + + attribute::MeshAttributeHandle& todo_handle = todo_handles[attr_it]; + auto todo_acc = m_mesh.create_accessor(todo_handle); + + int64_t count_todos = 0; + for (const Tuple& t : m_mesh.get_all(ta.m_ptype)) { + if (ta.m_accessor.const_scalar_attribute(t) == ta.m_val) { + continue; // t is tagged + } + + const PrimitiveType face_ptype = + get_primitive_type_from_id(get_primitive_type_id(ta.m_ptype) - 1); + const auto faces = simplex::faces_single_dimension_tuples( + m_mesh, + simplex::Simplex(m_mesh, ta.m_ptype, t), + face_ptype); + + const TagAttribute& face_ta = tag_attributes[attr_it + 1]; + + bool all_faces_are_tagged = true; + + for (const Tuple& f : faces) { + if (face_ta.m_accessor.const_scalar_attribute(f) != face_ta.m_val) { + all_faces_are_tagged = false; + break; + } + } + + if (all_faces_are_tagged) { + todo_acc.scalar_attribute(t) = 1; + ++count_todos; + } + } + + // split simplex because all its faces are tagged + Scheduler scheduler; + int64_t count_done = 0; + switch (ta.m_ptype) { + case PrimitiveType::Edge: { // edge split + EdgeSplit op_split(m_mesh); + op_split.add_invariant(std::make_shared( + m_mesh, + std::get>(todo_handle.handle()))); + + // todos + for (const attribute::MeshAttributeHandle& h : todo_handles) { + op_split.set_new_attribute_strategy( + h, + SplitBasicStrategy::None, + SplitRibBasicStrategy::None); + } + // labels + for (const attribute::MeshAttributeHandle& h : m_label_attributes) { + op_split.set_new_attribute_strategy( + h, + SplitBasicStrategy::Copy, + SplitRibBasicStrategy::None); + } + // pass_through + for (const auto& attr : m_pass_through_attributes) { + op_split.set_new_attribute_strategy(attr); + } + + logger().info("split {} edges", count_todos); + while (true) { + const auto stats = scheduler.run_operation_on_all(op_split); + count_done += stats.number_of_successful_operations(); + if (stats.number_of_successful_operations() == 0) { + break; + } + } + + break; + } + case PrimitiveType::Triangle: { // face split + composite::TriFaceSplit op_face_split(m_mesh); + op_face_split.add_invariant(std::make_shared( + m_mesh, + std::get>(todo_handle.handle()))); + + // todos + for (const attribute::MeshAttributeHandle& h : todo_handles) { + op_face_split.split().set_new_attribute_strategy( + h, + SplitBasicStrategy::None, + SplitRibBasicStrategy::None); + op_face_split.collapse().set_new_attribute_strategy(h, CollapseBasicStrategy::None); + } + // labels + for (const attribute::MeshAttributeHandle& h : m_label_attributes) { + op_face_split.split().set_new_attribute_strategy( + h, + SplitBasicStrategy::Copy, + SplitRibBasicStrategy::None); + op_face_split.collapse().set_new_attribute_strategy(h, CollapseBasicStrategy::None); + } + // pass_through + for (const auto& attr : m_pass_through_attributes) { + op_face_split.split().set_new_attribute_strategy(attr); + op_face_split.collapse().set_new_attribute_strategy( + attr, + CollapseBasicStrategy::None); + } + + + logger().info("split {} faces", count_todos); + while (true) { + const auto stats = scheduler.run_operation_on_all(op_face_split); + count_done += stats.number_of_successful_operations(); + if (stats.number_of_successful_operations() == 0) { + break; + } + } + + break; + } + case PrimitiveType::Tetrahedron: { // tet split + composite::TetCellSplit op_tet_split(m_mesh); + op_tet_split.add_invariant(std::make_shared( + m_mesh, + std::get>(todo_handle.handle()))); + + // todos + for (const attribute::MeshAttributeHandle& h : todo_handles) { + op_tet_split.split().set_new_attribute_strategy( + h, + SplitBasicStrategy::None, + SplitRibBasicStrategy::None); + op_tet_split.collapse().set_new_attribute_strategy(h, CollapseBasicStrategy::None); + } + // labels + for (const attribute::MeshAttributeHandle& h : m_label_attributes) { + op_tet_split.split().set_new_attribute_strategy( + h, + SplitBasicStrategy::Copy, + SplitRibBasicStrategy::None); + op_tet_split.collapse().set_new_attribute_strategy(h, CollapseBasicStrategy::None); + } + // pass_through + for (const auto& attr : m_pass_through_attributes) { + op_tet_split.split().set_new_attribute_strategy(attr); + op_tet_split.collapse().set_new_attribute_strategy( + attr, + CollapseBasicStrategy::None); + } + + logger().info("split {} tetrahedra", count_todos); + while (true) { + const auto stats = scheduler.run_operation_on_all(op_tet_split); + count_done += stats.number_of_successful_operations(); + if (stats.number_of_successful_operations() == 0) { + break; + } + } + + break; + } + default: log_and_throw_error("unknown primitive type"); break; + } + + if (count_done == count_todos) { + ++attr_it; + } else { + logger().info( + "{} remain. Regularize same primitive type once more.", + count_todos - count_done); + } + } + + // check for inversions + if (m_mesh.has_attribute("vertices", PrimitiveType::Vertex)) { + const auto pos_handle = + m_mesh.get_attribute_handle("vertices", PrimitiveType::Vertex); + SimplexInversionInvariant inv(pos_handle.mesh(), pos_handle.as()); + + if (pos_handle.dimension() == m_mesh.top_cell_dimension() - 1) { + SimplexInversionInvariant inv(pos_handle.mesh(), pos_handle.as()); + + if (!inv.after({}, m_mesh.get_all(m_mesh.top_simplex_type()))) { + logger().error("Mesh is not inversion free after simplicial embedding!"); + } + } + } +} + + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.hpp b/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.hpp new file mode 100644 index 0000000000..ff334e6a89 --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/internal/SimplicialEmbedding.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include + +namespace wmtk::components::internal { + +/* + * This class is used to seperate mesh and make sure there are no direct connection + * between independent simplicity collection + */ +class SimplicialEmbedding +{ +public: + SimplicialEmbedding( + Mesh& mesh, + const std::vector& label_attributes, + const int64_t& value, + const std::vector& pass_through_attributes); + + /** + * @brief Regularize tags in mesh. + * + * First, it is made sure that tags represent a simplicial complex, i.e., all faces of a tagged + * simplex are also tagged. + * Second (optional), split simplices who's faces are tagged but that are not tagged themselves. + */ + void regularize_tags(bool generate_simplicial_embedding = true); + +private: + Mesh& m_mesh; + + std::vector m_label_attributes; + int64_t m_value; + + std::vector m_pass_through_attributes; +}; + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.cpp b/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.cpp new file mode 100644 index 0000000000..b1dc9538f4 --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.cpp @@ -0,0 +1,35 @@ +#include "simplicial_embedding.hpp" + +#include +#include +#include +#include +#include + +#include "internal/SimplicialEmbedding.hpp" + +namespace wmtk::components { + +void simplicial_embedding(Mesh& mesh, const SimplicialEmbeddingOptions& options) +{ + using namespace internal; + + std::vector tag_attr_vec; + for (const PrimitiveType& ptype : wmtk::utils::primitive_below(mesh.top_simplex_type())) { + tag_attr_vec.emplace_back(options.tag_attributes.at(ptype)); + } + + SimplicialEmbedding rs(mesh, tag_attr_vec, options.value, options.pass_through_attributes); + rs.regularize_tags(options.generate_simplicial_embedding); + + // clean up attributes + { + std::vector keeps = options.pass_through_attributes; + keeps.insert(keeps.end(), tag_attr_vec.begin(), tag_attr_vec.end()); + mesh.clear_attributes(keeps); + } + + mesh.consolidate(); +} + +} // namespace wmtk::components \ No newline at end of file diff --git a/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.hpp b/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.hpp new file mode 100644 index 0000000000..3294d1a20f --- /dev/null +++ b/components/simplicial_embedding/wmtk/components/simplicial_embedding/simplicial_embedding.hpp @@ -0,0 +1,10 @@ +#pragma once +#include + +#include "SimplicialEmbeddingOptions.hpp" + +namespace wmtk::components { + +void simplicial_embedding(Mesh& mesh, const SimplicialEmbeddingOptions& options); + +} // namespace wmtk::components diff --git a/components/tests/simplicial_embedding/test_simplicial_embedding.cpp b/components/tests/simplicial_embedding/test_simplicial_embedding.cpp new file mode 100644 index 0000000000..127af771a1 --- /dev/null +++ b/components/tests/simplicial_embedding/test_simplicial_embedding.cpp @@ -0,0 +1,210 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace wmtk; +using namespace components; + +const std::filesystem::path data_dir = WMTK_DATA_DIR; + +TEST_CASE( + "simplicial_embedding_component_tri", + "[components][simplicial_embedding][trimesh][2D][scheduler]") +{ + tests::DEBUG_TriMesh m = wmtk::tests::hex_plus_two_with_position(); + + SimplicialEmbeddingOptions options; + options.value = 1; + options.tag_attributes[PrimitiveType::Vertex] = + m.register_attribute("vertex_tag", wmtk::PrimitiveType::Vertex, 1); + options.tag_attributes[PrimitiveType::Edge] = + m.register_attribute("edge_tag", wmtk::PrimitiveType::Edge, 1); + options.tag_attributes[PrimitiveType::Triangle] = + m.register_attribute("face_tag", wmtk::PrimitiveType::Triangle, 1); + options.pass_through_attributes.emplace_back( + m.get_attribute_handle("vertices", PrimitiveType::Vertex)); + + SECTION("points_in_2d_case") + { + // 0---1---2 + // / \ / \ / \ . + // 3---4---5---6 + // \ / \ / . + // 7---8 + // set 0 1 4 5 6 + { + const std::vector& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex); + attribute::Accessor acc_vertex_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Vertex]); + acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[1]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[6]) = options.value; + } + + simplicial_embedding(m, options); + + CHECK(m.get_all(PrimitiveType::Vertex).size() == 15); + + if (false) { + wmtk::io::ParaviewWriter writer( + data_dir / "simplicial_embedding_result_0d_case", + "vertices", + m, + true, + true, + true, + false); + m.serialize(writer); + } + } + SECTION("edges_in_2d_case") + { + // 0---1---2 + // / \ / \ / \ . + // 3---4---5---6 + // \ / \ / . + // 7---8 + // set vertex 0 1 4 5 6 7 + // set edge 4-5 5-1 1-4 4-7 7-3 + { + const std::vector& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex); + attribute::Accessor acc_vertex_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Vertex]); + acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[1]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[3]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[4]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[6]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[7]) = options.value; + attribute::Accessor acc_edge_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Edge]); + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(4, 5, 2)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(5, 1, 2)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(1, 4, 2)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(7, 4, 5)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(7, 3, 5)) = options.value; + } + + simplicial_embedding(m, options); + + CHECK(m.get_all(PrimitiveType::Triangle).size() == 17); + CHECK(m.get_all(PrimitiveType::Vertex).size() == 15); + + if (false) { + ParaviewWriter writer( + data_dir / "simplicial_embedding_result_1d_case", + "vertices", + m, + true, + true, + true, + false); + m.serialize(writer); + } + } +} + +TEST_CASE( + "simplicial_embedding_component_tet", + "[components][simplicial_embedding][tetmesh][3D][scheduler][.]") +{ + using namespace tests_3d; + // 0 ---------- 4 + // / \\ // \ . + // / \ \ // \ . + // / \ \ // \ . + // / \ \3 \ . + // 1 --------- 2/ -------- 5 tuple edge 2-3 + // \ / /\ \ / . + // \ / / \\ / . + // \ // \\ / . + // \ // \ / . + // 6 -----------7 + // + DEBUG_TetMesh m = six_cycle_tets(); + + SimplicialEmbeddingOptions options; + options.value = 1; + options.tag_attributes[PrimitiveType::Vertex] = + m.register_attribute("vertex_tag", wmtk::PrimitiveType::Vertex, 1); + options.tag_attributes[PrimitiveType::Edge] = + m.register_attribute("edge_tag", wmtk::PrimitiveType::Edge, 1); + options.tag_attributes[PrimitiveType::Triangle] = + m.register_attribute("face_tag", wmtk::PrimitiveType::Triangle, 1); + options.pass_through_attributes.emplace_back( + m.get_attribute_handle("vertices", PrimitiveType::Vertex)); + + SECTION("points_in_3d_case") + { + { + const std::vector& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex); + attribute::Accessor acc_vertex_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Vertex]); + acc_vertex_tag.scalar_attribute(vertex_tuples[1]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[2]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[3]) = options.value; + } + simplicial_embedding(m, options); + + CHECK(false); // TODO add real checks + + if (false) { + ParaviewWriter writer( + data_dir / "simplicial_embedding_result_points_3d_case", + "vertices", + m, + true, + true, + true, + true); + m.serialize(writer); + } + } + SECTION("edges_in_3d_case") + { + { + const std::vector& vertex_tuples = m.get_all(wmtk::PrimitiveType::Vertex); + attribute::Accessor acc_vertex_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Vertex]); + acc_vertex_tag.scalar_attribute(vertex_tuples[0]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[1]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[2]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[3]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[5]) = options.value; + acc_vertex_tag.scalar_attribute(vertex_tuples[7]) = options.value; + attribute::Accessor acc_edge_tag = + m.create_accessor(options.tag_attributes[PrimitiveType::Edge]); + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(0, 1, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(0, 2, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(0, 3, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(1, 2, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(2, 3, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(3, 1, 0)) = options.value; + acc_edge_tag.scalar_attribute(m.edge_tuple_with_vs_and_t(2, 5, 2)) = options.value; + } + simplicial_embedding(m, options); + + CHECK(false); // TODO add real checks + + if (false) { + ParaviewWriter writer( + data_dir / "simplicial_embedding_result_edges_3d_case", + "vertices", + m, + true, + true, + true, + true); + m.serialize(writer); + } + } +}