diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 2d0b9b3b58..3112f228d5 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -31,6 +31,7 @@ add_application(isotropic_remeshing OFF) add_application(tetwild_msh_converter ON) add_application(tetwild_simplification ON) add_application(triwild ON) +add_application(triwild_submesh ON) add_application(tetwild ON) add_application(cdt_opt ON) add_application(shortest_edge_collapse ON) diff --git a/applications/triwild_submesh/CMakeLists.txt b/applications/triwild_submesh/CMakeLists.txt new file mode 100644 index 0000000000..ffd800d9d8 --- /dev/null +++ b/applications/triwild_submesh/CMakeLists.txt @@ -0,0 +1,19 @@ +wmtk_add_application(triwild_submesh_app + triwild_submesh_main.cpp + triwild_submesh_spec.hpp + triwild_grid.hpp + triwild_grid.cpp + ) + +target_link_libraries(triwild_submesh_app PRIVATE +wmtk::input +wmtk::procedural +wmtk::edge_insertion +wmtk::wildmeshing +wmtk::winding_number +wmtk::output) + +# wmtk_register_integration_test(EXEC_NAME triwild_submesh_app +# CONFIG_FILE ${CMAKE_CURRENT_SOURCE_DIR}/triwild_test_config.json +# GIT_REPOSITORY "https://github.com/wildmeshing/data.git" +# GIT_TAG 82ea9d55c901bbc86d48de39383e9c85d8f67686) diff --git a/applications/triwild_submesh/triwild_grid.cpp b/applications/triwild_submesh/triwild_grid.cpp new file mode 100644 index 0000000000..ccf52c50ca --- /dev/null +++ b/applications/triwild_submesh/triwild_grid.cpp @@ -0,0 +1,66 @@ +#include "triwild_grid.hpp" + +#include +#include + +namespace wmtk::triwild { + +wmtk::TriMesh generate_bg_grid( + const double x_min, + const double y_min, + const double x_max, + const double y_max, + const double length_rel, + const double margin_eps) +{ + const double input_diag = + (Eigen::Vector2d(x_min, y_min) - Eigen::Vector2d(x_max, y_max)).norm(); + const double target_length = length_rel * input_diag; + const double y_start = y_min - input_diag * margin_eps; + const double x_start = x_min - input_diag * margin_eps; + const double y_end = y_max + input_diag * margin_eps; + const double x_end = x_max + input_diag * margin_eps; + + const int64_t x_grid_size = floor((x_end - x_start) / target_length); + const int64_t y_grid_size = floor((y_end - y_start) / target_length); + const double x_space = (x_end - x_start) / x_grid_size; + const double y_space = (y_end - y_start) / y_grid_size; + + Eigen::MatrixX V; + V.resize((x_grid_size + 1) * (y_grid_size + 1), 2); + + for (int64_t i = 0; i < y_grid_size + 1; ++i) { + for (int64_t j = 0; j < x_grid_size + 1; ++j) { + // V(i * (x_grid_size + 1) + j, 0) = x_min + j * x_space; + // V(i * (x_grid_size + 1) + j, 1) = y_min + i * y_space; + V.row(i * (x_grid_size + 1) + j) << x_start + j * x_space, y_start + i * y_space; + } + } + + RowVectors3l F; + F.resize(x_grid_size * y_grid_size * 2, 3); + + for (int64_t i = 0; i < y_grid_size; ++i) { + for (int64_t j = 0; j < x_grid_size; ++j) { + F.row((i * x_grid_size + j) * 2) << i * (x_grid_size + 1) + j, + i * (x_grid_size + 1) + j + 1, (i + 1) * (x_grid_size + 1) + j + 1; + + // std::cout << (i * x_grid_size + j) * 2 << ": " << i * x_grid_size + j << " " + // << i * x_grid_size + j + 1 << " " << (i + 1) * x_grid_size + j + 1 + // << std::endl; + + F.row((i * x_grid_size + j) * 2 + 1) << i * (x_grid_size + 1) + j, + (i + 1) * (x_grid_size + 1) + j + 1, (i + 1) * (x_grid_size + 1) + j; + } + } + + // std::cout << F << std::endl; + + wmtk::TriMesh mesh; + mesh.initialize(F, false); + wmtk::mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, mesh); + + return mesh; +} + +} // namespace wmtk::triwild \ No newline at end of file diff --git a/applications/triwild_submesh/triwild_grid.hpp b/applications/triwild_submesh/triwild_grid.hpp new file mode 100644 index 0000000000..3e2147f940 --- /dev/null +++ b/applications/triwild_submesh/triwild_grid.hpp @@ -0,0 +1,14 @@ +#include +#include + +namespace wmtk::triwild { + +wmtk::TriMesh generate_bg_grid( + const double x_min, + const double y_min, + const double x_max, + const double y_max, + const double length_rel, + const double margin_eps = 0.1); + +} // namespace wmtk::triwild \ No newline at end of file diff --git a/applications/triwild_submesh/triwild_submesh_main.cpp b/applications/triwild_submesh/triwild_submesh_main.cpp new file mode 100644 index 0000000000..d0371db8e5 --- /dev/null +++ b/applications/triwild_submesh/triwild_submesh_main.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "triwild_grid.hpp" +#include "triwild_submesh_spec.hpp" + +using namespace wmtk; +using namespace wmtk::components; +namespace fs = std::filesystem; + +using wmtk::components::utils::resolve_paths; + +int main(int argc, char* argv[]) +{ + opt_logger().set_level(spdlog::level::off); + CLI::App app{argv[0]}; + + app.ignore_case(); + + fs::path json_input_file; + app.add_option("-j, --json", json_input_file, "json specification file") + ->required(true) + ->check(CLI::ExistingFile); + CLI11_PARSE(app, argc, argv); + + nlohmann::json j; + { + std::ifstream ifs(json_input_file); + j = nlohmann::json::parse(ifs); + + jse::JSE spec_engine; + bool r = spec_engine.verify_json(j, triwild_submesh_spec); + if (!r) { + wmtk::logger().error("{}", spec_engine.log2str()); + return 1; + } else { + j = spec_engine.inject_defaults(j, triwild_submesh_spec); + } + } + + nlohmann::json out_json; + + wmtk::utils::StopWatch sw_total("Triwild Total"); + + fs::path input_file = resolve_paths(json_input_file, {j["root"], j["input"]}); + + if (!fs::exists(input_file)) { + log_and_throw_error("File {} does not exist.", input_file.string()); + } + + auto mesh = wmtk::components::input::input(input_file, true); + wmtk::logger().info( + "mesh has {} vertices and {} edges", + mesh->get_all(PrimitiveType::Vertex).size(), + mesh->get_all(PrimitiveType::Edge).size()); + + // get bbox; + double x_min, y_min, x_max, y_max; + x_min = std::numeric_limits::max(); + y_min = std::numeric_limits::max(); + x_max = std::numeric_limits::lowest(); + y_max = std::numeric_limits::lowest(); + + auto mesh_pt_handle = mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); + auto mesh_pt_accessor = mesh->create_const_accessor(mesh_pt_handle); + + for (const auto& v : mesh->get_all(PrimitiveType::Vertex)) { + x_min = std::min(x_min, mesh_pt_accessor.const_vector_attribute(v)[0]); + x_max = std::max(x_max, mesh_pt_accessor.const_vector_attribute(v)[0]); + y_min = std::min(y_min, mesh_pt_accessor.const_vector_attribute(v)[1]); + y_max = std::max(y_max, mesh_pt_accessor.const_vector_attribute(v)[1]); + } + + wmtk::utils::StopWatch sw_bg_mesh("Background Mesh"); + + auto bg_mesh = + wmtk::triwild::generate_bg_grid(x_min, y_min, x_max, y_max, j["target_edge_length"]); + + sw_bg_mesh.stop(); + out_json["runtime_seconds"]["background_mesh"] = sw_bg_mesh.getElapsedTime(); + + if (j["intermediate_output"]) { + wmtk::components::output::output(bg_mesh, "bg_mesh", "vertices"); + } + + wmtk::logger().info("generated bg mesh"); + + wmtk::utils::StopWatch sw_edge_insertion("Edge Insertion"); + + wmtk::components::EdgeInsertionMeshes eim = + wmtk::components::edge_insertion(static_cast(*mesh), bg_mesh); + + sw_edge_insertion.stop(); + out_json["runtime_seconds"]["edge_insertion"] = sw_edge_insertion.getElapsedTime(); + + wmtk::logger().info("finised edge insertion"); + + auto trimesh = eim.tri_mesh; + auto edgemesh = eim.inserted_edge_mesh; + auto bboxmesh = eim.bbox_mesh; + + std::string output_file = j["output"]; + + if (j["intermediate_output"]) { + wmtk::components::output::output(*trimesh, output_file + "_after_insertion", "vertices"); + wmtk::components::output::output( + *edgemesh, + output_file + "_after_insertion_edge_mesh", + "vertices"); + } + + // clean up + { + auto pos_handle = + trimesh->get_attribute_handle("vertices", PrimitiveType::Vertex); + trimesh->clear_attributes({pos_handle}); + } + + std::vector enves; + { + wmtk::components::EnvelopeOptions e; + e.envelope_name = "input"; + e.envelope_constrained_mesh = edgemesh; + e.envelope_geometry_mesh = edgemesh; + e.constrained_position_name = "vertices"; + e.geometry_position_name = "vertices"; + e.thickness = j["envelope_size"]; + + // if (e.envelope_name == "input") { + // e.envelope_geometry_mesh = mesh; // set as input + //} + + enves.push_back(e); + + wmtk::components::EnvelopeOptions e2; + e2.envelope_name = "bbox"; + e2.envelope_constrained_mesh = bboxmesh; + e2.envelope_geometry_mesh = bboxmesh; + e2.constrained_position_name = "vertices"; + e2.geometry_position_name = "vertices"; + e2.thickness = 0.0001; + + enves.push_back(e2); + } + + + wmtk::components::WildMeshingOptions wmo; + wmo.input_mesh = trimesh; + wmo.input_mesh_position = "vertices"; + wmo.target_edge_length = j["target_edge_length"]; + wmo.target_max_amips = j["target_max_amips"]; + wmo.max_passes = j["max_passes"]; + wmo.replace_double_coordinate = false; + wmo.scheduler_update_frequency = 0; + wmo.intermediate_output = j["intermediate_output"]; + wmo.intermediate_output_path = ""; + wmo.intermediate_output_name = j["output"]; + wmo.envelopes = enves; + // wmo.pass_through = pass_through; + wmo.skip_split = j["skip_split"]; + wmo.skip_collapse = j["skip_collapse"]; + wmo.skip_swap = j["skip_swap"]; + wmo.skip_smooth = j["skip_smooth"]; + wmo.use_embedding = true; + + + wmtk::utils::StopWatch sw_wildmeshing("Wildmeshing"); + + auto meshes_after_tetwild = wildmeshing(wmo); + + sw_wildmeshing.stop(); + out_json["runtime_seconds"]["wildmeshing"] = sw_wildmeshing.getElapsedTime(); + + assert(meshes_after_tetwild.size() == 1); + auto main_mesh = meshes_after_tetwild[0].first; + + wmtk::components::output::output(*main_mesh, output_file, "vertices"); + + sw_total.stop(); + out_json["runtime_seconds"]["total"] = sw_total.getElapsedTime(); + + const std::string report = j["report"]; + if (!report.empty()) { + out_json["vertices"] = main_mesh->get_all(PrimitiveType::Vertex).size(); + out_json["edges"] = main_mesh->get_all(PrimitiveType::Edge).size(); + out_json["cells"] = main_mesh->get_all(PrimitiveType::Triangle).size(); + + out_json["input"] = j; + + std::ofstream ofs(report); + ofs << std::setw(4) << out_json; + } + + return 0; +} \ No newline at end of file diff --git a/applications/triwild_submesh/triwild_submesh_spec.hpp b/applications/triwild_submesh/triwild_submesh_spec.hpp new file mode 100644 index 0000000000..8d4659062a --- /dev/null +++ b/applications/triwild_submesh/triwild_submesh_spec.hpp @@ -0,0 +1,96 @@ +#pragma once +#include +namespace { + +nlohmann::json triwild_submesh_spec = R"( +[ + { + "pointer": "/", + "type": "object", + "required": ["input", "output"], + "optional": [ + "root", + "report", + "envelope_size", + "target_edge_length", + "target_max_amips", + "max_passes", + "intermediate_output", + "skip_split", + "skip_collapse", + "skip_swap", + "skip_smooth" + ] +}, +{ + "pointer": "/input", + "type": "string" +}, +{ + "pointer": "/envelope_size", + "type": "float", + "default": 0.001 +}, +{ + "pointer": "/length_rel", + "type": "float", + "default": 0.1 +}, +{ + "pointer": "/intermediate_output", + "type": "bool", + "default": false +}, +{ + "pointer": "/skip_split", + "type": "bool", + "default": false +}, +{ + "pointer": "/skip_collapse", + "type": "bool", + "default": false +}, +{ + "pointer": "/skip_swap", + "type": "bool", + "default": false +}, +{ + "pointer": "/skip_smooth", + "type": "bool", + "default": false +}, +{ + "pointer": "/target_max_amips", + "type": "float", + "default": 10.0 +}, +{ + "pointer": "/target_edge_length", + "type": "float", + "default": 0.05 +}, +{ + "pointer": "/max_passes", + "type": "int", + "default": 10 +}, +{ + "pointer": "/output", + "type": "string" +}, +{ + "pointer": "/report", + "type": "string", + "default": "" +}, +{ + "pointer": "/root", + "type": "string", + "default": "" +} +] +)"_json; + +} \ No newline at end of file diff --git a/applications/triwild_submesh/triwild_submesh_test_config.json b/applications/triwild_submesh/triwild_submesh_test_config.json new file mode 100644 index 0000000000..10d836ab5a --- /dev/null +++ b/applications/triwild_submesh/triwild_submesh_test_config.json @@ -0,0 +1,9 @@ +{ + "test_directory": "unit_test", + "tests": ["triwild_siggraph_icon_out.json"], + "input_tag": "input", + "oracle_tag": "report", + "input_directory_tag": "root", + "platform": "Linux", + "checks": [] +} \ No newline at end of file diff --git a/components/tests/wildmeshing/test_wildmeshing.cpp b/components/tests/wildmeshing/test_wildmeshing.cpp new file mode 100644 index 0000000000..39e9ed584c --- /dev/null +++ b/components/tests/wildmeshing/test_wildmeshing.cpp @@ -0,0 +1,125 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +using namespace wmtk; +using namespace components; + +constexpr PrimitiveType PV = PrimitiveType::Vertex; +constexpr PrimitiveType PE = PrimitiveType::Edge; +constexpr PrimitiveType PF = PrimitiveType::Triangle; + +const std::filesystem::path data_dir = WMTK_DATA_DIR; + +TEST_CASE("test_wildmeshing_submesh", "[wildmeshing][.]") +{ + // logger().set_level(spdlog::level::off); + logger().set_level(spdlog::level::debug); + opt_logger().set_level(spdlog::level::warn); + + std::shared_ptr mesh_ptr = + std::make_shared(tests::edge_region_with_position_2D()); + tests::DEBUG_TriMesh& m = *mesh_ptr; + + // add input + std::shared_ptr child_input_ptr; + { + child_input_ptr = std::make_shared(tests::single_line()); + + const Tuple child_tuple = child_input_ptr->edge_tuple_from_vids(0, 1); + const Tuple parent_tuple = m.edge_tuple_from_vids(4, 5); + std::vector> map_tuples; + map_tuples.emplace_back(std::array{child_tuple, parent_tuple}); + m.register_child_mesh(child_input_ptr, map_tuples); + } + // add boundary + std::shared_ptr child_bnd_ptr; + { + child_bnd_ptr = std::make_shared(); + tests::DEBUG_EdgeMesh& c = *child_bnd_ptr; + + RowVectors2l edges; + edges.resize(8, 2); + edges.row(0) << 0, 1; + edges.row(1) << 1, 2; + edges.row(2) << 2, 3; + edges.row(3) << 3, 4; + edges.row(4) << 4, 5; + edges.row(5) << 5, 6; + edges.row(6) << 6, 7; + edges.row(7) << 7, 0; + c.initialize(edges); + + std::vector> map_tuples; + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(0, 1), m.edge_tuple_from_vids(0, 1)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(1, 2), m.edge_tuple_from_vids(1, 2)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(2, 3), m.edge_tuple_from_vids(2, 6)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(3, 4), m.edge_tuple_from_vids(6, 9)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(4, 5), m.edge_tuple_from_vids(9, 8)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(5, 6), m.edge_tuple_from_vids(8, 7)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(6, 7), m.edge_tuple_from_vids(7, 3)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(7, 0), m.edge_tuple_from_vids(3, 0)}); + m.register_child_mesh(child_bnd_ptr, map_tuples); + } + + // envelopes + std::vector enves; + { + wmtk::components::EnvelopeOptions e; + e.envelope_name = "input"; + e.envelope_constrained_mesh = child_input_ptr; + e.envelope_geometry_mesh = child_input_ptr; + e.thickness = 1e-3; + + enves.push_back(e); + + wmtk::components::EnvelopeOptions e2; + e2.envelope_name = "bbox"; + e2.envelope_constrained_mesh = child_bnd_ptr; + e2.envelope_geometry_mesh = child_bnd_ptr; + e2.thickness = 1e-4; + + enves.push_back(e2); + } + + wmtk::components::WildMeshingOptions wmo; + wmo.input_mesh = mesh_ptr; + wmo.target_edge_length = 0.05; + wmo.max_passes = 5; + wmo.replace_double_coordinate = true; + wmo.intermediate_output_path = ""; + wmo.intermediate_output_name = "out"; + wmo.intermediate_output = true; + wmo.envelopes = enves; + wmo.use_embedding = true; + + auto meshes_after_tetwild = wildmeshing(wmo); + + // CHECK(emb.has_child_mesh()); + // CHECK(emb.get_child_meshes().size() == 2); + //{ + // using submesh::utils::write; + // CHECK_NOTHROW(write("wildmeshing_submesh", "vertices", emb, true, true, true, false)); + // } +} \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/CMakeLists.txt b/components/wildmeshing/wmtk/components/wildmeshing/CMakeLists.txt index 150297afae..0d52a73ffd 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/CMakeLists.txt +++ b/components/wildmeshing/wmtk/components/wildmeshing/CMakeLists.txt @@ -14,6 +14,14 @@ set(SRC_FILES internal/wildmeshing3d.cpp internal/wildmeshing_utils.hpp internal/wildmeshing_utils.cpp - internal/WildmeshingOptions.hpp) + internal/WildmeshingOptions.hpp + internal/wildmeshing_embedding_2d.hpp + internal/wildmeshing_embedding_2d.cpp + utils/IntermediateWrite.hpp + utils/IntermediateWrite.cpp +) target_sources(wmtk_${COMPONENT_NAME} PRIVATE ${SRC_FILES}) + +add_component_test(wmtk::${COMPONENT_NAME} + ${PROJECT_SOURCE_DIR}/components/tests/wildmeshing/test_wildmeshing.cpp) \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/WildmeshingOptions.hpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/WildmeshingOptions.hpp index 539667f031..a906a39db0 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/internal/WildmeshingOptions.hpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/WildmeshingOptions.hpp @@ -5,35 +5,91 @@ namespace wmtk::components { struct EnvelopeOptions { + /** + * Name of the envelope. Mostly used for debugging. + */ std::string envelope_name; + /** + * The mesh that is used to construct the envelope. + */ std::shared_ptr envelope_constrained_mesh; + /** + * The mesh that is checked if it is within the envelope. + */ std::shared_ptr envelope_geometry_mesh; - std::string constrained_position_name; - std::string geometry_position_name; - double thickness; + /** + * The position attribute name of the envelope-construction mesh. + */ + std::string constrained_position_name = "vertices"; + /** + * The position attribute name for the mesh that must remain within the envelope. + */ + std::string geometry_position_name = "vertices"; + /** + * The envelope thickness, relative to the bounding box of the embedding. + */ + double thickness = 1e-6; }; struct WildMeshingOptions { + /** + * The mesh that is embedding the submeshes (represented by the envelope options). + */ std::shared_ptr input_mesh; - std::string input_mesh_position; - double target_edge_length; - double target_max_amips; - double max_passes; - bool intermediate_output; - bool replace_double_coordinate; - size_t scheduler_update_frequency; - std::string intermediate_output_path; + /** + * The position attribute of the input mesh, i.e., the embedding. + */ + std::string input_mesh_position = "vertices"; + /** + * Target edge length, relative to the bounding box of the input mesh. + */ + double target_edge_length = 0.05; + /** + * The targeted max amips. The optimization will continue until every cell has an amips below + * the max or until it runs out of iterations. + */ + double target_max_amips = 10.0; + /** + * The maximum number of iterations of the optimization. + */ + int64_t max_passes = 10; + /** + * Debugging output. + */ + bool intermediate_output = false; + /** + * Replace the `double` vertex positions by `Rational` positions. Use this if your mesh has its + * positions stored as `double`. + */ + bool replace_double_coordinate = false; + /** + * TODO document (I don't know what this is good for) + */ + size_t scheduler_update_frequency = 0; + std::string intermediate_output_path = ""; std::string intermediate_output_name; - bool skip_split; - bool skip_collapse; - bool skip_swap; - bool skip_smooth; + bool skip_split = false; + bool skip_collapse = false; + bool skip_swap = false; + bool skip_smooth = false; + /** + * The envelopes ensure that substructures remain in their position up to a given error. + */ std::vector envelopes; + /** + * Further attributes that will be handled with default behavior. + */ std::vector pass_through; + + /** + * Convert the multimesh into tags und use Embedding/SubMesh. + * For now, no topologica checks on the submesh will be performed. + */ + bool use_embedding = false; }; } // namespace wmtk::components \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.cpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.cpp new file mode 100644 index 0000000000..2c83bfbb8e --- /dev/null +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.cpp @@ -0,0 +1,880 @@ +#include "wildmeshing2d.hpp" + +#include +#include "WildmeshingOptions.hpp" + +#include +#include +#include +#include + +#include +#include + + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +namespace wmtk::components::internal { + +using namespace simplex; +using namespace operations; +using namespace operations::tri_mesh; +using namespace operations::tet_mesh; +using namespace operations::composite; +using namespace function; +using namespace invariants; + +double compute_max_amips( + const attribute::MeshAttributeHandle& amips_handle, + double& last_max_amips, + const bool verbose = true) +{ + const Mesh& mesh = amips_handle.mesh(); + const auto amips_accessor = mesh.create_const_accessor(amips_handle); + + // compute max energy + double max_energy = std::numeric_limits::lowest(); + double min_energy = std::numeric_limits::max(); + double avg_energy = 0; + for (const Tuple& t : mesh.get_all(mesh.top_simplex_type())) { + double e = amips_accessor.const_scalar_attribute(t); + max_energy = std::max(max_energy, e); + min_energy = std::min(min_energy, e); + avg_energy += e; + } + + avg_energy = avg_energy / mesh.get_all(mesh.top_simplex_type()).size(); + + if (verbose) { + logger().info( + "Max AMIPS Energy: {:>6.2f}, Min AMIPS Energy: {:>6.2f}, Avg AMIPS Energy: {:>6.2f}", + max_energy, + min_energy, + avg_energy); + + if (max_energy > last_max_amips) { + logger().info("Max amips increased from {} to {}", last_max_amips, max_energy); + } + } + + last_max_amips = max_energy; + + return max_energy; +} + +void print_stats(const wmtk::SchedulerStats& stats, const std::string& name = "") +{ + if (name.empty()) { + logger().info( + "Executed {} ops (S/F) {}/{}. Time executing: {:.4f}", + stats.number_of_performed_operations(), + stats.number_of_successful_operations(), + stats.number_of_failed_operations(), + stats.executing_time); + } else { + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time executing: {:.4f}", + name, + stats.number_of_performed_operations(), + stats.number_of_successful_operations(), + stats.number_of_failed_operations(), + stats.executing_time); + } +} + +int64_t count_unrounded_vertices( + const attribute::MeshAttributeHandle& pt_handle, + int64_t& last_num_unrounded, + const bool verbose = true) +{ + const Mesh& mesh = pt_handle.mesh(); + const auto pt_accessor = mesh.create_const_accessor(pt_handle); + + // verbose logger, can be removed + int64_t unrounded = 0; + for (const auto& v : mesh.get_all(PrimitiveType::Vertex)) { + const auto p = pt_accessor.const_vector_attribute(v); + for (int64_t d = 0; d < 2; ++d) { + if (!p[d].is_rounded()) { + ++unrounded; + break; + } + } + } + + if (unrounded > 0 && verbose) { + logger().info("Mesh has {} unrounded vertices", unrounded); + } + + if (unrounded > last_num_unrounded) { + logger().error( + "Number of unrounded vertices increased from {} to {}", + last_num_unrounded, + unrounded); + } + last_num_unrounded = unrounded; + + return unrounded; +} + +void test_flipped_triangles(const Mesh& mesh, const attribute::MeshAttributeHandle& pt_handle) +{ + const auto pt_accessor = mesh.create_const_accessor(pt_handle); + + for (const Tuple& t : mesh.get_all(mesh.top_simplex_type())) { + const auto vertices = mesh.orient_vertices(t); + std::vector pos; + for (int i = 0; i < vertices.size(); ++i) { + pos.push_back(pt_accessor.const_vector_attribute(vertices[i])); + } + if (wmtk::utils::wmtk_orient2d(pos[0], pos[1], pos[2]) <= 0) { + logger().error("Flipped triangle!"); + } + } +} + +std::vector, std::string>> wildmeshing_embedding_2d( + const WildMeshingOptions& options) +{ + auto& mesh = *options.input_mesh; + + if (!mesh.is_connectivity_valid()) { + log_and_throw_error("input mesh for wildmeshing connectivity invalid"); + } + + logger().trace("Getting rational point handle"); + + ////////////////////////////////// + // Retriving vertices + // + if (options.replace_double_coordinate) { + logger().trace("Found double attribute"); + auto pt_double_attribute = + mesh.get_attribute_handle(options.input_mesh_position, PrimitiveType::Vertex); + + if (!mesh.has_attribute(options.input_mesh_position, PrimitiveType::Vertex)) { + wmtk::utils::cast_attribute( + pt_double_attribute, + mesh, + options.input_mesh_position); + + + } else { + auto pt_attribute = mesh.get_attribute_handle( + options.input_mesh_position, + PrimitiveType::Vertex); + wmtk::utils::cast_attribute(pt_double_attribute, pt_attribute); + } + mesh.delete_attribute(pt_double_attribute); + } + auto pt_attribute = + mesh.get_attribute_handle(options.input_mesh_position, PrimitiveType::Vertex); + + ////////////////////////////////// + // computing bbox diagonal + const double bbdiag = wmtk::utils::bbox_diagonal_from_mesh(pt_attribute); + logger().info("bbox diag = {}", bbdiag); + + const double target_edge_length = options.target_edge_length * bbdiag; + logger().info("target edge length: {}", target_edge_length); + + ////////////////////////////////// + // substructures + ////////////////////////////////// + std::map, std::shared_ptr> mm_to_submesh_map; + auto emb_ptr = submesh::utils::submesh_from_multimesh(options.input_mesh, mm_to_submesh_map); + submesh::Embedding& emb = *emb_ptr; + // deregister all child meshes + for (const auto& [child, sub] : mm_to_submesh_map) { + mesh.deregister_child_mesh(child); + } + emb.update_tag_attribute_handles(); + + ////////////////////////////////// + // store amips + auto amips_attribute = + mesh.register_attribute("wildmeshing_amips", mesh.top_simplex_type(), 1); + auto amips_accessor = mesh.create_accessor(amips_attribute.as()); + // amips update + auto compute_amips = [](const Eigen::MatrixX& P) -> Eigen::VectorXd { + assert(P.rows() == 2); // rows --> attribute dimension + assert(P.cols() == 3); + // triangle + std::array pts; + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 2; ++j) { + pts[2 * i + j] = P(j, i).to_double(); + } + } + const double a = Tri_AMIPS_energy(pts); + return Eigen::VectorXd::Constant(1, a); + }; + auto amips_update = + std::make_shared>( + amips_attribute, + pt_attribute, + compute_amips); + amips_update->run_on_all(); + + double last_max_amips = std::numeric_limits::max(); + compute_max_amips(amips_attribute, last_max_amips); + + ////////////////////////////////// + // Storing target edge length + auto target_edge_length_attribute = mesh.register_attribute( + "wildmeshing_target_edge_length", + PrimitiveType::Edge, + 1, + false, + target_edge_length); // defaults to target edge length + + // Target edge length update + const double min_edge_length = [&]() -> double { + if (options.envelopes.empty()) { + return 1e-6; // some default value if no envelope exists + } else { + // use envelope thickness if available + double r = 0; + for (const EnvelopeOptions& e : options.envelopes) { + r = std::max(r, e.thickness * bbdiag); + } + assert(r > 0); + return r; + } + }(); + const double& target_max_amips = options.target_max_amips; + + + ////////////////////////////////// + // Storing edge lengths + auto edge_length_attribute = + mesh.register_attribute("edge_length", PrimitiveType::Edge, 1); + auto edge_length_accessor = mesh.create_accessor(edge_length_attribute.as()); + // Edge length update + auto compute_edge_length = [](const Eigen::MatrixX& P) -> Eigen::VectorXd { + assert(P.cols() == 2); + assert(P.rows() == 2 || P.rows() == 3); + return Eigen::VectorXd::Constant(1, sqrt((P.col(0) - P.col(1)).squaredNorm().to_double())); + }; + auto edge_length_update = + std::make_shared>( + edge_length_attribute, + pt_attribute, + compute_edge_length); + edge_length_update->run_on_all(); + + ////////////////////////////////// + // compute frozen vertices + ////////////////////////////////// + auto frozen_vertex_attribute = + mesh.register_attribute("frozen_vertex", PrimitiveType::Vertex, 1); + auto frozen_vertex_accessor = mesh.create_accessor(frozen_vertex_attribute.as()); + + for (const auto& sub_ptr : emb.get_child_meshes()) { + const submesh::SubMesh& sub = *sub_ptr; + for (const simplex::IdSimplex& v : sub.get_all_id_simplex(PrimitiveType::Vertex)) { + int64_t counter = 0; + for (const Tuple& cof : cofaces_single_dimension_iterable( + mesh, + mesh.get_simplex(v), + sub.top_simplex_type())) { + if (sub.contains(cof, sub.top_simplex_type())) { + ++counter; + } + } + if (counter < 2) { + frozen_vertex_accessor.scalar_attribute(v) = 1; + } + } + } + + // auto input_ptr = mesh.get_child_meshes().front(); + // + // for (const Tuple& v : input_ptr->get_all(PrimitiveType::Vertex)) { + // if (input_ptr->is_boundary(PrimitiveType::Vertex, v)) { + // const auto& parent_v = + // input_ptr->map_to_parent(simplex::Simplex::vertex(*input_ptr, v)); + // frozen_vertex_accessor.scalar_attribute(parent_v) = 1; + // } + // } + + { + int64_t frozen_cnt = 0; + for (const auto& v : mesh.get_all(PrimitiveType::Vertex)) { + if (frozen_vertex_accessor.scalar_attribute(v) == 1) { + frozen_cnt++; + } + } + logger().info("mesh has {} frozen vertices", frozen_cnt); + } + + ////////////////////////////////// + // default transfer + auto pass_through_attributes = options.pass_through; + pass_through_attributes.push_back(edge_length_attribute); + pass_through_attributes.push_back(amips_attribute); + + + ////////////////////////////////// + // Lambdas for priority + ////////////////////////////////// + auto long_edges_first = [&](const simplex::Simplex& s) { + assert(s.primitive_type() == PrimitiveType::Edge); + return -edge_length_accessor.scalar_attribute(s.tuple()); + }; + auto short_edges_first = [&](const simplex::Simplex& s) { + assert(s.primitive_type() == PrimitiveType::Edge); + return edge_length_accessor.scalar_attribute(s.tuple()); + }; + + + ////////////////////////////////// + // envelopes + ////////////////////////////////// + + using MeshConstrainPair = ProjectOperation::MeshConstrainPair; + + auto envelope_invariant = std::make_shared(mesh); + + for (const EnvelopeOptions& e : options.envelopes) { + Mesh& constrained_mesh = *e.envelope_constrained_mesh; + Mesh& geometry_mesh = *e.envelope_geometry_mesh; + + logger().info("wildmeshing2d: registered {} mesh as envelope constraints", e.envelope_name); + + if (&constrained_mesh != &geometry_mesh) { + log_and_throw_error( + "wildmeshing2d: constrained and geometry mesh must be the same for submesh usage."); + } + + submesh::SubMesh& sub = *mm_to_submesh_map.at(e.envelope_geometry_mesh); + + envelope_invariant->add( + std::make_shared(pt_attribute, e.thickness * bbdiag, sub)); + } + + + ////////////////////////////////// + // Invariants + ////////////////////////////////// + + logger().trace("Going through invariants"); + auto inversion_invariant = + std::make_shared>(mesh, pt_attribute.as()); + + std::shared_ptr amips = + std::make_shared(mesh, pt_attribute); + // auto function_invariant = std::make_shared(mesh.top_simplex_type(), + // amips); + auto function_invariant = + std::make_shared(mesh.top_simplex_type(), amips); + + auto link_condition = std::make_shared(mesh); + + auto todo_larger = std::make_shared( + mesh, + edge_length_attribute.as(), + target_edge_length_attribute.as(), + 4.0 / 3.0); + + auto todo_smaller = std::make_shared( + mesh, + edge_length_attribute.as(), + target_edge_length_attribute.as(), + 4.0 / 5.0); + + + auto interior_edge = std::make_shared(mesh); + + // for (const auto& em : multimesh_meshes) { + // interior_edge->add_boundary(*(em.first)); + // } + + // auto invariant_separate_substructures = + // std::make_shared( + // mesh); // TODO remove for submesh + + auto frozen_vertex_invariant = std::make_shared( + mesh, + frozen_vertex_attribute.as()); + auto frozen_opp_vertex_invariant = std::make_shared( + mesh, + frozen_vertex_attribute.as()); + + ////////////////////////////////// + // renew flags + ////////////////////////////////// + auto visited_edge_flag = + mesh.register_attribute("visited_edge", PrimitiveType::Edge, 1, false, char(1)); + + auto update_flag_func = [](const Eigen::MatrixX& P) -> Eigen::VectorX { + assert(P.cols() == 2); + assert(P.rows() == 2 || P.rows() == 3); + return Eigen::VectorX::Constant(1, char(1)); + }; + auto tag_update = + std::make_shared>( + visited_edge_flag, + pt_attribute, + update_flag_func); + + + ////////////////////////////////// + // Creation of the 4 ops + ////////////////////////////////// + std::vector> ops; + std::vector ops_name; + + ////////////////////////////////// + // 0) Rounding + ////////////////////////////////// + auto rounding_pt_attribute = mesh.get_attribute_handle_typed( + options.input_mesh_position, + PrimitiveType::Vertex); + auto rounding = std::make_shared(mesh, rounding_pt_attribute); + rounding->add_invariant( + std::make_shared(mesh, pt_attribute.as(), true)); + rounding->add_invariant(inversion_invariant); + + + ////////////////////////////////// + // 1) EdgeSplit + ////////////////////////////////// + auto split = std::make_shared(emb); + split->set_priority(long_edges_first); + + split->add_invariant(todo_larger); + split->add_invariant(inversion_invariant); + + split->set_new_attribute_strategy(pt_attribute); + split->set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + + split->set_new_attribute_strategy( + frozen_vertex_attribute, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + split->set_new_attribute_strategy( + target_edge_length_attribute, + wmtk::operations::SplitBasicStrategy::Copy, + wmtk::operations::SplitRibBasicStrategy::Mean); + + for (const auto& attr : pass_through_attributes) { + split->set_new_attribute_strategy( + attr, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + } + + auto split_then_round = std::make_shared(mesh); + split_then_round->add_operation(split); + split_then_round->add_operation(rounding); + + // split unrounded + auto split_unrounded = std::make_shared(emb); + split_unrounded->set_priority(long_edges_first); + + split_unrounded->add_invariant(todo_larger); + + auto split_unrounded_transfer_strategy = + std::make_shared>(pt_attribute); + split_unrounded_transfer_strategy->set_strategy( + [](const Eigen::VectorX& a, const std::bitset<2>&) { + return std::array, 2>{{a, a}}; + }); + split_unrounded_transfer_strategy->set_rib_strategy( + [](const Eigen::VectorX& p0_d, + const Eigen::VectorX& p1_d, + const std::bitset<2>& bs) -> Eigen::VectorX { + Eigen::VectorX p0(p0_d.size()); + Eigen::VectorX p1(p1_d.size()); + for (int i = 0; i < p0_d.size(); ++i) { + p0[i] = Rational(p0_d[i], false); + p1[i] = Rational(p1_d[i], false); + } + if (bs[0] == bs[1]) { + return (p0 + p1) / Rational(2, false); + } else if (bs[0]) { + return p0; + + } else { + return p1; + } + }); + + split_unrounded->set_new_attribute_strategy(pt_attribute, split_unrounded_transfer_strategy); + split_unrounded->set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + split_unrounded->set_new_attribute_strategy( + frozen_vertex_attribute, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + for (const auto& attr : pass_through_attributes) { + split_unrounded->set_new_attribute_strategy( + attr, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + } + split_unrounded->set_new_attribute_strategy( + target_edge_length_attribute, + wmtk::operations::SplitBasicStrategy::Copy, + wmtk::operations::SplitRibBasicStrategy::Mean); + + auto split_sequence = std::make_shared(mesh); + split_sequence->add_operation(split_then_round); + split_sequence->add_operation(split_unrounded); + + split_sequence->set_priority(long_edges_first); + + split_sequence->add_transfer_strategy(amips_update); + split_sequence->add_transfer_strategy(edge_length_update); + split_sequence->add_transfer_strategy(tag_update); // for renew the queue + + + if (!options.skip_split) { + ops.emplace_back(split_sequence); + ops_name.emplace_back("SPLIT"); + + ops.emplace_back(rounding); // TODO split tries to round already + ops_name.emplace_back("rounding"); + } + + + ////////////////////////////////// + // collapse transfer + ////////////////////////////////// + auto clps_strat1 = std::make_shared>(pt_attribute); + clps_strat1->set_strategy(CollapseBasicStrategy::CopyOther); + + auto clps_strat2 = std::make_shared>(pt_attribute); + clps_strat2->set_strategy(CollapseBasicStrategy::CopyTuple); + + ////////////////////////////////// + // 2) EdgeCollapse + ////////////////////////////////// + + auto setup_collapse = [&](std::shared_ptr& collapse) { + // collapse->add_invariant(invariant_separate_substructures); + // collapse->add_invariant(std::make_shared(mesh)); + collapse->add_invariant(link_condition); + collapse->add_invariant(inversion_invariant); + collapse->add_invariant(function_invariant); + collapse->add_invariant(envelope_invariant); + + collapse->set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::CollapseBasicStrategy::None); + + for (const auto& attr : pass_through_attributes) { + collapse->set_new_attribute_strategy( + attr, + wmtk::operations::CollapseBasicStrategy::None); + } + collapse->set_new_attribute_strategy( + target_edge_length_attribute, + wmtk::operations::CollapseBasicStrategy::None); + + collapse->add_transfer_strategy(tag_update); + collapse->add_transfer_strategy(amips_update); + collapse->add_transfer_strategy(edge_length_update); + }; + + auto collapse1 = std::make_shared(emb); + + collapse1->add_invariant(frozen_vertex_invariant); + collapse1->set_new_attribute_strategy(pt_attribute, clps_strat1); + collapse1->set_new_attribute_strategy( + frozen_vertex_attribute, + CollapseBasicStrategy::CopyOther); + setup_collapse(collapse1); + + auto collapse2 = std::make_shared(emb); + + collapse2->add_invariant(frozen_opp_vertex_invariant); + collapse2->set_new_attribute_strategy(pt_attribute, clps_strat2); + collapse2->set_new_attribute_strategy( + frozen_vertex_attribute, + CollapseBasicStrategy::CopyTuple); + setup_collapse(collapse2); + + auto collapse = std::make_shared(mesh); + collapse->add_operation(collapse1); + collapse->add_operation(collapse2); + collapse->add_invariant(todo_smaller); + + auto collapse_then_round = std::make_shared(mesh); + collapse_then_round->add_operation(collapse); + collapse_then_round->add_operation(rounding); + + collapse_then_round->set_priority(short_edges_first); + + if (!options.skip_collapse) { + ops.emplace_back(collapse_then_round); + ops_name.emplace_back("COLLAPSE"); + + ops.emplace_back(rounding); // TODO collapse rounds already, why a second rounding? + ops_name.emplace_back("rounding"); + } + + ////////////////////////////////// + // 3) Swap + ////////////////////////////////// + + auto setup_swap = [&](Operation& op, + EdgeCollapse& collapse, + EdgeSplit& split, + std::shared_ptr simplex_invariant, + bool is_edge = true) { + if (is_edge) op.set_priority(long_edges_first); + + op.add_invariant(simplex_invariant); + op.add_invariant( + std::make_shared(mesh, pt_attribute.as())); + op.add_invariant(inversion_invariant); + op.add_invariant(function_invariant); + + op.add_transfer_strategy(amips_update); + op.add_transfer_strategy(edge_length_update); + op.add_transfer_strategy(tag_update); + + collapse.add_invariant(link_condition); + collapse.add_invariant(envelope_invariant); + + + collapse.set_new_attribute_strategy(pt_attribute, CollapseBasicStrategy::CopyOther); + split.set_new_attribute_strategy(pt_attribute); + + split.set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + + collapse.set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::CollapseBasicStrategy::None); + + split.set_new_attribute_strategy( + frozen_vertex_attribute, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + + collapse.set_new_attribute_strategy( + frozen_vertex_attribute, + CollapseBasicStrategy::CopyOther); + + split.set_new_attribute_strategy( + target_edge_length_attribute, + wmtk::operations::SplitBasicStrategy::Copy, + wmtk::operations::SplitRibBasicStrategy::Mean); + collapse.set_new_attribute_strategy( + target_edge_length_attribute, + wmtk::operations::CollapseBasicStrategy::None); + + for (const auto& attr : pass_through_attributes) { + split.set_new_attribute_strategy( + attr, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + collapse.set_new_attribute_strategy( + attr, + wmtk::operations::CollapseBasicStrategy::None); + } + }; + + + auto swap = std::make_shared(emb); + setup_swap(*swap, swap->collapse(), swap->split(), interior_edge); + + if (!options.skip_swap) { + ops.push_back(swap); + ops_name.push_back("swap"); + + ops.emplace_back(rounding); + ops_name.emplace_back("rounding"); + } + + ///////////////////////////////////////// + // 4) Smoothing + ///////////////////////////////////////// + auto smoothing = std::make_shared(mesh, pt_attribute); + smoothing->add_invariant(std::make_shared(mesh, pt_attribute.as())); + smoothing->add_invariant(frozen_vertex_invariant); + smoothing->add_invariant(inversion_invariant); + + auto proj_smoothing = std::make_shared(smoothing, emb, pt_attribute); + // proj_smoothing->use_random_priority() = true; + proj_smoothing->add_invariant(frozen_vertex_invariant); + proj_smoothing->add_invariant(envelope_invariant); + proj_smoothing->add_invariant(inversion_invariant); + //proj_smoothing->add_invariant(function_invariant); // TODO confirm that this should be here + + proj_smoothing->add_transfer_strategy(amips_update); + proj_smoothing->add_transfer_strategy(edge_length_update); + + + if (!options.skip_smooth) { + for (int i = 0; i < 1; ++i) { + ops.push_back(proj_smoothing); + ops_name.push_back("SMOOTHING"); + } + + ops.emplace_back(rounding); + ops_name.emplace_back("rounding"); + } + + wildmeshing::utils::IntermediateWrite interm_writer( + mesh, + options.intermediate_output_path, + options.intermediate_output_name, + options.input_mesh_position, + options.intermediate_output); + interm_writer.disable_vertex_write(); + interm_writer.write(); + + ////////////////////////////////// + // Running all ops in order n times + Scheduler scheduler; + { + const size_t freq = options.scheduler_update_frequency; + scheduler.set_update_frequency(freq == 0 ? std::optional{} : freq); + } + + ////////////////////////////////// + // preprocessing + ////////////////////////////////// + + // debug code + test_flipped_triangles(mesh, pt_attribute); + + int64_t last_num_unrounded = std::numeric_limits::max(); + { + logger().info("----------------------- Preprocess Collapse -----------------------"); + + SchedulerStats pre_stats = + scheduler.run_operation_on_all(*collapse_then_round, visited_edge_flag.as()); + print_stats(pre_stats, "preprocessing collapse"); + + count_unrounded_vertices(pt_attribute, last_num_unrounded, false); + // test_flipped_triangles(mesh, pt_attribute); + compute_max_amips(amips_attribute, last_max_amips, false); + } + + interm_writer.write(); + + + for (int64_t i = 0; i < options.max_passes; ++i) { + logger().info("--------------------------- Pass {} ---------------------------", i); + + SchedulerStats pass_stats; + for (int64_t op_it = 0; op_it != ops.size(); ++op_it) { + auto& op = ops[op_it]; + auto& op_name = ops_name[op_it]; + + logger().debug("Executing {} ...", op_name); + + SchedulerStats stats; + if (op->primitive_type() == PrimitiveType::Edge) { + stats = scheduler.run_operation_on_all(*op, visited_edge_flag.as()); + } else { + stats = scheduler.run_operation_on_all(*op); + } + pass_stats += stats; + print_stats(stats, op_name); + + // count_unrounded_vertices(pt_attribute, last_num_unrounded, false); + // test_flipped_triangles(mesh, pt_attribute); + // compute_max_amips(amips_attribute, last_max_amips, false); + } + + print_stats(pass_stats); + + multimesh::consolidate(mesh); + + interm_writer.write(); + + assert(mesh.is_connectivity_valid()); + + const int64_t num_unrounded = count_unrounded_vertices(pt_attribute, last_num_unrounded); + const double max_amips = compute_max_amips(amips_attribute, last_max_amips); + + // stop at good quality + if (max_amips <= target_max_amips && num_unrounded == 0) { + break; + } + } + + logger().info("----------------------- Postprocess Collapse -----------------------"); + + // logger().info("Executing collapse ..."); + + auto post_stats = + scheduler.run_operation_on_all(*collapse_then_round, visited_edge_flag.as()); + print_stats(post_stats, "preprocessing collapse"); + + compute_max_amips(amips_attribute, last_max_amips); + + multimesh::consolidate(mesh); + + std::vector, std::string>> all_meshes; + all_meshes.push_back(std::make_pair(options.input_mesh, "main")); + + return all_meshes; +} + +} // namespace wmtk::components::internal diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.hpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.hpp new file mode 100644 index 0000000000..f9c91ab8e4 --- /dev/null +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_embedding_2d.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include +#include "WildmeshingOptions.hpp" + +namespace wmtk::components::internal { + +std::vector, std::string>> wildmeshing_embedding_2d( + const WildMeshingOptions& options); + +// void adjust_sizing_field( +// Mesh& m, +// const TypedAttributeHandle& coordinate_handle, +// const TypedAttributeHandle& edge_length_handle, +// const TypedAttributeHandle& sizing_field_scalar_handle, +// const TypedAttributeHandle& energy_handle, +// const TypedAttributeHandle& target_edge_length_handle, +// const TypedAttributeHandle& visited_handle, +// const double stop_energy, +// const double current_max_energy, +// const double initial_target_edge_length, +// const double min_target_edge_length); + +// void set_operation_energy_filter( +// Mesh& m, +// const TypedAttributeHandle& coordinate_handle, +// const TypedAttributeHandle& energy_handle, +// const TypedAttributeHandle& energy_filter_handle, +// const TypedAttributeHandle& visited_handle, +// const double stop_energy, +// const double current_max_energy, +// const double initial_target_edge_length); + +// void set_operation_energy_filter_after_sizing_field( +// Mesh& m, +// const TypedAttributeHandle& coordinate_handle, +// const TypedAttributeHandle& energy_handle, +// const TypedAttributeHandle& energy_filter_handle, +// const TypedAttributeHandle& visited_handle, +// const double stop_energy, +// const double current_max_energy, +// const double initial_target_edge_length); + +} // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.cpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.cpp index 65fc759895..b15b9217ee 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.cpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.cpp @@ -4,7 +4,7 @@ namespace wmtk::components::internal { void write( - const std::shared_ptr& mesh, + const Mesh& mesh, const std::string& out_dir, const std::string& name, const std::string& vname, @@ -12,43 +12,54 @@ void write( const bool intermediate_output) { if (intermediate_output) { - if (mesh->top_simplex_type() == PrimitiveType::Triangle) { + if (mesh.top_simplex_type() == PrimitiveType::Triangle) { // write trimesh const std::filesystem::path data_dir = ""; wmtk::io::ParaviewWriter writer( data_dir / (name + "_" + std::to_string(index)), vname, - *mesh, + mesh, true, true, true, false); - mesh->serialize(writer); - } else if (mesh->top_simplex_type() == PrimitiveType::Tetrahedron) { + mesh.serialize(writer); + } else if (mesh.top_simplex_type() == PrimitiveType::Tetrahedron) { // write tetmesh const std::filesystem::path data_dir = ""; wmtk::io::ParaviewWriter writer( data_dir / (name + "_" + std::to_string(index)), vname, - *mesh, + mesh, true, true, true, true); - mesh->serialize(writer); - } else if (mesh->top_simplex_type() == PrimitiveType::Edge) { + mesh.serialize(writer); + } else if (mesh.top_simplex_type() == PrimitiveType::Edge) { // write edgemesh const std::filesystem::path data_dir = ""; wmtk::io::ParaviewWriter writer( data_dir / (name + "_" + std::to_string(index)), vname, - *mesh, + mesh, true, true, false, false); - mesh->serialize(writer); + mesh.serialize(writer); } } } + +void write( + const std::shared_ptr& mesh, + const std::string& out_dir, + const std::string& name, + const std::string& vname, + const int64_t index, + const bool intermediate_output) +{ + write(*mesh, out_dir, name, vname, index, intermediate_output); +} } // namespace wmtk::components::internal \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.hpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.hpp index 4ad4961dfc..4ea06b4c69 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.hpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing_utils.hpp @@ -6,6 +6,14 @@ namespace wmtk::components::internal { +void write( + const Mesh& mesh, + const std::string& out_dir, + const std::string& name, + const std::string& vname, + const int64_t index, + const bool intermediate_output); + void write( const std::shared_ptr& mesh, const std::string& out_dir, diff --git a/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.cpp b/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.cpp new file mode 100644 index 0000000000..72b7730213 --- /dev/null +++ b/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.cpp @@ -0,0 +1,59 @@ +#include "IntermediateWrite.hpp" + +#include + +wmtk::components::wildmeshing::utils::IntermediateWrite::IntermediateWrite( + const Mesh& mesh, + const std::filesystem::path& out_dir, + const std::filesystem::path& name, + const std::string& vname, + const bool intermediate_output) + : m_mesh(mesh) + , m_name(out_dir / name) + , m_pos_attribute_name(vname) + , m_write(intermediate_output) +{ + switch (m_mesh.top_simplex_type()) { + case PrimitiveType::Tetrahedron: m_t = true; [[fallthrough]]; + case PrimitiveType::Triangle: m_f = true; [[fallthrough]]; + case PrimitiveType::Edge: m_e = true; [[fallthrough]]; + case PrimitiveType::Vertex: m_v = true; break; + default: log_and_throw_error("Unkown primitive type in IntermediateWrite"); + } +} + +void wmtk::components::wildmeshing::utils::IntermediateWrite::write() +{ + if (!m_write) { + return; + } + wmtk::io::ParaviewWriter writer( + fmt::format("{}_{}", m_name.string(), m_iter++), + m_pos_attribute_name, + m_mesh, + m_v, + m_e, + m_f, + m_t); + m_mesh.serialize(writer); +} + +void wmtk::components::wildmeshing::utils::IntermediateWrite::disable_vertex_write() +{ + m_v = false; +} + +void wmtk::components::wildmeshing::utils::IntermediateWrite::disable_edge_write() +{ + m_e = false; +} + +void wmtk::components::wildmeshing::utils::IntermediateWrite::disable_face_write() +{ + m_f = false; +} + +void wmtk::components::wildmeshing::utils::IntermediateWrite::disable_tetrahedron_write() +{ + m_t = false; +} diff --git a/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.hpp b/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.hpp new file mode 100644 index 0000000000..ad2a33c1ad --- /dev/null +++ b/components/wildmeshing/wmtk/components/wildmeshing/utils/IntermediateWrite.hpp @@ -0,0 +1,40 @@ +#pragma once + +#include +#include + +namespace wmtk::components::wildmeshing::utils { + +class IntermediateWrite +{ +public: + IntermediateWrite( + const Mesh& mesh, + const std::filesystem::path& out_dir, + const std::filesystem::path& name, + const std::string& vname, + const bool intermediate_output); + + IntermediateWrite(IntermediateWrite&) = delete; + + void write(); + + void disable_vertex_write(); + void disable_edge_write(); + void disable_face_write(); + void disable_tetrahedron_write(); + +private: + const Mesh& m_mesh; + const std::filesystem::path m_name; + const std::string m_pos_attribute_name; + const bool m_write; + int64_t m_iter = 0; + + bool m_v = false; + bool m_e = false; + bool m_f = false; + bool m_t = false; +}; + +} // namespace wmtk::components::wildmeshing::utils \ No newline at end of file diff --git a/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.cpp b/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.cpp index 274db376fe..47edf2f9eb 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.cpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.cpp @@ -2,18 +2,30 @@ #include "internal/wildmeshing2d.hpp" #include "internal/wildmeshing3d.hpp" +#include "internal/wildmeshing_embedding_2d.hpp" + +#include namespace wmtk::components { using namespace internal; std::vector, std::string>> wildmeshing( - const WildMeshingOptions& option) + const WildMeshingOptions& options) { - if (option.input_mesh->top_simplex_type() == PrimitiveType::Triangle) { - return wildmeshing2d(option); + if (options.use_embedding) { + if (options.input_mesh->top_simplex_type() == PrimitiveType::Triangle) { + return wildmeshing_embedding_2d(options); + } else { + log_and_throw_error("3D embedding not implemented yet."); + // return wildmeshing3d(option); + } + } + + if (options.input_mesh->top_simplex_type() == PrimitiveType::Triangle) { + return wildmeshing2d(options); } else { - return wildmeshing3d(option); + return wildmeshing3d(options); } assert(false); diff --git a/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.hpp b/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.hpp index df1ac057e8..2b5ea50dbb 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.hpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/wildmeshing.hpp @@ -5,10 +5,9 @@ #include "internal/WildmeshingOptions.hpp" - namespace wmtk::components { std::vector, std::string>> wildmeshing( - const WildMeshingOptions& option); + const WildMeshingOptions& options); } // namespace wmtk::components diff --git a/src/wmtk/CMakeLists.txt b/src/wmtk/CMakeLists.txt index a419a30f4e..6ecac967b9 100644 --- a/src/wmtk/CMakeLists.txt +++ b/src/wmtk/CMakeLists.txt @@ -2,6 +2,7 @@ set(SRC_FILES Cell.cpp Cell.hpp + MeshBase.hpp Mesh.cpp Mesh_attributes.cpp #the attribute parts of Mesh class Mesh_construction.cpp #the construction/destruction parts of Mesh class @@ -47,3 +48,4 @@ add_subdirectory(multimesh) # add_subdirectory(function) +add_subdirectory(submesh) \ No newline at end of file diff --git a/src/wmtk/EdgeMesh.hpp b/src/wmtk/EdgeMesh.hpp index d459a166ca..fae1b77af7 100644 --- a/src/wmtk/EdgeMesh.hpp +++ b/src/wmtk/EdgeMesh.hpp @@ -56,7 +56,7 @@ class EdgeMesh : public MeshCRTP Tuple switch_edge(const Tuple& tuple) const; std::vector orient_vertices(const Tuple& tuple) const override; - int64_t id(const Tuple& tuple, PrimitiveType type) const; + int64_t id(const Tuple& tuple, PrimitiveType type) const final override; using MeshCRTP::id; // getting the (simplex) prototype int64_t id_vertex(const Tuple& tuple) const { return id(tuple, PrimitiveType::Vertex); } diff --git a/src/wmtk/Mesh.hpp b/src/wmtk/Mesh.hpp index 1b556122e9..9444785fb2 100644 --- a/src/wmtk/Mesh.hpp +++ b/src/wmtk/Mesh.hpp @@ -18,6 +18,7 @@ // basic data for the class #include +#include "MeshBase.hpp" #include "Tuple.hpp" #include "Types.hpp" #include "attribute/AttributeManager.hpp" @@ -85,13 +86,19 @@ class TupleTag; } // namespace utils } // namespace multimesh +namespace submesh { +class Embedding; +} + // NOTE: the implementation of this class is split into several files to improve clang-format // performance // * Mesh.cpp // * Mesh_attributes.cpp // * Mesh_construction.cpp -class Mesh : public std::enable_shared_from_this, public wmtk::utils::MerkleTreeInteriorNode +class Mesh : public std::enable_shared_from_this, + public MeshBase, + public wmtk::utils::MerkleTreeInteriorNode { public: friend class tests::DEBUG_Mesh; @@ -113,6 +120,7 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk friend class multimesh::MultiMeshVisitorExecutor; friend class multimesh::attribute::AttributeScopeHandle; friend class multimesh::utils::internal::TupleTag; + friend class submesh::Embedding; friend class operations::utils::UpdateEdgeOperationMultiMeshMapFunctor; friend class operations::Operation; friend class operations::EdgeCollapse; @@ -121,10 +129,12 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk friend class operations::internal::CollapseAlternateFacetData; - int64_t top_cell_dimension() const; + int64_t top_cell_dimension() const override; PrimitiveType top_simplex_type() const; bool is_free() const; + MeshType mesh_type() const override; + // attribute directly hashes its "children" components so it overrides "child_hashes" std::map child_hashables() const override; std::map child_hashes() const override; @@ -147,9 +157,9 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk * @param type the type of tuple, can be vertex/edge/triangle/tetrahedron * @return vector of Tuples referring to each type */ - std::vector get_all(PrimitiveType type) const; + std::vector get_all(PrimitiveType type) const override; - std::vector get_all_id_simplex(PrimitiveType type) const; + std::vector get_all_id_simplex(PrimitiveType type) const override; /** * @brief Retrieve the IdSimplex that is represented by the tuple and primitive type. */ @@ -340,7 +350,7 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk public: /** * @brief switch the orientation of the Tuple of the given dimension - * @note this is not doen in place. Return a new Tuple of the switched state + * @note this is not done in place. Return a new Tuple of the switched state * * @param m * @param type d-0 -> switch vertex @@ -348,7 +358,7 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk d-2 -> switch face d-3 -> switch tetrahedron */ - virtual Tuple switch_tuple(const Tuple& tuple, PrimitiveType type) const = 0; + virtual Tuple switch_tuple(const Tuple& tuple, PrimitiveType type) const override = 0; // NOTE: adding per-simplex utility functions here is _wrong_ and will be removed @@ -411,7 +421,7 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk * @return true if this simplex lies on the boundary of the mesh * @return false otherwise */ - virtual bool is_boundary(PrimitiveType, const Tuple& tuple) const = 0; + virtual bool is_boundary(PrimitiveType, const Tuple& tuple) const override = 0; bool is_hash_valid(const Tuple& tuple, const attribute::Accessor& hash_accessor) const; @@ -775,6 +785,10 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk **/ bool is_from_same_multi_mesh_structure(const Mesh& other) const; + bool has_embedding() const; + + submesh::Embedding& get_embedding() const; + protected: // creates a scope as int64_t as the AttributeScopeHandle exists [[nodiscard]] attribute::AttributeScopeHandle create_single_mesh_scope(); @@ -792,13 +806,13 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk d-3 -> tetrahedron * @return int64_t id of the entity */ - int64_t id(const Tuple& tuple, PrimitiveType type) const; + int64_t id(const Tuple& tuple, PrimitiveType type) const override; int64_t id(const simplex::NavigatableSimplex& s) const { return s.index(); } int64_t id(const simplex::IdSimplex& s) const { return s.index(); } /// Forwarding version of id on simplices that does id caching - virtual int64_t id(const simplex::Simplex& s) const = 0; + virtual int64_t id(const simplex::Simplex& s) const override = 0; protected: /// Internal utility to allow id to be virtual with a non-virtual overload in derived -Mesh classes. @@ -828,6 +842,8 @@ class Mesh : public std::enable_shared_from_this, public wmtk::utils::Merk multimesh::MultiMeshManager m_multi_mesh_manager; + submesh::Embedding* m_embedding = nullptr; // a pointer to the embedding (if there is any) + int64_t m_top_cell_dimension = -1; // assumes no adjacency data exists @@ -973,6 +989,11 @@ inline bool Mesh::is_free() const return m_is_free; } +inline MeshType Mesh::mesh_type() const +{ + return MeshType::Mesh; +} + inline int64_t Mesh::top_cell_dimension() const { return m_top_cell_dimension; diff --git a/src/wmtk/MeshBase.hpp b/src/wmtk/MeshBase.hpp new file mode 100644 index 0000000000..0692291cb8 --- /dev/null +++ b/src/wmtk/MeshBase.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include +#include +#include "PrimitiveType.hpp" +#include "Tuple.hpp" + +namespace wmtk { + + +namespace simplex { +class Simplex; +class IdSimplex; +} // namespace simplex + +/** + * This enum is used to perform casts at runtime from the MeshBase type to the specific mesh type. + * + * Something like: + * if(emb_ptr->mesh_type() == MeshType::Embedding) { + * Embedding& emb = static_cast(*emb_ptr); + * } + */ +enum class MeshType { Mesh, Embedding, SubMesh }; + +class MeshBase +{ +public: + virtual ~MeshBase() = default; + + virtual std::vector get_all(PrimitiveType type) const = 0; + virtual std::vector get_all_id_simplex(PrimitiveType type) const = 0; + + virtual int64_t top_cell_dimension() const = 0; + + virtual MeshType mesh_type() const = 0; + + /** + * @brief switch the orientation of the Tuple of the given dimension + * @note this is not done in place. Return a new Tuple of the switched state + * + * @param m + * @param type d-0 -> switch vertex + d-1 -> switch edge + d-2 -> switch face + d-3 -> switch tetrahedron + */ + virtual Tuple switch_tuple(const Tuple& tuple, PrimitiveType type) const = 0; + + /** + * @brief check if a simplex (encoded as a tuple/primitive pair) lies on a boundary or not + * + * @param simplex + * @return true if this simplex lies on the boundary of the mesh + * @return false otherwise + */ + virtual bool is_boundary(PrimitiveType, const Tuple& tuple) const = 0; + + virtual int64_t id(const simplex::Simplex& s) const = 0; + virtual int64_t id(const Tuple& tuple, PrimitiveType pt) const = 0; +}; + +} // namespace wmtk \ No newline at end of file diff --git a/src/wmtk/MeshCRTP.hpp b/src/wmtk/MeshCRTP.hpp index ddfae5d357..2f022b9d40 100644 --- a/src/wmtk/MeshCRTP.hpp +++ b/src/wmtk/MeshCRTP.hpp @@ -100,7 +100,10 @@ class MeshCRTP : public Mesh public: /// Returns the id of a simplex encoded in a tuple - int64_t id(const Tuple& tuple, PrimitiveType type) const { return derived().id(tuple, type); } + int64_t id(const Tuple& tuple, PrimitiveType type) const override + { + return derived().id(tuple, type); + } /// variant of id that can cache internally held values int64_t id(const simplex::Simplex& s) const final override diff --git a/src/wmtk/Mesh_multimesh.cpp b/src/wmtk/Mesh_multimesh.cpp index 166f2ec389..932d96716e 100644 --- a/src/wmtk/Mesh_multimesh.cpp +++ b/src/wmtk/Mesh_multimesh.cpp @@ -1,6 +1,6 @@ #include "Mesh.hpp" -#include #include +#include #include #include @@ -36,6 +36,17 @@ bool Mesh::is_from_same_multi_mesh_structure(const Mesh& other) const return &get_multi_mesh_root() == &other.get_multi_mesh_root(); } +bool Mesh::has_embedding() const +{ + return m_embedding; +} + +submesh::Embedding& Mesh::get_embedding() const +{ + assert(has_embedding()); + return *m_embedding; +} + bool Mesh::can_map(const Mesh& other_mesh, const simplex::Simplex& my_simplex) const { if (!is_from_same_multi_mesh_structure(other_mesh)) { @@ -244,15 +255,15 @@ std::vector> Mesh::get_all_meshes() const auto meshes2 = get_all_child_meshes(); std::vector> meshes; meshes.emplace_back(shared_from_this()); - for(const auto& m: meshes2) { + for (const auto& m : meshes2) { meshes.emplace_back(m); } return meshes; - //std::queue> queue; + // std::queue> queue; ////std::queue queue; - //meshes.emplace_back(this); - //while(!queue.empty()) { - // const auto& cur = queue.front(); + // meshes.emplace_back(this); + // while(!queue.empty()) { + // const auto& cur = queue.front(); // //Mesh const* cur = queue.front(); // queue.pop(); // meshes.emplace_back(cur->shared_from_this()); @@ -260,6 +271,6 @@ std::vector> Mesh::get_all_meshes() const // queue.emplace(m.get()); // } //} - //return meshes; + // return meshes; } } // namespace wmtk diff --git a/src/wmtk/PointMesh.hpp b/src/wmtk/PointMesh.hpp index 87e5bcedc7..35c4684800 100644 --- a/src/wmtk/PointMesh.hpp +++ b/src/wmtk/PointMesh.hpp @@ -44,7 +44,7 @@ class PointMesh : public MeshCRTP std::vector orient_vertices(const Tuple& tuple) const override; using MeshCRTP::id; // getting the (simplex) prototype - int64_t id(const Tuple& tuple, PrimitiveType type) const; + int64_t id(const Tuple& tuple, PrimitiveType type) const override; protected: /** diff --git a/src/wmtk/TetMesh.hpp b/src/wmtk/TetMesh.hpp index 4a367a9b57..94022f8ad2 100644 --- a/src/wmtk/TetMesh.hpp +++ b/src/wmtk/TetMesh.hpp @@ -61,7 +61,7 @@ class TetMesh : public MeshCRTP std::vector orient_vertices(const Tuple& t) const override; - int64_t id(const Tuple& tuple, PrimitiveType type) const; + int64_t id(const Tuple& tuple, PrimitiveType type) const override; using MeshCRTP::id; // getting the (simplex) prototype diff --git a/src/wmtk/TriMesh.hpp b/src/wmtk/TriMesh.hpp index a003312156..fbbc5717e3 100644 --- a/src/wmtk/TriMesh.hpp +++ b/src/wmtk/TriMesh.hpp @@ -67,7 +67,7 @@ class TriMesh : public MeshCRTP std::vector orient_vertices(const Tuple& t) const override; - int64_t id(const Tuple& tuple, PrimitiveType type) const; + int64_t id(const Tuple& tuple, PrimitiveType type) const override; using MeshCRTP::id; // getting the (simplex) prototype int64_t id_vertex(const Tuple& tuple) const { return id(tuple, PrimitiveType::Vertex); } diff --git a/src/wmtk/attribute/AttributePrototypeDZ.hpp b/src/wmtk/attribute/AttributePrototypeDZ.hpp new file mode 100644 index 0000000000..d666c66d51 --- /dev/null +++ b/src/wmtk/attribute/AttributePrototypeDZ.hpp @@ -0,0 +1,116 @@ +template +class Attribute +{ +public: + using MapResult = internal::MapResult; + using ConstMapResult = internal::ConstMapResult; + + template + friend class AccessorBase; + friend class internal::AttributeTransactionStack; + + void serialize(const std::string& name, const int dim, MeshWriter& writer) const; + + Attribute(const std::string& name, int64_t dimension, T default_value = T(0), int64_t size = 0); + + Attribute(Attribute&& o); + ~Attribute(); + Attribute& operator=(Attribute&& o); + + ConstMapResult vector_attribute(const int64_t index) const; + MapResult vector_attribute(const int64_t index); + + T scalar_attribute(const int64_t index) const; + T& scalar_attribute(const int64_t index); + + int64_t reserved_size() const; + + int64_t dimension() const; + void reserve(const int64_t size); + + const T& default_value() const; + + bool operator==(const Attribute& o) const; + + void push_scope(); + void pop_scope(bool apply_updates); + void rollback_current_scope(); + + const AttributeTransactionStack& get_local_scope_stack() const; + AttributeTransactionStack& get_local_scope_stack(); + + void consolidate(const std::vector& new2old); + + void index_remap(const std::vector& old2new); + void index_remap(const std::vector& old2new, const std::vector& cols); + + + /////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////// + /////////////////////////////////////////////////////////// + // new stuff + + std::string name() const; + + struct SplitStrategies + { + SplitFuncType split_spine_strategy; + SplitRibFuncType split_rib_strategy; + }; + + void set_split_strategies(const SplitStrategies& strategies); + const SplitStrategies& get_split_strategies() const; + + struct CollapseStrategies + { + CollapseFuncType collapse_strategy; + }; + + void set_collapse_strategies(const CollapseStrategies& strategies); + const CollapseStrategies& get_collapse_strategies() const; + + struct SwapStrategies + { + SplitFuncType split_spine_strategy; + SplitRibFuncType split_rib_strategy; + CollapseFuncType collapse_strategy; + }; + + void set_swap_strategies(const SwapStrategies& strategies); + const SwapStrategies& get_swap_strategies() const; + + struct FaceSplitStrategies + { + SplitFuncType split_1_spine_strategy; + SplitRibFuncType split_1_rib_strategy; + SplitFuncType split_2_spine_strategy; + SplitRibFuncType split_2_rib_strategy; + CollapseFuncType collapse_strategy; + }; + + void set_face_split_strategies(const FaceSplitStrategies& strategies); + const FaceSplitStrategies& get_face_split_strategies() const; + + // set strategies to make it behave like a position + void set_default_position_strategies(); + // set strategies to make it behave like a tag + void set_default_tag_strategies(); + // set none strategies + void set_none_strategies(); + + // check if all strategies were set and print info + bool validate_strategies(const spdlog::log_level& l = spdlog::log_level::NONE); + + +private: + std::vector m_data; + PerThreadAttributeScopeStacks m_scope_stacks; + int64_t m_dimension = -1; + T m_default_value = T(0); + std::string m_name; + + SplitStrategies m_split_strategies; + CollapseStrategies m_collapse_strategies; + SwapStrategies m_swap_strategies; + FaceSplitStrategies m_face_split_strategies; +}; \ No newline at end of file diff --git a/src/wmtk/invariants/EnvelopeInvariant.cpp b/src/wmtk/invariants/EnvelopeInvariant.cpp index 98a0ea9ce9..95fb2c5240 100644 --- a/src/wmtk/invariants/EnvelopeInvariant.cpp +++ b/src/wmtk/invariants/EnvelopeInvariant.cpp @@ -9,6 +9,7 @@ #include "EnvelopeInvariant.hpp" #include +#include #include #include #include @@ -20,413 +21,613 @@ namespace wmtk { - -constexpr PrimitiveType PV = PrimitiveType::Vertex; +constexpr PrimitiveType PF = PrimitiveType::Triangle; constexpr PrimitiveType PE = PrimitiveType::Edge; -} // namespace wmtk -namespace wmtk::invariants { +constexpr PrimitiveType PV = PrimitiveType::Vertex; -EnvelopeInvariant::EnvelopeInvariant( +template +std::shared_ptr init_fast_envelope( const attribute::MeshAttributeHandle& envelope_mesh_coordinate, - double envelope_size, - const attribute::MeshAttributeHandle& coordinate) - : Invariant(coordinate.mesh()) - , m_coordinate_handle(coordinate) - , m_envelope_size(envelope_size) + const double envelope_size) { const auto& envelope_mesh = envelope_mesh_coordinate.mesh(); - assert(envelope_mesh_coordinate.holds() || envelope_mesh_coordinate.holds()); + const attribute::Accessor accessor = + envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as()); - if (envelope_mesh_coordinate.holds()) { - // for rational - const attribute::Accessor accessor = - envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as()); - - - if (envelope_mesh.top_simplex_type() == PrimitiveType::Triangle) { - std::vector vertices; - std::vector faces; - - int count = 0; - assert(accessor.dimension() == 3); - - const std::vector& facest = envelope_mesh.get_all(wmtk::PrimitiveType::Triangle); - for (const auto& f : facest) { - Eigen::Vector3d p0 = accessor.const_vector_attribute(f).cast(); - Eigen::Vector3d p1 = - accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV)) - .cast(); - Eigen::Vector3d p2 = - accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV})) - .cast(); - - faces.emplace_back(count, count + 1, count + 2); - vertices.push_back(p0); - vertices.push_back(p1); - vertices.push_back(p2); - - count += 3; - } - - m_envelope = - std::make_shared(vertices, faces, envelope_size); - - } else if (envelope_mesh.top_simplex_type() == PrimitiveType::Edge) { - logger().warn("Envelope for edge mesh is using sampling"); + std::vector vertices; + std::vector faces; - int64_t count = 0; - int64_t index = 0; + int count = 0; + assert(accessor.dimension() == 3); - const std::vector& edgest = envelope_mesh.get_all(wmtk::PrimitiveType::Edge); + const std::vector& facest = envelope_mesh.get_all(PF); + for (const auto& f : facest) { + if constexpr (std::is_same()) { + Eigen::Vector3d p0 = accessor.const_vector_attribute(f).template cast(); + Eigen::Vector3d p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV)) + .template cast(); + Eigen::Vector3d p2 = + accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV})) + .template cast(); - Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension()); - Eigen::MatrixXi edges(edgest.size(), 2); + vertices.push_back(p0); + vertices.push_back(p1); + vertices.push_back(p2); + } else { + Eigen::Vector3d p0 = accessor.const_vector_attribute(f); + Eigen::Vector3d p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV)); + Eigen::Vector3d p2 = + accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV})); + + vertices.push_back(p0); + vertices.push_back(p1); + vertices.push_back(p2); + } + faces.emplace_back(count, count + 1, count + 2); - for (const auto& e : edgest) { - auto p0 = accessor.const_vector_attribute(e).cast(); - auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV)) - .cast(); + count += 3; + } - edges.row(index) << count, count + 1; - vertices.row(2 * index) = p0; - vertices.row(2 * index + 1) = p1; + return std::make_shared(vertices, faces, envelope_size); +} - count += 2; - ++index; - } +template +std::shared_ptr init_bvh( + const attribute::MeshAttributeHandle& envelope_mesh_coordinate) +{ + const auto& envelope_mesh = envelope_mesh_coordinate.mesh(); - m_bvh = std::make_shared(); - m_bvh->init(vertices, edges, 1e-10); - } else { - throw std::runtime_error("Envelope works only for tri/edges meshes"); - } - } else if (envelope_mesh_coordinate.holds()) { - // for double - const attribute::Accessor accessor = - envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as()); + const attribute::Accessor accessor = + envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as()); + logger().warn("Envelope for edge mesh is using sampling"); - if (envelope_mesh.top_simplex_type() == PrimitiveType::Triangle) { - std::vector vertices; - std::vector faces; + int64_t count = 0; + int64_t index = 0; - int count = 0; - assert(accessor.dimension() == 3); + const std::vector& edgest = envelope_mesh.get_all(PE); - const std::vector& facest = envelope_mesh.get_all(wmtk::PrimitiveType::Triangle); - for (const auto& f : facest) { - Eigen::Vector3d p0 = accessor.const_vector_attribute(f); - Eigen::Vector3d p1 = - accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV)); - Eigen::Vector3d p2 = - accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV})); + Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension()); + Eigen::MatrixXi edges(edgest.size(), 2); - faces.emplace_back(count, count + 1, count + 2); - vertices.push_back(p0); - vertices.push_back(p1); - vertices.push_back(p2); + for (const Tuple& e : edgest) { + if constexpr (std::is_same()) { + const auto p0 = accessor.const_vector_attribute(e).template cast(); + const auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV)) + .template cast(); - count += 3; - } + vertices.row(2 * index) = p0; + vertices.row(2 * index + 1) = p1; + } else { + const auto p0 = accessor.const_vector_attribute(e); + const auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV)); - m_envelope = - std::make_shared(vertices, faces, envelope_size); + vertices.row(2 * index) = p0; + vertices.row(2 * index + 1) = p1; + } + edges.row(index) << count, count + 1; - } else if (envelope_mesh.top_simplex_type() == PrimitiveType::Edge) { - logger().warn("Envelope for edge mesh is using sampling"); + count += 2; + ++index; + } - int64_t count = 0; - int64_t index = 0; + std::shared_ptr bvh = std::make_shared(); + bvh->init(vertices, edges, 1e-10); + return bvh; +} - const std::vector& edgest = envelope_mesh.get_all(wmtk::PrimitiveType::Edge); +template +std::shared_ptr init_fast_envelope( + const attribute::MeshAttributeHandle& envelope_mesh_coordinate, + const double envelope_size, + const submesh::SubMesh& sub) +{ + const Mesh& m = envelope_mesh_coordinate.mesh(); - Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension()); - Eigen::MatrixXi edges(edgest.size(), 2); + const attribute::Accessor accessor = + m.create_const_accessor(envelope_mesh_coordinate.as()); - for (const auto& e : edgest) { - auto p0 = accessor.const_vector_attribute(e); - auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV)); + std::vector vertices; + std::vector faces; - edges.row(index) << count, count + 1; - vertices.row(2 * index) = p0; - vertices.row(2 * index + 1) = p1; + int count = 0; + assert(accessor.dimension() == 3); - count += 2; - ++index; - } + const std::vector& facest = sub.get_all(PF); + for (const auto& f : facest) { + if constexpr (std::is_same()) { + Eigen::Vector3d p0 = accessor.const_vector_attribute(f).template cast(); + Eigen::Vector3d p1 = + accessor.const_vector_attribute(sub.switch_tuple(f, PV)).template cast(); + Eigen::Vector3d p2 = accessor.const_vector_attribute(sub.switch_tuples(f, {PE, PV})) + .template cast(); - m_bvh = std::make_shared(); - m_bvh->init(vertices, edges, 1e-10); + vertices.push_back(p0); + vertices.push_back(p1); + vertices.push_back(p2); } else { - throw std::runtime_error("Envelope works only for tri/edges meshes"); + Eigen::Vector3d p0 = accessor.const_vector_attribute(f); + Eigen::Vector3d p1 = accessor.const_vector_attribute(sub.switch_tuple(f, PV)); + Eigen::Vector3d p2 = accessor.const_vector_attribute(sub.switch_tuples(f, {PE, PV})); + + vertices.push_back(p0); + vertices.push_back(p1); + vertices.push_back(p2); } - } else { - throw std::runtime_error("Envelope mesh handle type invlid"); + faces.emplace_back(count, count + 1, count + 2); + + count += 3; } + + return std::make_shared(vertices, faces, envelope_size); } -bool EnvelopeInvariant::after( - const std::vector& top_dimension_tuples_before, - const std::vector& top_dimension_tuples_after) const +template +std::shared_ptr init_bvh( + const attribute::MeshAttributeHandle& envelope_mesh_coordinate, + const submesh::SubMesh& sub) { - if (top_dimension_tuples_after.empty()) return true; + const Mesh& m = envelope_mesh_coordinate.mesh(); - assert(m_coordinate_handle.holds() || m_coordinate_handle.holds()); + const attribute::Accessor accessor = + m.create_const_accessor(envelope_mesh_coordinate.as()); - if (m_coordinate_handle.holds()) { - const attribute::Accessor accessor = - mesh().create_const_accessor(m_coordinate_handle.as()); - const PrimitiveType type = mesh().top_simplex_type(); + logger().warn("Envelope for edge mesh is using sampling"); - if (m_envelope) { - assert(accessor.dimension() == 3); + int64_t count = 0; + int64_t index = 0; - std::vector faces; + const std::vector& edgest = sub.get_all(PE); - if (type == PrimitiveType::Triangle) { - std::array triangle; + Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension()); + Eigen::MatrixXi edges(edgest.size(), 2); - for (const Tuple& tuple : top_dimension_tuples_after) { - faces = faces_single_dimension_tuples( - mesh(), - simplex::Simplex(mesh(), type, tuple), - PrimitiveType::Vertex); + for (const Tuple& e : edgest) { + if constexpr (std::is_same()) { + const auto p0 = accessor.const_vector_attribute(e).template cast(); + const auto p1 = + accessor.const_vector_attribute(sub.switch_tuple(e, PV)).template cast(); - triangle[0] = accessor.const_vector_attribute(faces[0]).cast(); - triangle[1] = accessor.const_vector_attribute(faces[1]).cast(); - triangle[2] = accessor.const_vector_attribute(faces[2]).cast(); + vertices.row(2 * index) = p0; + vertices.row(2 * index + 1) = p1; + } else { + const auto p0 = accessor.const_vector_attribute(e); + const auto p1 = accessor.const_vector_attribute(sub.switch_tuple(e, PV)); - if (m_envelope->is_outside(triangle)) { - wmtk::logger().debug("fail envelope check 1"); - return false; - } - } + vertices.row(2 * index) = p0; + vertices.row(2 * index + 1) = p1; + } + edges.row(index) << count, count + 1; - return true; - } else if (type == PrimitiveType::Edge) { - for (const Tuple& tuple : top_dimension_tuples_after) { - faces = faces_single_dimension_tuples( - mesh(), - simplex::Simplex(mesh(), type, tuple), - PrimitiveType::Vertex); - - Eigen::Vector3d p0 = accessor.const_vector_attribute(faces[0]).cast(); - Eigen::Vector3d p1 = accessor.const_vector_attribute(faces[1]).cast(); - - if (m_envelope->is_outside(p0, p1)) { - wmtk::logger().debug("fail envelope check 2"); - return false; - } - } + count += 2; + ++index; + } - return true; - } else if (type == PrimitiveType::Vertex) { - for (const Tuple& tuple : top_dimension_tuples_after) { - Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast(); + std::shared_ptr bvh = std::make_shared(); + bvh->init(vertices, edges, 1e-10); + return bvh; +} - if (m_envelope->is_outside(p)) { - wmtk::logger().debug("fail envelope check 3"); - return false; - } - } +} // namespace wmtk +namespace wmtk::invariants { - return true; - } else { - throw std::runtime_error("Invalid mesh type"); - } - return true; - } else { - assert(m_bvh); +EnvelopeInvariant::EnvelopeInvariant( + const attribute::MeshAttributeHandle& envelope_mesh_coordinate, + double envelope_size, + const attribute::MeshAttributeHandle& coordinate) + : Invariant(coordinate.mesh(), false, false, true) + , m_coordinate_handle(coordinate) + , m_envelope_size(envelope_size) +{ + const auto& envelope_mesh = envelope_mesh_coordinate.mesh(); - SimpleBVH::VectorMax3d nearest_point; - double sq_dist; + // log_assert( + // envelope_mesh_coordinate.holds() || envelope_mesh_coordinate.holds(), + // "Envelope mesh handle type invalid"); - const double d = m_envelope_size; - const double real_envelope = m_envelope_size - d / sqrt(accessor.dimension()); - const double real_envelope_2 = real_envelope * real_envelope; + if (!(envelope_mesh_coordinate.holds() || envelope_mesh_coordinate.holds())) { + log_and_throw_error("Envelope mesh handle type invalid"); + } - if (type == PrimitiveType::Edge) { - std::vector pts; + // log_assert( + // envelope_mesh.top_simplex_type() == PF || envelope_mesh.top_simplex_type() == PE, + // "Envelope works only for tri/edges meshes"); - for (const Tuple& tuple : top_dimension_tuples_after) { - const auto p0 = accessor.const_vector_attribute(tuple).cast().eval(); - const auto p1 = accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV)) - .cast() - .eval(); + if (!(envelope_mesh.top_simplex_type() == PF || envelope_mesh.top_simplex_type() == PE)) { + log_and_throw_error("Envelope works only for tri/edges meshes"); + } - const int64_t N = (p0 - p1).norm() / d + 1; - pts.reserve(pts.size() + N); + if (envelope_mesh_coordinate.holds()) { + if (envelope_mesh.top_simplex_type() == PF) { + m_envelope = init_fast_envelope(envelope_mesh_coordinate, envelope_size); + } else if (envelope_mesh.top_simplex_type() == PE) { + m_bvh = init_bvh(envelope_mesh_coordinate); + } + } else if (envelope_mesh_coordinate.holds()) { + if (envelope_mesh.top_simplex_type() == PF) { + m_envelope = init_fast_envelope(envelope_mesh_coordinate, envelope_size); + } else if (envelope_mesh.top_simplex_type() == PE) { + m_bvh = init_bvh(envelope_mesh_coordinate); + } + } +} - for (int64_t n = 0; n <= N; n++) { - auto tmp = p0 * (double(n) / N) + p1 * (N - double(n)) / N; - pts.push_back(tmp); - } - } +EnvelopeInvariant::EnvelopeInvariant( + const attribute::MeshAttributeHandle& pt_attribute, + double envelope_size, + const submesh::SubMesh& sub) + : Invariant(pt_attribute.mesh(), false, false, true) + , m_coordinate_handle(pt_attribute) + , m_submesh(&sub) + , m_envelope_size(envelope_size) +{ + // log_assert( + // envelope_mesh_coordinate.holds() || envelope_mesh_coordinate.holds(), + // "Envelope mesh handle type invalid"); - auto current_point = pts[0]; + if (!(pt_attribute.holds() || pt_attribute.holds())) { + log_and_throw_error("Envelope mesh handle type invalid"); + } - int prev_facet = m_bvh->nearest_facet(current_point, nearest_point, sq_dist); - if (sq_dist > real_envelope_2) { - wmtk::logger().debug("fail envelope check 4"); - return false; - } + // log_assert( + // envelope_mesh.top_simplex_type() == PF || envelope_mesh.top_simplex_type() == PE, + // "Envelope works only for tri/edges meshes"); - for (const auto& v : pts) { - sq_dist = (v - nearest_point).squaredNorm(); - m_bvh->nearest_facet_with_hint(v, prev_facet, nearest_point, sq_dist); - if (sq_dist > real_envelope_2) { - wmtk::logger().debug("fail envelope check 5"); - return false; - } - } + if (!(sub.top_simplex_type() == PF || sub.top_simplex_type() == PE)) { + log_and_throw_error("Envelope works only for tri/edges meshes"); + } - return true; - } else if (type == PrimitiveType::Vertex) { - for (const Tuple& tuple : top_dimension_tuples_after) { - Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast(); - m_bvh->nearest_facet(p, nearest_point, sq_dist); - if (sq_dist > m_envelope_size * m_envelope_size) { - wmtk::logger().debug("fail envelope check 6"); - return false; - } + if (pt_attribute.holds()) { + if (sub.top_simplex_type() == PF) { + m_envelope = init_fast_envelope(pt_attribute, envelope_size, sub); + } else if (sub.top_simplex_type() == PE) { + m_bvh = init_bvh(pt_attribute, sub); + } + } else if (pt_attribute.holds()) { + if (sub.top_simplex_type() == PF) { + m_envelope = init_fast_envelope(pt_attribute, envelope_size, sub); + } else if (sub.top_simplex_type() == PE) { + m_bvh = init_bvh(pt_attribute, sub); + } + } +} + +bool EnvelopeInvariant::after( + const std::vector& top_dimension_tuples_before, + const std::vector& top_dimension_tuples_after) const +{ + if (top_dimension_tuples_after.empty()) { + return true; + } + + assert(m_coordinate_handle.holds() || m_coordinate_handle.holds()); + assert(m_envelope || m_bvh); + + /* + Modification for submesh: + - Get faces of top_dimension_tuples_after that are of type sub.top_simplex_type() + - Filter by sub.contains(tuple, pt_top) + */ + + if (m_submesh) { + const submesh::SubMesh& sub = *m_submesh; + + const PrimitiveType pt_top = sub.top_simplex_type(); + + simplex::IdSimplexCollection tops(mesh()); + for (const Tuple& t : top_dimension_tuples_after) { + if (mesh().top_simplex_type() == pt_top) { + if (sub.contains(t, pt_top)) { + tops.add(pt_top, t); } + continue; + } + + const auto faces = simplex::faces_single_dimension_tuples( + mesh(), + simplex::Simplex(mesh().top_simplex_type(), t), + pt_top); - return true; - } else { - throw std::runtime_error("Invalid mesh type"); + for (const Tuple& f : faces) { + if (sub.contains(f, pt_top)) { + tops.add(pt_top, f); + } } + } + tops.sort_and_clean(); + + if (tops.size() == 0) { return true; } - } else if (m_coordinate_handle.holds()) { - const attribute::Accessor accessor = - mesh().create_const_accessor(m_coordinate_handle.as()); - const auto type = mesh().top_simplex_type(); if (m_envelope) { - assert(accessor.dimension() == 3); + return after_with_envelope(tops.simplex_vector_tuples()); + } + if (m_bvh) { + return after_with_bvh(tops.simplex_vector_tuples()); + } + assert(false); // this code should be unreachable + return false; + } - std::vector faces; + if (m_envelope) { + return after_with_envelope(top_dimension_tuples_after); + } + if (m_bvh) { + return after_with_bvh(top_dimension_tuples_after); + } + assert(false); // this code should be unreachable + return false; +} + + +bool EnvelopeInvariant::after_with_envelope( + const std::vector& top_dimension_tuples_after) const +{ + assert(m_envelope); + assert(m_coordinate_handle.dimension() == 3); + + const PrimitiveType pt_top = + m_submesh ? m_submesh->top_simplex_type() : mesh().top_simplex_type(); + + if (pt_top == PrimitiveType::Triangle) { + return after_with_envelope_triangle(top_dimension_tuples_after); + } + if (pt_top == PrimitiveType::Edge) { + return after_with_envelope_edge(top_dimension_tuples_after); + } + if (pt_top == PrimitiveType::Vertex) { + return after_with_envelope_vertex(top_dimension_tuples_after); + } + log_and_throw_error("Invalid mesh type"); +} + +bool EnvelopeInvariant::after_with_envelope_triangle( + const std::vector& top_dimension_tuples_after) const +{ + assert(m_envelope); - if (type == PrimitiveType::Triangle) { - std::array triangle; + const bool res = std::visit( + [&](auto&& tah) noexcept { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; - for (const Tuple& tuple : top_dimension_tuples_after) { - faces = faces_single_dimension_tuples( - mesh(), - simplex::Simplex(mesh(), type, tuple), - PrimitiveType::Vertex); + const auto accessor = mesh().create_const_accessor(tah); + std::array triangle; + + for (const Tuple& tuple : top_dimension_tuples_after) { + const auto faces = faces_single_dimension_tuples( + mesh(), + simplex::Simplex(PrimitiveType::Triangle, tuple), + PrimitiveType::Vertex); + + if constexpr (std::is_same_v) { triangle[0] = accessor.const_vector_attribute(faces[0]); triangle[1] = accessor.const_vector_attribute(faces[1]); triangle[2] = accessor.const_vector_attribute(faces[2]); + } + if constexpr (std::is_same_v) { + triangle[0] = accessor.const_vector_attribute(faces[0]).template cast(); + triangle[1] = accessor.const_vector_attribute(faces[1]).template cast(); + triangle[2] = accessor.const_vector_attribute(faces[2]).template cast(); + } + + if (m_envelope->is_outside(triangle)) { + wmtk::logger().debug("fail envelope check 1"); + return false; + } + } + return true; + }, + m_coordinate_handle.handle()); + + return res; +} - if (m_envelope->is_outside(triangle)) { - wmtk::logger().debug("fail envelope check 7"); - return false; - } +bool EnvelopeInvariant::after_with_envelope_edge( + const std::vector& top_dimension_tuples_after) const +{ + assert(m_envelope); + assert(mesh().top_simplex_type() == PrimitiveType::Edge); + + const bool res = std::visit( + [&](auto&& tah) noexcept { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; + + const auto accessor = mesh().create_const_accessor(tah); + + std::array triangle; + const PrimitiveType pt_top = mesh().top_simplex_type(); + + for (const Tuple& tuple : top_dimension_tuples_after) { + const auto faces = faces_single_dimension_tuples( + mesh(), + simplex::Simplex(mesh(), pt_top, tuple), + PrimitiveType::Vertex); + + Eigen::Vector3d p0; + Eigen::Vector3d p1; + + if constexpr (std::is_same_v) { + p0 = accessor.const_vector_attribute(faces[0]); + p1 = accessor.const_vector_attribute(faces[1]); + } else if constexpr (std::is_same_v) { + p0 = accessor.const_vector_attribute(faces[0]).template cast(); + p1 = accessor.const_vector_attribute(faces[1]).template cast(); + } else { + log_and_throw_error("Unknown attribute type"); } - return true; - } else if (type == PrimitiveType::Edge) { - for (const Tuple& tuple : top_dimension_tuples_after) { - faces = faces_single_dimension_tuples( - mesh(), - simplex::Simplex(mesh(), type, tuple), - PrimitiveType::Vertex); - - Eigen::Vector3d p0 = accessor.const_vector_attribute(faces[0]); - Eigen::Vector3d p1 = accessor.const_vector_attribute(faces[1]); - - if (m_envelope->is_outside(p0, p1)) { - wmtk::logger().debug("fail envelope check 8"); - return false; - } + if (m_envelope->is_outside(p0, p1)) { + wmtk::logger().debug("fail envelope check 2"); + return false; } + } - return true; - } else if (type == PrimitiveType::Vertex) { - for (const Tuple& tuple : top_dimension_tuples_after) { - Eigen::Vector3d p = accessor.const_vector_attribute(tuple); + return true; + }, + m_coordinate_handle.handle()); - if (m_envelope->is_outside(p)) { - wmtk::logger().debug("fail envelope check 9"); - return false; - } + return res; +} + +bool EnvelopeInvariant::after_with_envelope_vertex( + const std::vector& top_dimension_tuples_after) const +{ + assert(m_envelope); + assert(mesh().top_simplex_type() == PrimitiveType::Vertex); + + const bool res = std::visit( + [&](auto&& tah) noexcept { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; + + const auto accessor = mesh().create_const_accessor(tah); + + for (const Tuple& tuple : top_dimension_tuples_after) { + Eigen::Vector3d p; + if constexpr (std::is_same_v) { + p = accessor.const_vector_attribute(tuple); + } else if constexpr (std::is_same_v) { + p = accessor.const_vector_attribute(tuple).template cast(); + } else { + log_and_throw_error("Unknown attribute type"); } - return true; - } else { - throw std::runtime_error("Invalid mesh type"); + if (m_envelope->is_outside(p)) { + wmtk::logger().debug("fail envelope check 3"); + return false; + } } + return true; - } else { - assert(m_bvh); + }, + m_coordinate_handle.handle()); - SimpleBVH::VectorMax3d nearest_point; - double sq_dist; + return res; +} - const double d = m_envelope_size; - const double real_envelope = m_envelope_size - d / sqrt(accessor.dimension()); - const double real_envelope_2 = real_envelope * real_envelope; +bool EnvelopeInvariant::after_with_bvh(const std::vector& top_dimension_tuples_after) const +{ + assert(m_bvh); - if (type == PrimitiveType::Edge) { - std::vector pts; + const PrimitiveType pt_top = + m_submesh ? m_submesh->top_simplex_type() : mesh().top_simplex_type(); - for (const Tuple& tuple : top_dimension_tuples_after) { - SimpleBVH::VectorMax3d p0 = - accessor.const_vector_attribute(tuple).cast(); - SimpleBVH::VectorMax3d p1 = - accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV)) - .cast(); + if (pt_top == PrimitiveType::Edge) { + return after_with_bvh_edge(top_dimension_tuples_after); + } + if (pt_top == PrimitiveType::Vertex) { + return after_with_bvh_vertex(top_dimension_tuples_after); + } - const int64_t N = (p0 - p1).norm() / d + 1; - pts.reserve(pts.size() + N); + log_and_throw_error("Invalid mesh type"); +} - for (int64_t n = 0; n <= N; n++) { - auto tmp = p0 * (double(n) / N) + p1 * (N - double(n)) / N; - pts.push_back(tmp); - } +bool EnvelopeInvariant::after_with_bvh_edge( + const std::vector& top_dimension_tuples_after) const +{ + SimpleBVH::VectorMax3d nearest_point; + double sq_dist; + + const double d = m_envelope_size; + const double real_envelope = m_envelope_size - d / sqrt(m_coordinate_handle.dimension()); + const double real_envelope_2 = real_envelope * real_envelope; + + std::vector pts; + + std::visit( + [&](auto&& tah) noexcept { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; + + const auto accessor = mesh().create_const_accessor(tah); + + for (const Tuple& tuple : top_dimension_tuples_after) { + SimpleBVH::VectorMax3d p0; + SimpleBVH::VectorMax3d p1; + if constexpr (std::is_same_v) { + p0 = accessor.const_vector_attribute(tuple); + p1 = accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV)); + } else if constexpr (std::is_same_v) { + p0 = accessor.const_vector_attribute(tuple).template cast(); + p1 = accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV)) + .template cast(); + } else { + log_and_throw_error("Unknown attribute type"); } - auto current_point = pts[0]; + const int64_t N = (p0 - p1).norm() / d + 1; + pts.reserve(pts.size() + N); - int prev_facet = m_bvh->nearest_facet(current_point, nearest_point, sq_dist); - if (sq_dist > real_envelope_2) { - wmtk::logger().debug("fail envelope check 10"); - return false; + for (int64_t n = 0; n <= N; n++) { + auto tmp = p0 * (double(n) / N) + p1 * (N - double(n)) / N; + pts.push_back(tmp); } + } + }, + m_coordinate_handle.handle()); - for (const auto& v : pts) { - sq_dist = (v - nearest_point).squaredNorm(); - m_bvh->nearest_facet_with_hint(v, prev_facet, nearest_point, sq_dist); - if (sq_dist > real_envelope_2) { - wmtk::logger().debug("fail envelope check 11"); - return false; - } - } + SimpleBVH::VectorMax3d current_point = pts[0]; - return true; - } else if (type == PrimitiveType::Vertex) { - for (const Tuple& tuple : top_dimension_tuples_after) { - Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast(); - m_bvh->nearest_facet(p, nearest_point, sq_dist); - if (sq_dist > m_envelope_size * m_envelope_size) { - wmtk::logger().debug("fail envelope check 12"); - return false; - } - } + int prev_facet = m_bvh->nearest_facet(current_point, nearest_point, sq_dist); + if (sq_dist > real_envelope_2) { + wmtk::logger().debug("fail envelope check 4"); + return false; + } - return true; - } else { - throw std::runtime_error("Invalid mesh type"); - } - return true; + for (const auto& v : pts) { + sq_dist = (v - nearest_point).squaredNorm(); + m_bvh->nearest_facet_with_hint(v, prev_facet, nearest_point, sq_dist); + if (sq_dist > real_envelope_2) { + wmtk::logger().debug("fail envelope check 5"); + return false; } - } else { - throw std::runtime_error("Envelope mesh handle type invlid"); } + + return true; } +bool EnvelopeInvariant::after_with_bvh_vertex( + const std::vector& top_dimension_tuples_after) const +{ + assert(m_bvh); + assert(mesh().top_simplex_type() == PrimitiveType::Vertex); + + SimpleBVH::VectorMax3d nearest_point; + double sq_dist; + + for (const Tuple& tuple : top_dimension_tuples_after) { + SimpleBVH::VectorMax3d p = std::visit( + [&](auto&& tah) noexcept { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; + + const auto accessor = mesh().create_const_accessor(tah); + + SimpleBVH::VectorMax3d _p; + if constexpr (std::is_same_v) { + _p = accessor.const_vector_attribute(tuple); + } else if constexpr (std::is_same_v) { + _p = accessor.const_vector_attribute(tuple).template cast(); + } else { + log_and_throw_error("Unknown attribute type"); + } + return _p; + }, + m_coordinate_handle.handle()); + + m_bvh->nearest_facet(p, nearest_point, sq_dist); + if (sq_dist > m_envelope_size * m_envelope_size) { + wmtk::logger().debug("fail envelope check 6"); + return false; + } + } + + return true; +} } // namespace wmtk::invariants diff --git a/src/wmtk/invariants/EnvelopeInvariant.hpp b/src/wmtk/invariants/EnvelopeInvariant.hpp index ea9f6b4fcf..570e1fb603 100644 --- a/src/wmtk/invariants/EnvelopeInvariant.hpp +++ b/src/wmtk/invariants/EnvelopeInvariant.hpp @@ -5,6 +5,7 @@ #include #include +#include #include namespace fastEnvelope { @@ -19,20 +20,47 @@ namespace wmtk::invariants { class EnvelopeInvariant : public Invariant { public: + /** + * @brief Creates an envelope for checking if vertices are inside or outside of it. + * + * The check is performed on the `coordinate`, the envelope is constructed from + * `envelope_mesh_coordinate`. + * + * @param envelope_mesh_coordinate Used for constructing the envelope. + * @param envelope_size + * @param coordinate This position handle represents the mesh the envelope is tested on. + */ EnvelopeInvariant( const attribute::MeshAttributeHandle& envelope_mesh_coordinate, double envelope_size, const attribute::MeshAttributeHandle& coordinate); + EnvelopeInvariant( + const attribute::MeshAttributeHandle& pt_attribute, + double envelope_size, + const submesh::SubMesh& sub); + bool after( const std::vector& top_dimension_tuples_before, const std::vector& top_dimension_tuples_after) const override; +private: + bool after_with_envelope(const std::vector& top_dimension_tuples_after) const; + bool after_with_envelope_triangle(const std::vector& top_dimension_tuples_after) const; + bool after_with_envelope_edge(const std::vector& top_dimension_tuples_after) const; + bool after_with_envelope_vertex(const std::vector& top_dimension_tuples_after) const; + + bool after_with_bvh(const std::vector& top_dimension_tuples_after) const; + bool after_with_bvh_edge(const std::vector& top_dimension_tuples_after) const; + bool after_with_bvh_vertex(const std::vector& top_dimension_tuples_after) const; + private: std::shared_ptr m_envelope = nullptr; std::shared_ptr m_bvh = nullptr; const attribute::MeshAttributeHandle m_coordinate_handle; + const submesh::SubMesh* m_submesh = nullptr; + const double m_envelope_size; }; } // namespace wmtk::invariants diff --git a/src/wmtk/io/ParaviewWriter.cpp b/src/wmtk/io/ParaviewWriter.cpp index 8f71e6c39a..c0f219173a 100644 --- a/src/wmtk/io/ParaviewWriter.cpp +++ b/src/wmtk/io/ParaviewWriter.cpp @@ -2,6 +2,7 @@ #include +#include #include #include @@ -96,7 +97,8 @@ ParaviewWriter::ParaviewWriter( bool write_points, bool write_edges, bool write_faces, - bool write_tetrahedra) + bool write_tetrahedra, + const std::function& filter) : m_vertices_name(vertices_name) { m_enabled[0] = write_points; @@ -106,41 +108,38 @@ ParaviewWriter::ParaviewWriter( std::array cells; - for (size_t i = 0; i < 4; ++i) { - const auto pt = PrimitiveType(i); - if (m_enabled[i]) { - // include deleted tuples so that attributes are aligned - const auto tuples = mesh.get_all(pt, true); - cells[i].resize(tuples.size(), i + 1); - - for (size_t j = 0; j < tuples.size(); ++j) { - const auto& t = tuples[j]; - if (t.is_null()) { - for (size_t d = 0; d < cells[i].cols(); ++d) { - cells[i](j, d) = 0; - } - } else { - int64_t vid = mesh.id(t, PrimitiveType::Vertex); - cells[i](j, 0) = vid; - if (i > 0) { - auto t1 = mesh.switch_tuple(t, PrimitiveType::Vertex); - - cells[i](j, 1) = mesh.id(t1, PrimitiveType::Vertex); - } - if (i > 1) { - auto t1 = mesh.switch_tuple(t, PrimitiveType::Edge); - auto t2 = mesh.switch_tuple(t1, PrimitiveType::Vertex); - - cells[i](j, 2) = mesh.id(t2, PrimitiveType::Vertex); - } - if (i > 2) { - auto t1 = mesh.switch_tuple(t, PrimitiveType::Triangle); - auto t2 = mesh.switch_tuple(t1, PrimitiveType::Edge); - auto t3 = mesh.switch_tuple(t2, PrimitiveType::Vertex); - - cells[i](j, 3) = mesh.id(t3, PrimitiveType::Vertex); - } - } + for (size_t dim = 0; dim < 4; ++dim) { + const PrimitiveType pt = PrimitiveType(dim); + if (!m_enabled[dim]) { + continue; + } + // include deleted simplices so that attributes are aligned + + const auto simplices = mesh.get_all_id_simplex(pt, true); + + cells[dim].resize(simplices.size(), dim + 1); + + for (int64_t s_id = 0; s_id < simplices.size(); ++s_id) { + const simplex::IdSimplex& s = simplices[s_id]; + + if (s.index() != s_id || (filter != nullptr && !filter(s))) { + // deleted simplex + cells[dim].row(s_id).setZero(); + continue; + } + + if (dim == 0) { + // vertex + cells[dim](s_id, 0) = s.index(); + continue; + } + + // edge, triangle, tetrahedron + const auto vertices = + simplex::faces_single_dimension(mesh, mesh.get_simplex(s), PrimitiveType::Vertex); + assert(vertices.size() == dim + 1); + for (int64_t i = 0; i < vertices.size(); ++i) { + cells[dim](s_id, i) = mesh.id(vertices.simplex_vector()[i]); } } } diff --git a/src/wmtk/io/ParaviewWriter.hpp b/src/wmtk/io/ParaviewWriter.hpp index 925bc2868d..9fb8c6cf41 100644 --- a/src/wmtk/io/ParaviewWriter.hpp +++ b/src/wmtk/io/ParaviewWriter.hpp @@ -4,6 +4,7 @@ #include #include +#include namespace paraviewo { class ParaviewWriter; @@ -11,6 +12,9 @@ class ParaviewWriter; namespace wmtk { class Mesh; +namespace simplex { +class IdSimplex; +} namespace io { class ParaviewWriter : public MeshWriter @@ -51,6 +55,20 @@ class ParaviewWriter : public MeshWriter }; public: + /** + * @brief Write in VTU format. + * + * The writer generates one file for each simplex dimension and attaches all the attributes for + * the corresponding simplex dimension. All higher dimensions also contain the simplex + * attributes. + * + * The writer stores ALL simplices, even those that are invalid. This helps with debugging, as + * the IDs in the VTU file correspond to those in the code. All invalid simplices will contain 0 + * vertex IDs, so in case one vertex looks strange in Paraview, that might be because of that. + * + * The filter function can be used to treat simplices as invalid ones. A simplex is treated as + * invalid if the filter function returns false. + */ ParaviewWriter( const std::filesystem::path& filename, const std::string& vertices_name, @@ -58,7 +76,8 @@ class ParaviewWriter : public MeshWriter bool write_points = true, bool write_edges = true, bool write_faces = true, - bool write_tetrahedra = true); + bool write_tetrahedra = true, + const std::function& filter = {}); bool write(const int dim) override { return dim == 0 || m_enabled[dim]; } diff --git a/src/wmtk/multimesh/utils/extract_child_mesh_from_tag.cpp b/src/wmtk/multimesh/utils/extract_child_mesh_from_tag.cpp index 3181630af1..f9637151c6 100644 --- a/src/wmtk/multimesh/utils/extract_child_mesh_from_tag.cpp +++ b/src/wmtk/multimesh/utils/extract_child_mesh_from_tag.cpp @@ -5,7 +5,6 @@ #include #include #include -#include #include #include #include diff --git a/src/wmtk/multimesh/utils/internal/TupleTag.cpp b/src/wmtk/multimesh/utils/internal/TupleTag.cpp index d07824280c..cf719e355b 100644 --- a/src/wmtk/multimesh/utils/internal/TupleTag.cpp +++ b/src/wmtk/multimesh/utils/internal/TupleTag.cpp @@ -1,6 +1,6 @@ #include "TupleTag.hpp" +#include #include -#include namespace wmtk::multimesh::utils::internal { TupleTag::TupleTag(Mesh& mesh, const std::set& critical_points) : m_mesh(mesh) diff --git a/src/wmtk/operations/AMIPSOptimizationSmoothing.cpp b/src/wmtk/operations/AMIPSOptimizationSmoothing.cpp index b12a76f377..847ade60dd 100644 --- a/src/wmtk/operations/AMIPSOptimizationSmoothing.cpp +++ b/src/wmtk/operations/AMIPSOptimizationSmoothing.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -259,7 +260,7 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: if (m_coordinate_handle.holds()) { auto accessor = mesh().create_accessor(m_coordinate_handle.as()); - const auto neighs = wmtk::simplex::cofaces_single_dimension_simplices( + auto neighs = wmtk::simplex::cofaces_single_dimension_iterable( mesh(), simplex, mesh().top_simplex_type()); @@ -267,7 +268,8 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: if (mesh().top_simplex_type() == PrimitiveType::Triangle) { std::vector> cells; - for (const simplex::Simplex& cell : neighs) { + for (const Tuple& cell_tuple : neighs) { + const simplex::Simplex cell(mesh().top_simplex_type(), cell_tuple); cells.emplace_back(m_amips.get_raw_coordinates<3, 2>(cell, simplex)); } WMTKAMIPSProblem<6> problem(cells); @@ -287,7 +289,8 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: std::vector> cells; - for (simplex::Simplex cell : neighs) { + for (const Tuple& cell_tuple : neighs) { + simplex::Simplex cell(mesh().top_simplex_type(), cell_tuple); // auto vertices = mesh().orient_vertices(cell.tuple()); // int idx = -1; // for (int i = 0; i < 4; ++i) { @@ -399,7 +402,7 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: } else { auto accessor = mesh().create_accessor(m_coordinate_handle.as()); - const auto neighs = wmtk::simplex::cofaces_single_dimension_simplices( + auto neighs = wmtk::simplex::cofaces_single_dimension_iterable( mesh(), simplex, mesh().top_simplex_type()); @@ -407,7 +410,8 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: if (mesh().top_simplex_type() == PrimitiveType::Triangle) { std::vector> cells; - for (const simplex::Simplex& cell : neighs) { + for (const Tuple& cell_tuple : neighs) { + const simplex::Simplex cell(mesh().top_simplex_type(), cell_tuple); cells.emplace_back(m_amips.get_raw_coordinates<3, 2>(cell, simplex)); } WMTKAMIPSProblem<6> problem(cells); @@ -427,7 +431,9 @@ std::vector AMIPSOptimizationSmoothing::execute(const simplex: std::vector> cells; - for (simplex::Simplex cell : neighs) { + for (const Tuple& cell_tuple : neighs) { + simplex::Simplex cell(mesh().top_simplex_type(), cell_tuple); + if (!mesh().is_ccw(cell.tuple())) { // switch any local id but NOT the vertex cell = simplex::Simplex( diff --git a/src/wmtk/operations/EdgeCollapse.cpp b/src/wmtk/operations/EdgeCollapse.cpp index 78ad76851f..f5f455d533 100644 --- a/src/wmtk/operations/EdgeCollapse.cpp +++ b/src/wmtk/operations/EdgeCollapse.cpp @@ -8,6 +8,8 @@ #include #include "utils/multi_mesh_edge_collapse.hpp" +#include + namespace wmtk::operations { @@ -53,6 +55,12 @@ EdgeCollapse::EdgeCollapse(Mesh& m) custom_attribute_collector.execute_from_root(m); } +EdgeCollapse::EdgeCollapse(submesh::Embedding& m) + : EdgeCollapse(m.mesh()) +{ + m.set_collapse_strategies(*this); +} + std::vector EdgeCollapse::execute(const simplex::Simplex& simplex) { return utils::multi_mesh_edge_collapse_with_modified_simplices( diff --git a/src/wmtk/operations/EdgeCollapse.hpp b/src/wmtk/operations/EdgeCollapse.hpp index 4b92136f59..da409afe2c 100644 --- a/src/wmtk/operations/EdgeCollapse.hpp +++ b/src/wmtk/operations/EdgeCollapse.hpp @@ -3,12 +3,17 @@ #include "Operation.hpp" #include "attribute_new/CollapseNewAttributeStrategy.hpp" +namespace wmtk::submesh { +class Embedding; +} + namespace wmtk::operations { class EdgeCollapse : public Operation { public: // constructor for default factory pattern construction EdgeCollapse(Mesh& m); + EdgeCollapse(submesh::Embedding& m); PrimitiveType primitive_type() const override { return PrimitiveType::Edge; } diff --git a/src/wmtk/operations/EdgeSplit.cpp b/src/wmtk/operations/EdgeSplit.cpp index eae7e549c6..d81615cfb8 100644 --- a/src/wmtk/operations/EdgeSplit.cpp +++ b/src/wmtk/operations/EdgeSplit.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -53,6 +54,12 @@ EdgeSplit::EdgeSplit(Mesh& m) custom_attribute_collector.execute_from_root(m); } +EdgeSplit::EdgeSplit(submesh::Embedding& m) + : EdgeSplit(m.mesh()) +{ + m.set_split_strategies(*this); +} + std::vector EdgeSplit::execute(const simplex::Simplex& simplex) { return utils::multi_mesh_edge_split_with_modified_simplices( diff --git a/src/wmtk/operations/EdgeSplit.hpp b/src/wmtk/operations/EdgeSplit.hpp index 15a0e35b27..75d99a86a8 100644 --- a/src/wmtk/operations/EdgeSplit.hpp +++ b/src/wmtk/operations/EdgeSplit.hpp @@ -4,11 +4,16 @@ #include "attribute_new/SplitNewAttributeStrategy.hpp" +namespace wmtk::submesh { +class Embedding; +} + namespace wmtk::operations { class EdgeSplit : public Operation { public: EdgeSplit(Mesh& m); + EdgeSplit(submesh::Embedding& m); PrimitiveType primitive_type() const override { return PrimitiveType::Edge; } diff --git a/src/wmtk/operations/composite/ProjectOperation.cpp b/src/wmtk/operations/composite/ProjectOperation.cpp index daf81341fe..ddf4067c39 100644 --- a/src/wmtk/operations/composite/ProjectOperation.cpp +++ b/src/wmtk/operations/composite/ProjectOperation.cpp @@ -2,6 +2,9 @@ #include #include +#include +#include +#include #include @@ -18,8 +21,15 @@ ProjectOperation::ProjectOperation( ProjectOperation::ProjectOperation( std::shared_ptr main_op, const std::vector& mesh_constaint_pairs) - : AttributesUpdate(main_op->mesh()) - , m_main_op(main_op) + : ProjectOperation(main_op->mesh(), mesh_constaint_pairs) +{ + m_main_op = main_op; +} + +ProjectOperation::ProjectOperation( + Mesh& mesh, + const std::vector& mesh_constaint_pairs) + : AttributesUpdate(mesh) { for (auto& pair : mesh_constaint_pairs) { int64_t count = 0; @@ -84,37 +94,120 @@ ProjectOperation::ProjectOperation( } } +ProjectOperation::ProjectOperation( + std::shared_ptr main_op, + submesh::Embedding& emb, + const attribute::MeshAttributeHandle& pos_handle) + : ProjectOperation(emb, pos_handle) +{ + m_main_op = main_op; +} + +ProjectOperation::ProjectOperation( + submesh::Embedding& emb, + const attribute::MeshAttributeHandle& pos_handle) + : AttributesUpdate(emb.mesh()) + , m_embedding_pos_handle(pos_handle) +{ + const Mesh& m = emb.mesh(); + assert(&m == &pos_handle.mesh()); + // Wrapper for the position accessor that works for double and Rational. Probably not the most + // efficient code but good enough for what is required here + auto get_pos = [&pos_handle, &m](const Tuple& _t) -> Eigen::VectorXd { + return std::visit( + [&m, &_t](auto&& tah) noexcept -> Eigen::VectorXd { + using HandleType = typename std::decay_t; + using AttributeType = typename HandleType::Type; + + const auto accessor = m.create_const_accessor(tah); + if constexpr (std::is_same_v) { + return accessor.const_vector_attribute(_t); + } + if constexpr (std::is_same_v) { + return accessor.const_vector_attribute(_t).template cast(); + } + log_and_throw_error("Position attribute must be double or rational"); + }, + pos_handle.handle()); + }; + + for (const auto& sub_ptr : emb.get_child_meshes()) { + const submesh::SubMesh& sub = *sub_ptr; + const PrimitiveType pt = sub.top_simplex_type(); + + if (pt == PrimitiveType::Vertex) { + logger().info("Ignoring vertex submeshes in ProjectOperation"); + continue; + } + + int64_t count = 0; + int64_t index = 0; + + const std::vector facest = sub.get_all_id_simplex(pt); + + const int64_t dim = int64_t(pt) + 1; + + Eigen::MatrixXd vertices(dim * facest.size(), pos_handle.dimension()); + Eigen::MatrixXi faces(facest.size(), dim); + + for (const simplex::IdSimplex& cell : facest) { + const auto tmp = + faces_single_dimension_tuples(m, m.get_simplex(cell), PrimitiveType::Vertex); + + assert(tmp.size() == dim); + for (int64_t j = 0; j < tmp.size(); ++j) { + Eigen::VectorXd p = get_pos(tmp[j]); + faces(index, j) = count; + vertices.row(dim * index + j) = p; + + ++count; + } + ++index; + } + + auto bvh = std::make_shared(); + bvh->init(vertices, faces, 1e-10); + + m_submesh_bvh.emplace_back(sub_ptr, bvh); + } +} + std::vector ProjectOperation::execute(const simplex::Simplex& simplex) { - // mesh has to be the same as the main_op mesh - assert(&m_main_op->mesh() == &mesh()); - const auto main_simplices = (*m_main_op)(simplex); - if (main_simplices.empty()) return {}; - assert(main_simplices.size() == 1); + std::vector main_simplices; + if (m_main_op) { + // mesh has to be the same as the main_op mesh + assert(&m_main_op->mesh() == &mesh()); + main_simplices = (*m_main_op)(simplex); + if (main_simplices.empty()) return {}; + assert(main_simplices.size() == 1); + } else { + main_simplices.emplace_back(simplex); + } const auto main_tup = main_simplices.front().tuple(); - for (auto& pair : m_bvh) { + for (auto& [pos_attribute, bvh] : m_bvh) { const std::vector mapped_tuples_after = - mesh().map_tuples(pair.first.mesh(), primitive_type(), {main_tup}); + mesh().map_tuples(pos_attribute.mesh(), primitive_type(), {main_tup}); if (mapped_tuples_after.empty()) continue; - if (pair.first.holds()) { + if (pos_attribute.holds()) { wmtk::attribute::Accessor accessor = - pair.first.mesh().create_accessor(pair.first.as()); + pos_attribute.mesh().create_accessor(pos_attribute.as()); for (const auto& t : mapped_tuples_after) { auto p = accessor.vector_attribute(t); SimpleBVH::VectorMax3d nearest_point; double sq_dist; - pair.second->nearest_facet(p, nearest_point, sq_dist); + bvh->nearest_facet(p, nearest_point, sq_dist); p = nearest_point; } } else { - assert((pair.first.holds())); + assert((pos_attribute.holds())); wmtk::attribute::Accessor accessor = - pair.first.mesh().create_accessor(pair.first.as()); + pos_attribute.mesh().create_accessor(pos_attribute.as()); for (const auto& t : mapped_tuples_after) { auto p_map = accessor.vector_attribute(t); @@ -123,16 +216,16 @@ std::vector ProjectOperation::execute(const simplex::Simplex& const Eigen::Vector3d p = p_map.cast(); SimpleBVH::VectorMax3d nearest_point; double sq_dist; - pair.second->nearest_facet(p, nearest_point, sq_dist); - for (int64_t d = 0; d < pair.first.dimension(); ++d) { + bvh->nearest_facet(p, nearest_point, sq_dist); + for (int64_t d = 0; d < pos_attribute.dimension(); ++d) { p_map(d) = Rational(nearest_point[d], true); } } else if (p_map.rows() == 2) { const Eigen::Vector2d p = p_map.cast(); SimpleBVH::VectorMax3d nearest_point; double sq_dist; - pair.second->nearest_facet(p, nearest_point, sq_dist); - for (int64_t d = 0; d < pair.first.dimension(); ++d) { + bvh->nearest_facet(p, nearest_point, sq_dist); + for (int64_t d = 0; d < pos_attribute.dimension(); ++d) { p_map(d) = Rational(nearest_point[d], true); } } else { @@ -142,6 +235,50 @@ std::vector ProjectOperation::execute(const simplex::Simplex& } } + for (auto& [sub_ptr, bvh] : m_submesh_bvh) { + const submesh::SubMesh& sub = *sub_ptr; + if (!sub.contains(main_tup, primitive_type())) { + continue; + } + + if (m_embedding_pos_handle.holds()) { + wmtk::attribute::Accessor accessor = + mesh().create_accessor(m_embedding_pos_handle.as()); + + auto p = accessor.vector_attribute(main_tup); + SimpleBVH::VectorMax3d nearest_point; + double sq_dist; + bvh->nearest_facet(p, nearest_point, sq_dist); + p = nearest_point; + } else { + assert((m_embedding_pos_handle.holds())); + wmtk::attribute::Accessor accessor = + mesh().create_accessor(m_embedding_pos_handle.as()); + + auto p_map = accessor.vector_attribute(main_tup); + + if (p_map.rows() == 3) { + const Eigen::Vector3d p = p_map.cast(); + SimpleBVH::VectorMax3d nearest_point; + double sq_dist; + bvh->nearest_facet(p, nearest_point, sq_dist); + for (int64_t d = 0; d < m_embedding_pos_handle.dimension(); ++d) { + p_map(d) = Rational(nearest_point[d], true); + } + } else if (p_map.rows() == 2) { + const Eigen::Vector2d p = p_map.cast(); + SimpleBVH::VectorMax3d nearest_point; + double sq_dist; + bvh->nearest_facet(p, nearest_point, sq_dist); + for (int64_t d = 0; d < m_embedding_pos_handle.dimension(); ++d) { + p_map(d) = Rational(nearest_point[d], true); + } + } else { + throw std::runtime_error("wrong vector dimension"); + } + } + } + return main_simplices; } diff --git a/src/wmtk/operations/composite/ProjectOperation.hpp b/src/wmtk/operations/composite/ProjectOperation.hpp index b6df673d38..bc5d880531 100644 --- a/src/wmtk/operations/composite/ProjectOperation.hpp +++ b/src/wmtk/operations/composite/ProjectOperation.hpp @@ -8,10 +8,16 @@ namespace SimpleBVH { class BVH; } +namespace wmtk::submesh { +class Embedding; +class SubMesh; +} // namespace wmtk::submesh + namespace wmtk::operations::composite { class ProjectOperation : public AttributesUpdate { public: + // first: construction, second: query using MeshConstrainPair = std::pair; @@ -24,15 +30,28 @@ class ProjectOperation : public AttributesUpdate std::shared_ptr main_op, const std::vector& mesh_constaint_pairs); + ProjectOperation(Mesh& mesh, const std::vector& mesh_constaint_pairs); + + ProjectOperation( + std::shared_ptr main_op, + submesh::Embedding& emb, + const attribute::MeshAttributeHandle& pos_handle); + + ProjectOperation(submesh::Embedding& emb, const attribute::MeshAttributeHandle& pos_handle); + std::vector execute(const simplex::Simplex& simplex) override; - PrimitiveType primitive_type() const override { return m_main_op->primitive_type(); } + PrimitiveType primitive_type() const override { return PrimitiveType::Vertex; } private: using BVHConstrainPair = std::pair>; + using BVHSubMeshConstraintPair = + std::pair, std::shared_ptr>; - const std::shared_ptr m_main_op; + std::shared_ptr m_main_op; std::vector m_bvh; + std::vector m_submesh_bvh; + attribute::MeshAttributeHandle m_embedding_pos_handle; }; } // namespace wmtk::operations::composite diff --git a/src/wmtk/operations/composite/TriEdgeSwap.cpp b/src/wmtk/operations/composite/TriEdgeSwap.cpp index 99935fb082..f49e4367a5 100644 --- a/src/wmtk/operations/composite/TriEdgeSwap.cpp +++ b/src/wmtk/operations/composite/TriEdgeSwap.cpp @@ -1,6 +1,7 @@ #include "TriEdgeSwap.hpp" #include +#include namespace wmtk::operations::composite { @@ -10,6 +11,12 @@ TriEdgeSwap::TriEdgeSwap(Mesh& m) , m_collapse(m) {} +TriEdgeSwap::TriEdgeSwap(submesh::Embedding& m) + : Operation(m.mesh()) + , m_split(m) + , m_collapse(m) +{} + std::vector TriEdgeSwap::execute(const simplex::Simplex& simplex) { diff --git a/src/wmtk/operations/composite/TriEdgeSwap.hpp b/src/wmtk/operations/composite/TriEdgeSwap.hpp index 2e35b040f4..9c77ca99fc 100644 --- a/src/wmtk/operations/composite/TriEdgeSwap.hpp +++ b/src/wmtk/operations/composite/TriEdgeSwap.hpp @@ -4,6 +4,10 @@ #include #include +namespace wmtk::submesh { +class Embedding; +} + namespace wmtk::operations::composite { /** * Performs an edge swap, implemented as a combination of swap and collapse. @@ -41,6 +45,7 @@ class TriEdgeSwap : public Operation { public: TriEdgeSwap(Mesh& m); + TriEdgeSwap(submesh::Embedding& m); PrimitiveType primitive_type() const override { return PrimitiveType::Edge; } diff --git a/src/wmtk/simplex/IdSimplexCollection.cpp b/src/wmtk/simplex/IdSimplexCollection.cpp index e5ace4a375..fea82a3d48 100644 --- a/src/wmtk/simplex/IdSimplexCollection.cpp +++ b/src/wmtk/simplex/IdSimplexCollection.cpp @@ -40,6 +40,19 @@ std::vector IdSimplexCollection::simplex_vector_tuples(PrimitiveType ptyp return tuples; } +std::vector IdSimplexCollection::simplex_vector_tuples() const +{ + std::vector tuples; + tuples.reserve(m_simplices.size()); + + // add simplices to the vector + for (const IdSimplex& s : m_simplices) { + tuples.emplace_back(m_mesh.get_tuple_from_id_simplex(s)); + } + + return tuples; +} + void IdSimplexCollection::add(const IdSimplex& simplex) { m_simplices.push_back(simplex); diff --git a/src/wmtk/simplex/IdSimplexCollection.hpp b/src/wmtk/simplex/IdSimplexCollection.hpp index d232608cb9..d16f25d3d5 100644 --- a/src/wmtk/simplex/IdSimplexCollection.hpp +++ b/src/wmtk/simplex/IdSimplexCollection.hpp @@ -30,6 +30,7 @@ class IdSimplexCollection * @brief Return vector of all simplices of the requested type, as tuples */ std::vector simplex_vector_tuples(PrimitiveType ptype) const; + std::vector simplex_vector_tuples() const; /** * @brief Add simplex to the collection. diff --git a/src/wmtk/submesh/CMakeLists.txt b/src/wmtk/submesh/CMakeLists.txt new file mode 100644 index 0000000000..db8bce5608 --- /dev/null +++ b/src/wmtk/submesh/CMakeLists.txt @@ -0,0 +1,12 @@ +set(SRC_FILES + Embedding.hpp + Embedding.cpp + SubMesh.hpp + SubMesh.cpp + + utils/write.hpp + utils/write.cpp + utils/submesh_from_multimesh.hpp + utils/submesh_from_multimesh.cpp +) +target_sources(wildmeshing_toolkit PRIVATE ${SRC_FILES}) \ No newline at end of file diff --git a/src/wmtk/submesh/Embedding.cpp b/src/wmtk/submesh/Embedding.cpp new file mode 100644 index 0000000000..a94a3b791c --- /dev/null +++ b/src/wmtk/submesh/Embedding.cpp @@ -0,0 +1,348 @@ +#include "Embedding.hpp" + +#include +#include +#include +#include +#include + +#include "SubMesh.hpp" + +namespace wmtk::submesh { + +Embedding::Embedding(const std::shared_ptr& mesh) + : m_mesh(mesh) +{ + m_mesh->m_embedding = this; + + m_tag_attribute_name[PrimitiveType::Vertex] = "WMTK_submesh_tag_v"; + m_tag_attribute_name[PrimitiveType::Edge] = "WMTK_submesh_tag_e"; + m_tag_attribute_name[PrimitiveType::Triangle] = "WMTK_submesh_tag_f"; + m_tag_attribute_name[PrimitiveType::Tetrahedron] = "WMTK_submesh_tag_t"; + + Mesh& m = *m_mesh; + + // register tag attributes + for (const PrimitiveType& pt : utils::primitive_below(m.top_simplex_type())) { + if (m.has_attribute(m_tag_attribute_name[pt], pt)) { + log_and_throw_error( + "Cannot create embedding. Mesh already has an attribute with name {}", + m_tag_attribute_name[pt]); + } + + m_tag_handle[pt] = m.register_attribute_typed(m_tag_attribute_name[pt], pt, 1); + + + attribute::MeshAttributeHandle h(m, m_tag_handle.at(pt)); + + m_split_new[pt] = std::make_shared>(h); + auto& split_new_strat = *(m_split_new[pt]); + + split_new_strat.set_strategy(operations::SplitBasicStrategy::Copy); + split_new_strat.set_rib_strategy(operations::SplitRibBasicStrategy::None); + + m_collapse_new[pt] = std::make_shared>(h); + auto& collapse_new_strat = *(m_collapse_new[pt]); + + auto collapse_new_func = [](const VectorX& a, + const VectorX& b, + const std::bitset<2>&) -> VectorX { + VectorX r(a.rows()); + + assert(a.rows() == b.rows()); + assert(a.rows() == r.rows()); + + for (int64_t i = 0; i < a.rows(); ++i) { + r[i] = a[i] | b[i]; + } + + return r; + }; + collapse_new_strat.set_strategy(collapse_new_func); + } + + + m_substructure_predicate = [this](const simplex::Simplex& s) -> bool { + return simplex_is_in_submesh(s); + }; + + auto update_tag_func = [this]( + const Eigen::MatrixX& P, + const std::vector& tuples) -> Eigen::VectorX { + // transfer from vertices (P.cols()) to top_simplex + assert(P.rows() == 1); // rows --> attribute dimension + // cols --> number of input simplices (vertices) + + const simplex::Simplex cell(m_mesh->top_simplex_type(), tuples[0]); + + assert(cell.primitive_type() != PrimitiveType::Vertex); + + // transfer from cell to facets + int64_t cell_tag; + { + auto tag_cell_acc = tag_accessor(cell.primitive_type()); + cell_tag = tag_cell_acc.const_scalar_attribute(cell); + + auto tag_facet_acc = tag_accessor(cell.primitive_type() - 1); + const auto facets = + simplex::faces_single_dimension_tuples(*m_mesh, cell, cell.primitive_type() - 1); + + for (const Tuple& f : facets) { + tag_facet_acc.scalar_attribute(f) |= cell_tag; + } + } + + if (cell.primitive_type() != PrimitiveType::Edge) { + // cell is triangle or tet + for (const PrimitiveType& pt : + utils::primitive_range(cell.primitive_type() - 1, PrimitiveType::Edge)) { + auto s_acc = tag_accessor(pt); + auto f_acc = tag_accessor(pt - 1); + const auto simplices = simplex::faces_single_dimension(*m_mesh, cell, pt); + for (const simplex::Simplex& s : simplices) { + const auto faces = simplex::faces_single_dimension_tuples(*m_mesh, s, pt - 1); + for (const Tuple& f : faces) { + f_acc.scalar_attribute(f) |= s_acc.const_scalar_attribute(s); + } + } + } + } + + return Eigen::VectorX::Constant(1, cell_tag); + }; + + attribute::MeshAttributeHandle h_v(m, m_tag_handle.at(PrimitiveType::Vertex)); + attribute::MeshAttributeHandle h_c(m, m_tag_handle.at(m.top_simplex_type())); + + m_transfer = + std::make_shared>( + h_c, + h_v, + update_tag_func); +} + +std::shared_ptr Embedding::add_submesh() +{ + if (m_submeshes.size() == 63) { + log_and_throw_error("An embedding can only hold up to 63 submeshes"); + } + + std::shared_ptr sub = std::make_shared(*this, m_submeshes.size()); + m_submeshes.emplace_back(sub); + + return sub; +} + +Mesh& Embedding::mesh() +{ + return *m_mesh; +} + +const Mesh& Embedding::mesh() const +{ + return *m_mesh; +} + +attribute::TypedAttributeHandle& Embedding::tag_handle(const PrimitiveType pt) +{ + return m_tag_handle.at(pt); +} + +attribute::Accessor Embedding::tag_accessor(const PrimitiveType pt) +{ + return mesh().create_accessor(m_tag_handle.at(pt)); +} + +const attribute::Accessor Embedding::tag_accessor(const PrimitiveType pt) const +{ + return mesh().create_const_accessor(m_tag_handle.at(pt)); +} + +std::vector Embedding::get_all(PrimitiveType type) const +{ + return mesh().get_all(type); +} + +std::vector Embedding::get_all_id_simplex(PrimitiveType type) const +{ + return mesh().get_all_id_simplex(type); +} + +int64_t Embedding::top_cell_dimension() const +{ + return mesh().top_cell_dimension(); +} + +MeshType Embedding::mesh_type() const +{ + return MeshType::Embedding; +} + +Tuple Embedding::switch_tuple(const Tuple& tuple, PrimitiveType type) const +{ + return mesh().switch_tuple(tuple, type); +} + +bool Embedding::is_boundary(PrimitiveType pt, const Tuple& tuple) const +{ + return mesh().is_boundary(pt, tuple); +} + +int64_t Embedding::id(const simplex::Simplex& s) const +{ + return mesh().id(s); +} + +int64_t Embedding::id(const Tuple& tuple, PrimitiveType pt) const +{ + return mesh().id(tuple, pt); +} + +std::vector> Embedding::get_child_meshes() const +{ + return m_submeshes; +} + +bool Embedding::has_child_mesh() const +{ + return !m_submeshes.empty(); +} + +bool Embedding::simplex_is_in_submesh(const simplex::Simplex& s) const +{ + const auto acc = tag_accessor(s.primitive_type()); + return acc.const_scalar_attribute(s.tuple()) > 0; +} + +void Embedding::set_split_strategies(operations::EdgeSplit& split) const +{ + for (const auto& [pt, strat] : m_split_new) { + assert(m_mesh->validate_handle(m_tag_handle.at(pt))); + attribute::MeshAttributeHandle h(*m_mesh, m_tag_handle.at(pt)); + split.set_new_attribute_strategy(h, strat); + } + + split.add_transfer_strategy(m_transfer); +} + +void Embedding::set_collapse_strategies(operations::EdgeCollapse& collapse) const +{ + for (const auto& [pt, strat] : m_collapse_new) { + attribute::MeshAttributeHandle h(*m_mesh, m_tag_handle.at(pt)); + collapse.set_new_attribute_strategy(h, strat); + } +} + +std::function Embedding::substructure_predicate() const +{ + return m_substructure_predicate; +} + +void Embedding::update_tag_attribute_handles() +{ + Mesh& m = *m_mesh; + + for (const PrimitiveType& pt : utils::primitive_below(m.top_simplex_type())) { + if (!m.has_attribute(m_tag_attribute_name[pt], pt)) { + log_and_throw_error( + "Cannot update handles. Mesh already no attribute with name {}", + m_tag_attribute_name[pt]); + } + const auto& name = m_tag_attribute_name[pt]; + m_tag_handle[pt] = m.get_attribute_handle_typed(name, pt); + + + attribute::MeshAttributeHandle h(m, m_tag_handle.at(pt)); + + m_split_new[pt] = std::make_shared>(h); + auto& split_new_strat = *(m_split_new[pt]); + + split_new_strat.set_strategy(operations::SplitBasicStrategy::Copy); + split_new_strat.set_rib_strategy(operations::SplitRibBasicStrategy::None); + + m_collapse_new[pt] = std::make_shared>(h); + auto& collapse_new_strat = *(m_collapse_new[pt]); + + auto collapse_new_func = [](const VectorX& a, + const VectorX& b, + const std::bitset<2>&) -> VectorX { + VectorX r(a.rows()); + + assert(a.rows() == b.rows()); + assert(a.rows() == r.rows()); + + for (int64_t i = 0; i < a.rows(); ++i) { + r[i] = a[i] | b[i]; + } + + return r; + }; + collapse_new_strat.set_strategy(collapse_new_func); + } + + + m_substructure_predicate = [this](const simplex::Simplex& s) -> bool { + return simplex_is_in_submesh(s); + }; + + auto update_tag_func = [this]( + const Eigen::MatrixX& P, + const std::vector& tuples) -> Eigen::VectorX { + // transfer from vertices (P.cols()) to top_simplex + assert(P.rows() == 1); // rows --> attribute dimension + // cols --> number of input simplices (vertices) + + const simplex::Simplex cell(m_mesh->top_simplex_type(), tuples[0]); + + assert(cell.primitive_type() != PrimitiveType::Vertex); + + // transfer from cell to facets + int64_t cell_tag; + { + auto tag_cell_acc = tag_accessor(cell.primitive_type()); + cell_tag = tag_cell_acc.const_scalar_attribute(cell); + + auto tag_facet_acc = tag_accessor(cell.primitive_type() - 1); + const auto facets = + simplex::faces_single_dimension_tuples(*m_mesh, cell, cell.primitive_type() - 1); + + for (const Tuple& f : facets) { + tag_facet_acc.scalar_attribute(f) |= cell_tag; + } + } + + if (cell.primitive_type() != PrimitiveType::Edge) { + // cell is triangle or tet + for (const PrimitiveType& pt : + utils::primitive_range(cell.primitive_type() - 1, PrimitiveType::Edge)) { + auto s_acc = tag_accessor(pt); + auto f_acc = tag_accessor(pt - 1); + const auto simplices = simplex::faces_single_dimension(*m_mesh, cell, pt); + for (const simplex::Simplex& s : simplices) { + const auto faces = simplex::faces_single_dimension_tuples(*m_mesh, s, pt - 1); + for (const Tuple& f : faces) { + f_acc.scalar_attribute(f) |= s_acc.const_scalar_attribute(s); + } + } + } + } + + return Eigen::VectorX::Constant(1, cell_tag); + }; + + attribute::MeshAttributeHandle h_v(m, m_tag_handle.at(PrimitiveType::Vertex)); + attribute::MeshAttributeHandle h_c(m, m_tag_handle.at(m.top_simplex_type())); + + m_transfer = + std::make_shared>( + h_c, + h_v, + update_tag_func); + + for (const auto& [pt, h] : m_tag_handle) { + assert(m.validate_handle(h)); + } +} + + +} // namespace wmtk::submesh diff --git a/src/wmtk/submesh/Embedding.hpp b/src/wmtk/submesh/Embedding.hpp new file mode 100644 index 0000000000..52f6a79eaa --- /dev/null +++ b/src/wmtk/submesh/Embedding.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace wmtk::operations { +class EdgeSplit; +class EdgeCollapse; +} // namespace wmtk::operations + +namespace wmtk::submesh { + +class SubMesh; + +/** + * The embedding is a wrapper for the embedding mesh for the submeshes. It contains a pointer to the + * mesh, the attribute for the submesh tags, and a factory for SubMeshes. + */ +class Embedding : public std::enable_shared_from_this, public MeshBase +{ +public: + Embedding(const std::shared_ptr& mesh); + Embedding(Embedding&) = delete; + + std::shared_ptr add_submesh(); + + Mesh& mesh(); + const Mesh& mesh() const; + + attribute::TypedAttributeHandle& tag_handle(const PrimitiveType pt); + attribute::Accessor tag_accessor(const PrimitiveType pt); + const attribute::Accessor tag_accessor(const PrimitiveType pt) const; + + std::vector get_all(PrimitiveType type) const override; + std::vector get_all_id_simplex(PrimitiveType type) const override; + + int64_t top_cell_dimension() const override; + + MeshType mesh_type() const override; + + Tuple switch_tuple(const Tuple& tuple, PrimitiveType type) const override; + bool is_boundary(PrimitiveType, const Tuple& tuple) const override; + int64_t id(const simplex::Simplex& s) const override; + int64_t id(const Tuple& tuple, PrimitiveType pt) const override; + + std::vector> get_child_meshes() const; + + bool has_child_mesh() const; + + bool simplex_is_in_submesh(const simplex::Simplex& s) const; + + void set_split_strategies(operations::EdgeSplit& split) const; + + void set_collapse_strategies(operations::EdgeCollapse& collapse) const; + + std::function substructure_predicate() const; + + /** + * @brief Update the tag attribute handles. + * + * This function must be called after removing attributes or deregistering child meshes. + */ + void update_tag_attribute_handles(); + +private: + std::shared_ptr m_mesh; + + std::map m_tag_attribute_name; + std::map> m_tag_handle; + + std::vector> m_submeshes; + + std::map>> + m_split_new; + std::map>> + m_collapse_new; + + std::shared_ptr> m_transfer; + + std::function m_substructure_predicate; +}; + +} // namespace wmtk::submesh \ No newline at end of file diff --git a/src/wmtk/submesh/SubMesh.cpp b/src/wmtk/submesh/SubMesh.cpp new file mode 100644 index 0000000000..04c898d9d1 --- /dev/null +++ b/src/wmtk/submesh/SubMesh.cpp @@ -0,0 +1,321 @@ +#include "SubMesh.hpp" + +#include "Embedding.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace wmtk::submesh { + +SubMesh::SubMesh(Embedding& embedding, int64_t submesh_id) + : m_embedding(embedding) + , m_submesh_id(submesh_id) +{} + +Mesh& SubMesh::mesh() +{ + return m_embedding.mesh(); +} + +const Mesh& SubMesh::mesh() const +{ + return m_embedding.mesh(); +} + +void SubMesh::add_simplex(const Tuple& tuple, PrimitiveType pt) +{ + const int64_t pt_dim = get_primitive_type_id(pt); + if (m_top_cell_dimension < pt_dim) { + logger().trace("Top cell dimension changed from {} to {}", m_top_cell_dimension, pt_dim); + m_top_cell_dimension = pt_dim; + } + + auto acc = m_embedding.tag_accessor(pt); + acc.scalar_attribute(tuple) |= (int64_t)1 << m_submesh_id; + + const simplex::Simplex s(mesh(), pt, tuple); + + auto faces = simplex::faces(mesh(), s); + + for (const simplex::Simplex& f : faces) { + auto a = m_embedding.tag_accessor(f.primitive_type()); + a.scalar_attribute(f.tuple()) |= (int64_t)1 << m_submesh_id; + } +} + +void SubMesh::add_simplex(const simplex::Simplex& s) +{ + add_simplex(s.tuple(), s.primitive_type()); +} + +void SubMesh::add_simplex(const simplex::IdSimplex& simplex) +{ + const PrimitiveType& pt = simplex.primitive_type(); + const Tuple t = mesh().get_tuple_from_id_simplex(simplex); + add_simplex(t, pt); +} + +std::vector SubMesh::get_all(const PrimitiveType pt) const +{ + const auto all = mesh().get_all(pt); + std::vector sub; + for (const Tuple& t : all) { + if (contains(t, pt)) { + sub.emplace_back(t); + } + } + return sub; +} + +std::vector SubMesh::get_all_id_simplex(PrimitiveType pt) const +{ + const auto all = mesh().get_all_id_simplex(pt); + std::vector sub; + for (const simplex::IdSimplex& s : all) { + if (contains(s)) { + sub.emplace_back(s); + } + } + return sub; +} + +PrimitiveType SubMesh::top_simplex_type(const Tuple& tuple) const +{ + const Mesh& m = mesh(); + + for (const PrimitiveType& pt : utils::primitive_below(m.top_simplex_type())) { + if (contains(tuple, pt)) { + return pt; + } + } + + log_and_throw_error("No simplex of the tuple contains the submesh tag."); +} + +MeshType SubMesh::mesh_type() const +{ + return MeshType::SubMesh; +} + +PrimitiveType SubMesh::top_simplex_type() const +{ + // const Mesh& m = mesh(); + // + // for (const PrimitiveType& pt : utils::primitive_below(m.top_simplex_type())) { + // const auto tuples = m.get_all(pt); + // for (const Tuple& t : tuples) { + // if (contains(t, pt)) { + // return pt; + // } + // } + // } + // + // log_and_throw_error("No simplex of the tuple contains the submesh tag."); + + if (m_top_cell_dimension < 0) { + log_and_throw_error("No simplex of the tuple contains the submesh tag."); + } + + return get_primitive_type_from_id(m_top_cell_dimension); +} + +int64_t SubMesh::top_cell_dimension() const +{ + assert(m_top_cell_dimension >= 0); + assert(m_top_cell_dimension < 4); + return m_top_cell_dimension; +} + +Tuple SubMesh::switch_tuple(const Tuple& tuple, PrimitiveType pt) const +{ + const int8_t pt_id = get_primitive_type_id(pt); + const int8_t max_pt_id = get_primitive_type_id(top_simplex_type(tuple)); + + if (pt_id > max_pt_id) { + // invalid switch + log_and_throw_error("Required PrimitiveType switch does not exist in submesh."); + } + if (pt_id == max_pt_id) { + // global switch + const PrimitiveType pt_face = get_primitive_type_from_id(pt_id - 1); + const simplex::Simplex s_face(mesh(), pt_face, tuple); + + Tuple other; + int64_t other_counter = 0; + for (const Tuple& t : simplex::cofaces_single_dimension_iterable(mesh(), s_face, pt)) { + if (contains(t, pt)) { + ++other_counter; + if (other_counter == 2) { + other = t; + } + } + } + + if (other_counter != 2) { + log_and_throw_error( + "SubMesh `switch_tuple` cannot be used on non-manifold or boundary simplices."); + } + return other; + } + + // local switch + return local_switch_tuple(tuple, pt); +} + +Tuple SubMesh::switch_tuples( + const Tuple& tuple, + const std::initializer_list& op_sequence) const +{ + Tuple r = tuple; + for (const PrimitiveType& primitive : op_sequence) { + r = switch_tuple(r, primitive); + } + return r; +} + +// std::vector SubMesh::switch_tuple_vector(const Tuple& tuple, PrimitiveType pt) const +//{ +// const int8_t pt_id = get_primitive_type_id(pt); +// const int8_t max_pt_id = get_primitive_type_id(top_simplex_type(tuple)); +// +// if (pt_id > max_pt_id) { +// log_and_throw_error("Required PrimitiveType switch does not exist in submesh."); +// } +// if (pt_id < max_pt_id) { +// // log_and_throw_error("Cannot perform global switches for vertex PrimitiveType. Use " +// // "`switch_tuple` instead of `switch_tuple_vector`."); +// return {local_switch_tuple(tuple, pt)}; +// } +// +// assert(pt_id <= max_pt_id); +// +// const PrimitiveType pt_face = get_primitive_type_from_id(pt_id - 1); +// +// const simplex::Simplex s_face(mesh(), pt_face, tuple); +// +// std::vector neighs; +// neighs.reserve(2); +// for (const Tuple& t : simplex::cofaces_single_dimension_iterable(mesh(), s_face, pt)) { +// if (contains(t, pt)) { +// neighs.emplace_back(t); +// } +// } +// +// assert(!neighs.empty()); +// assert(neighs[0] == tuple); +// +// return neighs; +//} + +bool SubMesh::is_boundary(PrimitiveType pt, const Tuple& tuple) const +{ + if (!contains(tuple, pt)) { + log_and_throw_error("Cannot check for boundary if simplex is not contained in submesh"); + } + + const simplex::Simplex s(pt, tuple); + + const PrimitiveType top_pt = top_simplex_type(); + if (top_pt == PrimitiveType::Vertex) { + return true; + } + const PrimitiveType face_pt = get_primitive_type_from_id(top_cell_dimension() - 1); + + const auto neighbors = simplex::neighbors_single_dimension(mesh(), s, face_pt); + + // check if any incident facet has less than two neighbors + int64_t n_neighbors = 0; + for (const simplex::Simplex& face_simplex : neighbors) { + if (!contains(face_simplex)) { + continue; + } + ++n_neighbors; + + int64_t cell_counter = 0; + for (const Tuple& cell_tuple : + simplex::cofaces_single_dimension_iterable(mesh(), face_simplex, top_pt)) { + if (contains(cell_tuple, top_pt)) { + ++cell_counter; + } + } + if (cell_counter < 2) { + return true; + } + } + + if (n_neighbors == 0) { + // this simplex has no cell incident and is therefore considered boundary + return true; + } + + return false; +} + +bool SubMesh::is_boundary(const Tuple& tuple, PrimitiveType pt) const +{ + return is_boundary(pt, tuple); +} + +bool SubMesh::contains(const Tuple& tuple, PrimitiveType pt) const +{ + const auto acc = m_embedding.tag_accessor(pt); + return (acc.const_scalar_attribute(tuple) & ((int64_t)1 << m_submesh_id)) != 0; +} + +bool SubMesh::contains(const simplex::IdSimplex& s) const +{ + const auto acc = m_embedding.tag_accessor(s.primitive_type()); + return (acc.const_scalar_attribute(s) & ((int64_t)1 << m_submesh_id)) != 0; +} + +bool SubMesh::contains(const simplex::Simplex& s) const +{ + return contains(s.tuple(), s.primitive_type()); +} + +int64_t SubMesh::id(const Tuple& tuple, PrimitiveType pt) const +{ + return mesh().id(tuple, pt); +} + +int64_t SubMesh::id(const simplex::Simplex& s) const +{ + return id(s.tuple(), s.primitive_type()); +} + +Tuple SubMesh::local_switch_tuple(const Tuple& tuple, PrimitiveType pt) const +{ + return mesh().switch_tuple(tuple, pt); +} + +template +void SubMesh::add_from_tag_attribute( + const attribute::TypedAttributeHandle& tag_attribute, + const T tag_value) +{ + const attribute::Accessor acc = mesh().create_const_accessor(tag_attribute); + const PrimitiveType pt = tag_attribute.primitive_type(); + + const auto tuples = mesh().get_all(pt); + for (const Tuple& t : tuples) { + if (acc.const_scalar_attribute(t) == tag_value) { + add_simplex(t, pt); + } + } +} + +template void SubMesh::add_from_tag_attribute( + const attribute::TypedAttributeHandle&, + const int64_t); +template void SubMesh::add_from_tag_attribute( + const attribute::TypedAttributeHandle&, + const char); + +} // namespace wmtk::submesh diff --git a/src/wmtk/submesh/SubMesh.hpp b/src/wmtk/submesh/SubMesh.hpp new file mode 100644 index 0000000000..f280640b5f --- /dev/null +++ b/src/wmtk/submesh/SubMesh.hpp @@ -0,0 +1,126 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace wmtk { +namespace simplex { +class IdSimplex; +class Simplex; +} // namespace simplex +} // namespace wmtk + +namespace wmtk::submesh { + +class Embedding; + +/** + * Embedding m; + * SubMesh s = m.add_sub_mesh(); + * s.add_simplex(IdSimplex); + * s.add_from_tag(tag_handle, tag_value); + */ +class SubMesh : public std::enable_shared_from_this, public MeshBase +{ +public: + SubMesh(Embedding& embedding, int64_t submesh_id); + SubMesh() = delete; + SubMesh(const SubMesh&) = delete; + SubMesh& operator=(const SubMesh&) = delete; + SubMesh(SubMesh&&) = delete; + SubMesh& operator=(SubMesh&&) = delete; + ~SubMesh() override = default; + + Mesh& mesh(); + const Mesh& mesh() const; + + /** + * @brief Adds a simpex to the submesh. + * + * The top simplex type of the mesh is updated if the added simplex is of higher dimension than + * any other simplex of the submesh. + */ + void add_simplex(const Tuple& tuple, PrimitiveType pt); + void add_simplex(const simplex::Simplex& s); + void add_simplex(const simplex::IdSimplex& simplex); + + std::vector get_all(const PrimitiveType pt) const override; + std::vector get_all_id_simplex(PrimitiveType type) const override; + + /** + * @brief Get the maximum primitive type that has a tag for a given tuple. + */ + PrimitiveType top_simplex_type(const Tuple& tuple) const; + + MeshType mesh_type() const override; + + /** + * @brief Get the maximum primitive type that has a tag in the entire mesh. + */ + PrimitiveType top_simplex_type() const; + + int64_t top_cell_dimension() const override; + + /** + * @brief Can only perform local switches! + * + * Throws if `pt` is larger or equal to top_simplex_type(tuple) + */ + Tuple switch_tuple(const Tuple& tuple, PrimitiveType pt) const override; + + Tuple switch_tuples(const Tuple& tuple, const std::initializer_list& op_sequence) + const; + + /** + * This function does not make sense. The proper way to navigate is using + * cofaces_single_dimension. + */ + // std::vector switch_tuple_vector(const Tuple& tuple, PrimitiveType pt) const; + + /** + * @brief Is the given simplex on the boundary? + * + * This check is more complex than the one of the mesh itself but follows a similar idea. First, + * the top simplex type in the open star is determined. Next, we check if any facet in the star + * has less than two top simplices incident. + * + * Note that the behavior might be unexpected for non-homogenuous or non-manifold simplicial + * complexes! + */ + bool is_boundary(PrimitiveType pt, const Tuple& tuple) const override; + bool is_boundary(const Tuple& tuple, PrimitiveType pt) const; + + // This is going to be some ugly recursive stuff I guess... + bool is_manifold(PrimitiveType pt, const Tuple& tuple) const; + + // Check if sub mesh contains the simplex + bool contains(const Tuple& tuple, PrimitiveType pt) const; + bool contains(const simplex::IdSimplex& s) const; + bool contains(const simplex::Simplex& s) const; + + template + void add_from_tag_attribute( + const attribute::TypedAttributeHandle& tag_attribute, + const T tag_value); + + /** + * Wrapping the simplex id function of Mesh. + */ + int64_t id(const Tuple& tuple, PrimitiveType pt) const override; + int64_t id(const simplex::Simplex& s) const override; + +private: + Embedding& m_embedding; + const int64_t m_submesh_id; + + int64_t m_top_cell_dimension = -1; + +private: + Tuple local_switch_tuple(const Tuple& tuple, PrimitiveType pt) const; +}; + +} // namespace wmtk::submesh diff --git a/src/wmtk/submesh/utils/submesh_from_multimesh.cpp b/src/wmtk/submesh/utils/submesh_from_multimesh.cpp new file mode 100644 index 0000000000..9748a87dd0 --- /dev/null +++ b/src/wmtk/submesh/utils/submesh_from_multimesh.cpp @@ -0,0 +1,40 @@ +#include "submesh_from_multimesh.hpp" + +#include +#include + +namespace wmtk::submesh::utils { + +std::shared_ptr submesh_from_multimesh(const std::shared_ptr& mesh) +{ + std::map, std::shared_ptr> smm; + return submesh_from_multimesh(mesh, smm); +} + +std::shared_ptr submesh_from_multimesh( + const std::shared_ptr& mesh, + std::map, std::shared_ptr>& submesh_map) +{ + // log_assert(mesh->is_multi_mesh_root(), "submesh_from_multimesh must be called on root mesh"); + + std::shared_ptr emb_ptr = std::make_shared(mesh); + Embedding& emb = *emb_ptr; + + for (const auto& child_mesh_ptr : mesh->get_child_meshes()) { + logger().trace("Add submesh from multimesh"); + auto sub_ptr = emb.add_submesh(); + + submesh_map[child_mesh_ptr] = sub_ptr; + + const Mesh& cm = *child_mesh_ptr; + for (const Tuple& child_tuple : cm.get_all(cm.top_simplex_type())) { + const simplex::Simplex child_simplex(cm.top_simplex_type(), child_tuple); + const simplex::Simplex parent_simplex = cm.map_to_parent(child_simplex); + sub_ptr->add_simplex(parent_simplex); + } + } + + return emb_ptr; +} + +} // namespace wmtk::submesh::utils \ No newline at end of file diff --git a/src/wmtk/submesh/utils/submesh_from_multimesh.hpp b/src/wmtk/submesh/utils/submesh_from_multimesh.hpp new file mode 100644 index 0000000000..2bb1d86c1c --- /dev/null +++ b/src/wmtk/submesh/utils/submesh_from_multimesh.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include +#include + +namespace wmtk::submesh::utils { + +/** + * @brief Make the mesh an Embedding and construct SubMeshes from all children. + * + * The method does not delete any child meshes. + * + * @param mesh The root mesh that becomes the embedding mesh. + * @return Shared pointer of the Embedding. + */ +std::shared_ptr submesh_from_multimesh(const std::shared_ptr& mesh); + +/** + * @brief Make the mesh an Embedding and construct SubMeshes from all children. + * + * The method does not delete any child meshes. + * + * @param mesh The root mesh that becomes the embedding mesh. + * @param submesh_map A map from child meshes to the new SubMeshes. + * @return Shared pointer of the Embedding. + */ +std::shared_ptr submesh_from_multimesh( + const std::shared_ptr& mesh, + std::map, std::shared_ptr>& submesh_map); + +} // namespace wmtk::submesh::utils diff --git a/src/wmtk/submesh/utils/write.cpp b/src/wmtk/submesh/utils/write.cpp new file mode 100644 index 0000000000..82850ca41b --- /dev/null +++ b/src/wmtk/submesh/utils/write.cpp @@ -0,0 +1,51 @@ +#include "write.hpp" + +#include +#include +#include +#include + +void wmtk::submesh::utils::write( + const std::filesystem::path& filename, + const std::string& vertices_name, + const SubMesh& sub, + bool write_points, + bool write_edges, + bool write_faces, + bool write_tetrahedra) +{ + const Mesh& m = sub.mesh(); + + ParaviewWriter writer( + filename, + vertices_name, + m, + write_points, + write_edges, + write_faces, + write_tetrahedra, + [&sub](const simplex::IdSimplex& s) { return sub.contains(s); }); + m.serialize(writer); +} + +void wmtk::submesh::utils::write( + const std::filesystem::path& filename, + const std::string& vertices_name, + const Embedding& emb, + bool write_points, + bool write_edges, + bool write_faces, + bool write_tetrahedra) +{ + const Mesh& m = emb.mesh(); + + ParaviewWriter writer( + filename, + vertices_name, + m, + write_points, + write_edges, + write_faces, + write_tetrahedra); + m.serialize(writer); +} diff --git a/src/wmtk/submesh/utils/write.hpp b/src/wmtk/submesh/utils/write.hpp new file mode 100644 index 0000000000..82e6e97f6e --- /dev/null +++ b/src/wmtk/submesh/utils/write.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include + +namespace wmtk::submesh { +class SubMesh; +class Embedding; +} // namespace wmtk::submesh + + +namespace wmtk::submesh::utils { + +void write( + const std::filesystem::path& filename, + const std::string& vertices_name, + const SubMesh& sub, + bool write_points = true, + bool write_edges = true, + bool write_faces = true, + bool write_tetrahedra = true); + +void write( + const std::filesystem::path& filename, + const std::string& vertices_name, + const Embedding& emb, + bool write_points = true, + bool write_edges = true, + bool write_faces = true, + bool write_tetrahedra = true); + +} // namespace wmtk::submesh::utils diff --git a/src/wmtk/utils/CMakeLists.txt b/src/wmtk/utils/CMakeLists.txt index a72bb1a053..825e2e1372 100644 --- a/src/wmtk/utils/CMakeLists.txt +++ b/src/wmtk/utils/CMakeLists.txt @@ -72,6 +72,11 @@ set(SRC_FILES DynamicArray.hpp DynamicArray.hxx + + bbox_from_mesh.hpp + bbox_from_mesh.cpp + bvh_from_mesh.hpp + bvh_from_mesh.cpp ) target_sources(wildmeshing_toolkit PRIVATE ${SRC_FILES}) diff --git a/src/wmtk/utils/Logger.cpp b/src/wmtk/utils/Logger.cpp index aa7fcc8214..c99a8b8a92 100644 --- a/src/wmtk/utils/Logger.cpp +++ b/src/wmtk/utils/Logger.cpp @@ -6,9 +6,10 @@ #include namespace wmtk { -bool has_user_overloaded_logger_level() { +bool has_user_overloaded_logger_level() +{ const char* val = std::getenv("WMTK_LOGGER_LEVEL"); - if(val == nullptr) { + if (val == nullptr) { return false; } std::string env_val = val; @@ -104,4 +105,11 @@ void log_and_throw_error(const std::string& msg) throw std::runtime_error(msg); } +void log_assert(const bool check, const std::string& msg) +{ + if (!check) { + log_and_throw_error(msg); + } +} + } // namespace wmtk diff --git a/src/wmtk/utils/Logger.hpp b/src/wmtk/utils/Logger.hpp index d8aa7498c4..5e8c5b5b7e 100644 --- a/src/wmtk/utils/Logger.hpp +++ b/src/wmtk/utils/Logger.hpp @@ -45,4 +45,20 @@ template { log_and_throw_error(fmt::format(fmt::runtime(msg), args...)); } + +/** + * @brief Calls log_and_throw_error if check is false. + */ +void log_assert(const bool check, const std::string& msg); + +/** + * @brief Calls log_and_throw_error if check is false. + */ +template +[[noreturn]] void log_assert(const bool check, const std::string& msg, const Args&... args) +{ + if (!check) { + log_and_throw_error(msg, args...); + } +} } // namespace wmtk diff --git a/src/wmtk/utils/bbox_from_mesh.cpp b/src/wmtk/utils/bbox_from_mesh.cpp new file mode 100644 index 0000000000..0e17c12140 --- /dev/null +++ b/src/wmtk/utils/bbox_from_mesh.cpp @@ -0,0 +1,59 @@ +#include "bbox_from_mesh.hpp" + +#include + +#include +#include + +namespace wmtk::utils { + + +Eigen::MatrixXd bbox_from_mesh(const attribute::MeshAttributeHandle& position_handle) +{ + const Mesh& m = position_handle.mesh(); + + const int64_t dim = position_handle.dimension(); + + Eigen::MatrixXd bbox; + bbox.resize(2, dim); + for (int64_t d = 0; d < dim; ++d) { + bbox(0, d) = std::numeric_limits::max(); + bbox(1, d) = std::numeric_limits::lowest(); + } + + std::visit( + [&bbox, &dim, &m](const auto& v) -> void { + using T = typename std::decay_t::Type; + const auto p_acc = m.create_const_accessor(v); + for (const Tuple& t : m.get_all(PrimitiveType::Vertex)) { + const auto p = p_acc.const_vector_attribute(t); + for (int64_t d = 0; d < dim; ++d) { + if constexpr (std::is_same_v) { + bbox(0, d) = std::min(bbox(0, d), p[d]); + bbox(1, d) = std::max(bbox(1, d), p[d]); + } else if constexpr (std::is_same_v) { + bbox(0, d) = std::min(bbox(0, d), p[d].to_double()); + bbox(1, d) = std::max(bbox(1, d), p[d].to_double()); + } else { + log_and_throw_error("Cannot compute bbox. Unknown position handle type."); + } + } + } + }, + position_handle.handle()); + + // std::cout << "BBOX:\n" << bbox << std::endl; + + return bbox; +} + +double bbox_diagonal_from_mesh(const attribute::MeshAttributeHandle& position_handle) +{ + const Eigen::MatrixXd bbox = bbox_from_mesh(position_handle); + + const double d = (bbox.row(0) - bbox.row(1)).norm(); + + return d; +} + +} // namespace wmtk::utils \ No newline at end of file diff --git a/src/wmtk/utils/bbox_from_mesh.hpp b/src/wmtk/utils/bbox_from_mesh.hpp new file mode 100644 index 0000000000..3d850cba73 --- /dev/null +++ b/src/wmtk/utils/bbox_from_mesh.hpp @@ -0,0 +1,23 @@ +#pragma once + +#include + +namespace wmtk::utils { + +/** + * @brief Compute the bounding box of a mesh. + * + * @param position_handle The position handle, which also contains a pointer to the mesh. + * @return A matrix where the first row contains all minimal coordinates and the second row contains + * all maximal coordinates. + */ +Eigen::MatrixXd bbox_from_mesh(const attribute::MeshAttributeHandle& position_handle); + +/** + * @brief Compute the diagonal of the bounding box of a mesh. + * @param position_handle The position handle, which also contains a pointer to the mesh. + * @return The diagonal of the bounding box. + */ +double bbox_diagonal_from_mesh(const attribute::MeshAttributeHandle& position_handle); + +} // namespace wmtk::utils \ No newline at end of file diff --git a/src/wmtk/utils/bvh_from_mesh.cpp b/src/wmtk/utils/bvh_from_mesh.cpp new file mode 100644 index 0000000000..ccb1866aa1 --- /dev/null +++ b/src/wmtk/utils/bvh_from_mesh.cpp @@ -0,0 +1,217 @@ +#include "bvh_from_mesh.hpp" + +#include +#include +#include + +namespace wmtk::components::internal::utils { + + +std::shared_ptr bvh_from_mesh(const attribute::MeshAttributeHandle& position_handle) +{ + const Mesh& mesh = position_handle.mesh(); + + constexpr PrimitiveType PV = PrimitiveType::Vertex; + constexpr PrimitiveType PE = PrimitiveType::Edge; + + const attribute::Accessor accessor = + mesh.create_const_accessor(position_handle.as()); + + Eigen::MatrixXd V; + Eigen::MatrixXi F; + + if (mesh.top_simplex_type() == PrimitiveType::Triangle) { + int64_t count = 0; + assert(accessor.dimension() <= 3); + + const std::vector& face_tuples = mesh.get_all(PrimitiveType::Triangle); + + V.resize(3 * face_tuples.size(), accessor.dimension()); + V.setZero(); + F.resize(face_tuples.size(), 3); + + for (const Tuple& f : face_tuples) { + const int64_t fid = f.global_cid(); + + auto p0 = accessor.const_vector_attribute(f); + auto p1 = accessor.const_vector_attribute(mesh.switch_tuple(f, PV)); + auto p2 = accessor.const_vector_attribute(mesh.switch_tuples(f, {PE, PV})); + + F.row(fid) = Eigen::Vector3i(count, count + 1, count + 2); + V.row(3 * fid) = p0; + V.row(3 * fid + 1) = p1; + V.row(3 * fid + 2) = p2; + + count += 3; + } + + + } else if (mesh.top_simplex_type() == PrimitiveType::Edge) { + int64_t count = 0; + + const std::vector& edge_tuples = mesh.get_all(PrimitiveType::Edge); + + V.resize(2 * edge_tuples.size(), accessor.dimension()); + F.resize(edge_tuples.size(), 2); + + for (const Tuple& e : edge_tuples) { + const int64_t eid = e.global_cid(); + + auto p0 = accessor.const_vector_attribute(e); + auto p1 = accessor.const_vector_attribute(mesh.switch_tuple(e, PV)); + + F.row(eid) = Eigen::Vector2i(count, count + 1); + V.row(2 * eid) = p0; + V.row(2 * eid + 1) = p1; + + count += 2; + } + + } else { + log_and_throw_error("bvh_from_mesh works only for tri/edges meshes"); + } + + + std::shared_ptr bvh = std::make_shared(); + bvh->init(V, F, 1e-10); + return bvh; +} + +std::shared_ptr bvh_from_mesh( + const attribute::MeshAttributeHandle& position_handle, + const PrimitiveType pt) + +{ + const Mesh& mesh = position_handle.mesh(); + + constexpr PrimitiveType PV = PrimitiveType::Vertex; + constexpr PrimitiveType PE = PrimitiveType::Edge; + + const attribute::Accessor accessor = + mesh.create_const_accessor(position_handle.as()); + + Eigen::MatrixXd V; + Eigen::MatrixXi F; + + if (pt == PrimitiveType::Triangle) { + int64_t count = 0; + assert(accessor.dimension() <= 3); + + const std::vector& face_tuples = mesh.get_all(PrimitiveType::Triangle); + + V.resize(3 * face_tuples.size(), accessor.dimension()); + V.setZero(); + F.resize(face_tuples.size(), 3); + + for (int64_t fid = 0; fid < face_tuples.size(); ++fid) { + const Tuple& f = face_tuples[fid]; + + auto p0 = accessor.const_vector_attribute(f); + auto p1 = accessor.const_vector_attribute(mesh.switch_tuple(f, PV)); + auto p2 = accessor.const_vector_attribute(mesh.switch_tuples(f, {PE, PV})); + + F.row(fid) = Eigen::Vector3i(count, count + 1, count + 2); + V.row(3 * fid) = p0; + V.row(3 * fid + 1) = p1; + V.row(3 * fid + 2) = p2; + + count += 3; + } + + + } else if (pt == PrimitiveType::Edge) { + int64_t count = 0; + + const std::vector& edge_tuples = mesh.get_all(PrimitiveType::Edge); + + V.resize(2 * edge_tuples.size(), accessor.dimension()); + F.resize(edge_tuples.size(), 2); + + for (int64_t eid = 0; eid < edge_tuples.size(); ++eid) { + const Tuple& e = edge_tuples[eid]; + + auto p0 = accessor.const_vector_attribute(e); + auto p1 = accessor.const_vector_attribute(mesh.switch_tuple(e, PV)); + + F.row(eid) = Eigen::Vector2i(count, count + 1); + V.row(2 * eid) = p0; + V.row(2 * eid + 1) = p1; + + count += 2; + } + + } else { + log_and_throw_error("bvh_from_mesh works only for tri/edges meshes"); + } + + + std::shared_ptr bvh = std::make_shared(); + bvh->init(V, F, 1e-10); + return bvh; +} + +std::shared_ptr bvh_from_mesh( + attribute::MeshAttributeHandle& position_handle, + attribute::MeshAttributeHandle& label_handle, + const int64_t label_value) +{ + const Mesh& mesh = position_handle.mesh(); + + std::shared_ptr bvh; + + const auto p_acc = mesh.create_const_accessor(position_handle); + const auto tag_acc = mesh.create_const_accessor(label_handle); + + std::vector pts; + + for (const Tuple& t : mesh.get_all(PrimitiveType::Triangle)) { + const simplex::Simplex tri(mesh, PrimitiveType::Triangle, t); + + if (tag_acc.primitive_type() == PrimitiveType::Tetrahedron) { + const auto tets = simplex::top_dimension_cofaces(mesh, tri); + if (tets.size() != 2) { + continue; + } + + // find tets that are on the boundary of the tagged region + if ((tag_acc.const_scalar_attribute(tets.simplex_vector()[0]) == label_value) == + (tag_acc.const_scalar_attribute(tets.simplex_vector()[1]) == label_value)) { + continue; + } + } else { + // find tagged triangles + if (tag_acc.const_scalar_attribute(tri) != label_value) { + continue; + } + } + + const auto vertices = simplex::faces_single_dimension(mesh, tri, PrimitiveType::Vertex); + pts.emplace_back(p_acc.const_vector_attribute(vertices.simplex_vector()[0])); + pts.emplace_back(p_acc.const_vector_attribute(vertices.simplex_vector()[1])); + pts.emplace_back(p_acc.const_vector_attribute(vertices.simplex_vector()[2])); + } + + if (pts.empty()) { + log_and_throw_error("Seems like the input mesh does not contain any tagged elements."); + } + + Eigen::MatrixXd V; + Eigen::MatrixXi F; + V.resize(pts.size(), 3); + F.resize(pts.size() / 3, 3); + for (size_t i = 0; i < F.rows(); ++i) { + const size_t v0 = 3 * i + 0; + const size_t v1 = 3 * i + 1; + const size_t v2 = 3 * i + 2; + V.row(v0) = pts[v0]; + V.row(v1) = pts[v1]; + V.row(v2) = pts[v2]; + F.row(i) = Eigen::Vector3i(v0, v1, v2); + } + bvh = std::make_shared(); + bvh->init(V, F, 1e-10); + + return bvh; +} + +} // namespace wmtk::components::internal::utils \ No newline at end of file diff --git a/src/wmtk/utils/bvh_from_mesh.hpp b/src/wmtk/utils/bvh_from_mesh.hpp new file mode 100644 index 0000000000..cf0e946b00 --- /dev/null +++ b/src/wmtk/utils/bvh_from_mesh.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include + +namespace wmtk::components::internal::utils { + +std::shared_ptr bvh_from_mesh( + const attribute::MeshAttributeHandle& position_handle); + +std::shared_ptr bvh_from_mesh( + const attribute::MeshAttributeHandle& position_handle, + const PrimitiveType pt); + +/** + * @brief Generate a bvh from tagged simplices. + * + * This method assumes that the bvh should contain triangles. If the labels are on tetrahedra, the + * boundary of those is used. + */ +std::shared_ptr bvh_from_mesh( + attribute::MeshAttributeHandle& position_handle, + attribute::MeshAttributeHandle& label_handle, + const int64_t label_value); + +} // namespace wmtk::components::internal::utils \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 30cb0cac6d..b9902bd273 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -103,6 +103,8 @@ add_subdirectory_with_source_group(multimesh) add_subdirectory_with_source_group(operations) add_subdirectory_with_source_group(attributes) add_subdirectory_with_source_group(autogen) +add_subdirectory_with_source_group(submesh) +add_subdirectory_with_source_group(invariants) diff --git a/tests/invariants/CMakeLists.txt b/tests/invariants/CMakeLists.txt new file mode 100644 index 0000000000..7f0ee256ba --- /dev/null +++ b/tests/invariants/CMakeLists.txt @@ -0,0 +1,5 @@ +# Sources +set(TEST_SOURCES + test_envelope_invariant.cpp +) +target_sources(wmtk_tests PRIVATE ${TEST_SOURCES}) diff --git a/tests/invariants/test_envelope_invariant.cpp b/tests/invariants/test_envelope_invariant.cpp new file mode 100644 index 0000000000..f56c420fc8 --- /dev/null +++ b/tests/invariants/test_envelope_invariant.cpp @@ -0,0 +1,371 @@ +#include + +#include +#include +#include +#include +#include +#include +#include "tools/DEBUG_EdgeMesh.hpp" +#include "tools/DEBUG_PointMesh.hpp" +#include "tools/DEBUG_TriMesh.hpp" +#include "tools/EdgeMesh_examples.hpp" +#include "tools/TriMesh_examples.hpp" + +#include + +using namespace wmtk; +using namespace tests; +using namespace invariants; + +constexpr PrimitiveType PV = PrimitiveType::Vertex; +constexpr PrimitiveType PE = PrimitiveType::Edge; +constexpr PrimitiveType PF = PrimitiveType::Triangle; + +void positions_as_rational(Mesh& mesh) +{ + auto pt_double_attribute = mesh.get_attribute_handle("vertices", PrimitiveType::Vertex); + wmtk::utils::cast_attribute(pt_double_attribute, mesh, "vertices"); + mesh.delete_attribute(pt_double_attribute); +} + +std::shared_ptr construct_edge_45(DEBUG_TriMesh& m) +{ + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto acc = m.create_const_accessor(env_pos_handle); + + std::shared_ptr em_ptr = std::make_shared(single_line()); + Eigen::MatrixXd V; + V.resize(2, 3); + V.row(0) = acc.const_vector_attribute(m.vertex_tuple_from_id(4)); + V.row(1) = acc.const_vector_attribute(m.vertex_tuple_from_id(5)); + mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, *em_ptr); + + return em_ptr; +} + +std::shared_ptr construct_point_4(DEBUG_TriMesh& m) +{ + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto acc = m.create_const_accessor(env_pos_handle); + + std::shared_ptr pm_ptr = std::make_shared(1); + Eigen::MatrixXd V; + V.resize(1, 3); + V.row(0) = acc.const_vector_attribute(m.vertex_tuple_from_id(4)); + mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, *pm_ptr); + + return pm_ptr; +} + +TEST_CASE("envelope_invariant_envelope", "[invariants][envelope]") +{ + logger().set_level(spdlog::level::off); + + std::shared_ptr mesh_in = + std::make_shared(edge_region_with_position()); + DEBUG_TriMesh& m = *mesh_in; + + + SECTION("double_triangle") + { + using T = double; + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, env_pos_handle); + + const Tuple v4 = m.vertex_tuple_from_id(4); + const simplex::Simplex s4(PV, v4); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s4); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + acc.vector_attribute(s4)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("rational_triangle") + { + using T = Rational; + positions_as_rational(m); + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, env_pos_handle); + + const Tuple v4 = m.vertex_tuple_from_id(4); + const simplex::Simplex s4(PV, v4); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s4); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + acc.vector_attribute(s4)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("double_edge") + { + std::shared_ptr em_ptr = construct_edge_45(m); + DEBUG_EdgeMesh& em = *em_ptr; + + using T = double; + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto query_pos_handle = em.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0 = em.tuple_from_edge_id(0); + const simplex::Simplex s0(PV, v0); + + + const auto tops = simplex::top_dimension_cofaces_tuples(em, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = em.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("rational_edge") + { + std::shared_ptr em_ptr = construct_edge_45(m); + DEBUG_EdgeMesh& em = *em_ptr; + + using T = Rational; + positions_as_rational(em); + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto query_pos_handle = em.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0 = em.tuple_from_edge_id(0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(em, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = em.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("double_point") + { + std::shared_ptr pm_ptr = construct_point_4(m); + DEBUG_PointMesh& pm = *pm_ptr; + + using T = double; + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto query_pos_handle = pm.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0(-1, -1, -1, 0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(pm, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = pm.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("rational_point") + { + std::shared_ptr pm_ptr = construct_point_4(m); + DEBUG_PointMesh& pm = *pm_ptr; + + using T = Rational; + positions_as_rational(pm); + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + auto query_pos_handle = pm.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0(-1, -1, -1, 0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(pm, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = pm.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } +} + +TEST_CASE("envelope_invariant_bvh", "[invariants][envelope]") +{ + logger().set_level(spdlog::level::off); + + std::shared_ptr em_ptr = std::make_shared(single_line()); + DEBUG_EdgeMesh& em = *em_ptr; + { + Eigen::MatrixXd V; + V.resize(2, 2); + V.row(0) << 0, 0; + V.row(1) << 1, 0; + mesh_utils::set_matrix_attribute(V, "vertices", PV, em); + } + + SECTION("double_edge") + { + using T = double; + + auto env_pos_handle = em.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, env_pos_handle); + + const Tuple v0 = em.tuple_from_edge_id(0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(em, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = em.create_accessor(env_pos_handle); + acc.vector_attribute(s0)[1] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("rational_edge") + { + using T = Rational; + positions_as_rational(em); + + auto env_pos_handle = em.get_attribute_handle("vertices", PV); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, env_pos_handle); + + const Tuple v0 = em.tuple_from_edge_id(0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(em, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = em.create_accessor(env_pos_handle); + acc.vector_attribute(s0)[1] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("double_point") + { + std::shared_ptr pm_ptr = std::make_shared(1); + DEBUG_PointMesh& pm = *pm_ptr; + + using T = double; + + auto env_pos_handle = em.get_attribute_handle("vertices", PV); + auto query_pos_handle = pm.register_attribute("vertices", PV, 2); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0(-1, -1, -1, 0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(pm, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = pm.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[1] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("rational_point") + { + std::shared_ptr pm_ptr = std::make_shared(1); + DEBUG_PointMesh& pm = *pm_ptr; + + using T = Rational; + + auto env_pos_handle = em.get_attribute_handle("vertices", PV); + auto query_pos_handle = pm.register_attribute("vertices", PV, 2); + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, query_pos_handle); + + const Tuple v0(-1, -1, -1, 0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(pm, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = pm.create_accessor(query_pos_handle); + acc.vector_attribute(s0)[1] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } +} + +TEST_CASE("envelope_invariant_submesh_edge", "[invariants][envelope]") +{ + std::shared_ptr mesh_in = + std::make_shared(edge_region_with_position()); + DEBUG_TriMesh& m = *mesh_in; + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + + submesh::Embedding emb(mesh_in); + + auto sub1_ptr = emb.add_submesh(); + submesh::SubMesh& sub1 = *sub1_ptr; + sub1.add_simplex(m.edge_tuple_from_vids(4, 5), PE); + + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, sub1); + + // SECTION("4") + { + const Tuple v4 = m.vertex_tuple_from_id(4); + const simplex::Simplex s4(PV, v4); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s4); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + acc.vector_attribute(s4)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + // SECTION("0") + { + const Tuple v0 = m.vertex_tuple_from_id(0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + + acc.vector_attribute(s0)[2] = 1; + CHECK(env_inv.after({}, tops)); + } +} + +TEST_CASE("envelope_invariant_submesh_triangle", "[invariants][envelope]") +{ + std::shared_ptr mesh_in = + std::make_shared(edge_region_with_position()); + DEBUG_TriMesh& m = *mesh_in; + + auto env_pos_handle = m.get_attribute_handle("vertices", PV); + + submesh::Embedding emb(mesh_in); + + auto sub_ptr = emb.add_submesh(); + submesh::SubMesh& sub = *sub_ptr; + sub.add_simplex(m.face_tuple_from_vids(4, 5, 1), PF); + + EnvelopeInvariant env_inv(env_pos_handle, 1e-2, sub); + + SECTION("4") + { + const Tuple v4 = m.vertex_tuple_from_id(4); + const simplex::Simplex s4(PV, v4); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s4); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + acc.vector_attribute(s4)[2] = 1; + CHECK_FALSE(env_inv.after({}, tops)); + } + SECTION("0") + { + const Tuple v0 = m.vertex_tuple_from_id(0); + const simplex::Simplex s0(PV, v0); + + const auto tops = simplex::top_dimension_cofaces_tuples(m, s0); + CHECK(env_inv.after({}, tops)); + + auto acc = m.create_accessor(env_pos_handle); + + acc.vector_attribute(s0)[2] = 1; + CHECK(env_inv.after({}, tops)); + } +} \ No newline at end of file diff --git a/tests/operations/CMakeLists.txt b/tests/operations/CMakeLists.txt index 8ee8ac1c34..51ef97a39f 100644 --- a/tests/operations/CMakeLists.txt +++ b/tests/operations/CMakeLists.txt @@ -11,5 +11,6 @@ set(TEST_SOURCES test_attribute_transfer.cpp vertex_optimization.cpp edge_operation_data.cpp + test_projection_operation.cpp ) target_sources(wmtk_tests PRIVATE ${TEST_SOURCES}) diff --git a/tests/operations/test_projection_operation.cpp b/tests/operations/test_projection_operation.cpp new file mode 100644 index 0000000000..9332281407 --- /dev/null +++ b/tests/operations/test_projection_operation.cpp @@ -0,0 +1,229 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace wmtk; +using namespace operations; +using namespace composite; + +void write(const Mesh& mesh, const std::string& name) +{ + if (mesh.top_simplex_type() == PrimitiveType::Triangle) { + wmtk::io::ParaviewWriter writer(name, "vertices", mesh, true, true, true, false); + mesh.serialize(writer); + } else if (mesh.top_simplex_type() == PrimitiveType::Edge) { + wmtk::io::ParaviewWriter writer(name, "vertices", mesh, true, true, false, false); + mesh.serialize(writer); + } +} + +attribute::MeshAttributeHandle register_and_transfer_positions( + attribute::MeshAttributeHandle parent_pos_handle, + Mesh& child) +{ + auto child_pos_handle = child.register_attribute( + "vertices", + PrimitiveType::Vertex, + parent_pos_handle.dimension()); + + attribute_update::make_cast_attribute_transfer_strategy( + /*source=*/parent_pos_handle, + /*target=*/child_pos_handle) + ->run_on_all(); + + return child_pos_handle; +} + +TEST_CASE("test_projection_operation", "[operations][.]") +{ + logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + opt_logger().set_level(spdlog::level::off); + + std::shared_ptr mesh_ptr = + std::make_shared(tests::edge_region_with_position_2D()); + tests::DEBUG_TriMesh& m = *mesh_ptr; + + // add input (edge 4-5) + std::shared_ptr child_input_ptr; + { + child_input_ptr = std::make_shared(tests::single_line()); + + const Tuple child_tuple = child_input_ptr->edge_tuple_from_vids(0, 1); + const Tuple parent_tuple = m.edge_tuple_from_vids(4, 5); + std::vector> map_tuples; + map_tuples.emplace_back(std::array{child_tuple, parent_tuple}); + m.register_child_mesh(child_input_ptr, map_tuples); + } + // add boundary + std::shared_ptr child_bnd_ptr; + { + child_bnd_ptr = std::make_shared(); + tests::DEBUG_EdgeMesh& c = *child_bnd_ptr; + + RowVectors2l edges; + edges.resize(8, 2); + edges.row(0) << 0, 1; + edges.row(1) << 1, 2; + edges.row(2) << 2, 3; + edges.row(3) << 3, 4; + edges.row(4) << 4, 5; + edges.row(5) << 5, 6; + edges.row(6) << 6, 7; + edges.row(7) << 7, 0; + c.initialize(edges); + + std::vector> map_tuples; + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(0, 1), m.edge_tuple_from_vids(0, 1)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(1, 2), m.edge_tuple_from_vids(1, 2)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(2, 3), m.edge_tuple_from_vids(2, 6)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(3, 4), m.edge_tuple_from_vids(6, 9)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(4, 5), m.edge_tuple_from_vids(9, 8)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(5, 6), m.edge_tuple_from_vids(8, 7)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(6, 7), m.edge_tuple_from_vids(7, 3)}); + map_tuples.emplace_back( + std::array{c.edge_tuple_from_vids(7, 0), m.edge_tuple_from_vids(3, 0)}); + m.register_child_mesh(child_bnd_ptr, map_tuples); + } + + attribute::MeshAttributeHandle pos_handle = + m.get_attribute_handle("vertices", PrimitiveType::Vertex); + + // move vertices around + { + auto acc = m.create_accessor(pos_handle); + const Tuple t4 = m.vertex_tuple_from_id(4); + acc.vector_attribute(t4)[0] = 0.75; + acc.vector_attribute(t4)[1] = 0.9; + } + + SECTION("multimesh") + { + auto write_all = [&]() { + static int64_t write_counter = 0; + write(m, fmt::format("test_proj_mm_parent_{}", write_counter)); + write(*child_input_ptr, fmt::format("test_proj_mm_input_{}", write_counter)); + write(*child_bnd_ptr, fmt::format("test_proj_mm_bnd_{}", write_counter)); + ++write_counter; + }; + + auto child_input_pos_handle = register_and_transfer_positions(pos_handle, *child_input_ptr); + auto child_bnd_pos_handle = register_and_transfer_positions(pos_handle, *child_bnd_ptr); + + std::vector> update_child_position; + update_child_position.emplace_back(attribute_update::make_cast_attribute_transfer_strategy( + /*source=*/pos_handle, + /*target=*/child_input_pos_handle)); + update_child_position.emplace_back(attribute_update::make_cast_attribute_transfer_strategy( + /*source=*/pos_handle, + /*target=*/child_bnd_pos_handle)); + std::vector> update_parent_position; + update_parent_position.emplace_back(attribute_update::make_cast_attribute_transfer_strategy( + /*source=*/child_input_pos_handle, + /*target=*/pos_handle)); + update_parent_position.emplace_back(attribute_update::make_cast_attribute_transfer_strategy( + /*source=*/child_bnd_pos_handle, + /*target=*/pos_handle)); + + std::vector mesh_constraint_pairs; + mesh_constraint_pairs.emplace_back(child_input_pos_handle, child_input_pos_handle); + mesh_constraint_pairs.emplace_back(child_bnd_pos_handle, child_bnd_pos_handle); + + auto smoothing = std::make_shared(m, pos_handle); + for (auto& s : update_child_position) { + smoothing->add_transfer_strategy(s); + } + auto proj_smoothing = std::make_shared(m, mesh_constraint_pairs); + for (auto& s : update_parent_position) { + proj_smoothing->add_transfer_strategy(s); + } + for (auto& s : update_child_position) { + proj_smoothing->add_transfer_strategy(s); + } + + const simplex::Simplex v1(PrimitiveType::Vertex, m.vertex_tuple_from_id(1)); + + write_all(); + (*smoothing)(v1); + write_all(); + CHECK_NOTHROW((*proj_smoothing)(v1)); + + write_all(); + } + SECTION("submesh") + { + std::map, std::shared_ptr> submesh_map; + auto emb_ptr = submesh::utils::submesh_from_multimesh(mesh_ptr, submesh_map); + submesh::Embedding& emb = *emb_ptr; + auto sub_input_ptr = submesh_map[child_input_ptr]; + auto sub_bnd_ptr = submesh_map[child_bnd_ptr]; + submesh::SubMesh& sub_input = *sub_input_ptr; + submesh::SubMesh& sub_bnd = *sub_bnd_ptr; + + auto write_all = [&]() { + static int64_t write_counter = 0; + submesh::utils::write( + fmt::format("test_proj_sub_parent_{}", write_counter), + "vertices", + emb, + false, + true, + true, + false); + submesh::utils::write( + fmt::format("test_proj_sub_input_{}", write_counter), + "vertices", + sub_input, + false, + true, + false, + false); + submesh::utils::write( + fmt::format("test_proj_sub_bnd_{}", write_counter), + "vertices", + sub_bnd, + false, + true, + false, + false); + ++write_counter; + }; + + // deregister all child meshes + for (const auto [child, sub] : submesh_map) { + m.deregister_child_mesh(child); + } + emb.update_tag_attribute_handles(); + + write_all(); + + auto smoothing = std::make_shared(m, pos_handle); + auto proj_smoothing = std::make_shared(emb, pos_handle); + + const simplex::Simplex v1(PrimitiveType::Vertex, m.vertex_tuple_from_id(1)); + + (*smoothing)(v1); + write_all(); + CHECK_NOTHROW((*proj_smoothing)(v1)); + + write_all(); + } +} \ No newline at end of file diff --git a/tests/submesh/CMakeLists.txt b/tests/submesh/CMakeLists.txt new file mode 100644 index 0000000000..f6c22241fe --- /dev/null +++ b/tests/submesh/CMakeLists.txt @@ -0,0 +1,5 @@ +# Sources +set(TEST_SOURCES + test_submesh.cpp +) +target_sources(wmtk_tests PRIVATE ${TEST_SOURCES}) diff --git a/tests/submesh/test_submesh.cpp b/tests/submesh/test_submesh.cpp new file mode 100644 index 0000000000..50bd65b021 --- /dev/null +++ b/tests/submesh/test_submesh.cpp @@ -0,0 +1,429 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "tools/DEBUG_TriMesh.hpp" +#include "tools/TriMesh_examples.hpp" + + +using namespace wmtk; +using namespace submesh; + +constexpr PrimitiveType PV = PrimitiveType::Vertex; +constexpr PrimitiveType PE = PrimitiveType::Edge; +constexpr PrimitiveType PF = PrimitiveType::Triangle; +constexpr PrimitiveType PT = PrimitiveType::Tetrahedron; + +TEST_CASE("submesh_init", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + const Tuple edge45 = m.edge_tuple_from_vids(4, 5); + + Embedding emb(mesh_in); + std::shared_ptr sub_ptr = emb.add_submesh(); + SubMesh& sub = *sub_ptr; + + CHECK_THROWS(sub.top_simplex_type()); + + sub.add_simplex(edge45, PE); + CHECK(sub.top_simplex_type() == PrimitiveType::Edge); + CHECK(sub.top_cell_dimension() == 1); + + CHECK(sub.contains(edge45, PE)); + CHECK(sub.contains(edge45, PV)); + CHECK(sub.contains(sub.switch_tuple(edge45, PV), PV)); + CHECK(sub.is_boundary(edge45, PE)); + CHECK(sub.is_boundary(edge45, PV)); + CHECK(sub.is_boundary(sub.switch_tuple(edge45, PV), PV)); + CHECK(sub.top_simplex_type(edge45) == PE); + + // Test switch_tuple_vector + const Tuple edge34 = m.edge_tuple_from_vids(3, 4); + sub.add_simplex(edge34, PE); + CHECK(sub.contains(edge34, PE)); + CHECK(sub.contains(edge34, PV)); + CHECK(sub.contains(sub.switch_tuple(edge34, PV), PV)); + CHECK(sub.is_boundary(edge34, PV)); + CHECK(!sub.is_boundary(sub.switch_tuple(edge34, PV), PV)); + CHECK(!sub.is_boundary(edge45, PV)); + CHECK(sub.top_simplex_type(edge34) == PE); + + // switch from edge 45 to 43 + { + CHECK_THROWS(sub.switch_tuple(edge34, PE)); + CHECK_NOTHROW(sub.switch_tuple(edge45, PE)); + Tuple sw_edge45; + REQUIRE_NOTHROW(sw_edge45 = sub.switch_tuple(edge45, PE)); // global switch + CHECK(m.get_id_simplex(sw_edge45, PE) == m.get_id_simplex(edge34, PE)); + CHECK(m.get_id_simplex(sw_edge45, PV) == m.get_id_simplex(edge45, PV)); + } + + const Tuple edge04 = m.edge_tuple_from_vids(0, 4); + sub.add_simplex(edge04, PE); + CHECK(sub.is_boundary(edge04, PV)); + CHECK(!sub.is_boundary(sub.switch_tuple(edge04, PV), PV)); + + const Tuple edge54 = sub.switch_tuple(edge45, PV); + CHECK_THROWS(sub.switch_tuple(edge54, PE)); + CHECK(m.get_id_simplex(sub.switch_tuple(edge54, PV), PV) == m.get_id_simplex(edge45, PV)); + + + const Tuple face596 = m.face_tuple_from_vids(5, 9, 6); + sub.add_simplex(face596, PF); + CHECK(sub.top_simplex_type() == PF); + CHECK(sub.id(face596, PF) == 9); + CHECK(sub.is_boundary(face596, PF)); + CHECK(sub.is_boundary(m.vertex_tuple_from_id(5), PV)); + // after adding a face, all edges and vertices that are not incident to a face, are boundary + CHECK(sub.is_boundary(edge45, PE)); + CHECK(sub.is_boundary(edge45, PV)); + const Tuple edge59 = m.edge_tuple_with_vs_and_t(5, 9, 9); + { + // local edge switch + CHECK_NOTHROW(sub.switch_tuple(edge59, PE)); + // global face switch + CHECK_THROWS(sub.switch_tuple(face596, PF)); + } + + + { + ParaviewWriter writer("submesh_init", "vertices", m, false, true, true, false); + CHECK_NOTHROW(m.serialize(writer)); + } +} + +TEST_CASE("submesh_init_from_tag", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + + Embedding emb(mesh_in); + std::shared_ptr sub_ptr = emb.add_submesh(); + SubMesh& sub = *sub_ptr; + + // register tag attribute + { + auto tag_handle = m.register_attribute("tag", PF, 1); + auto acc = m.create_accessor(tag_handle); + acc.scalar_attribute(m.tuple_from_face_id(0)) = 1; + acc.scalar_attribute(m.tuple_from_face_id(2)) = 1; + + sub.add_from_tag_attribute(tag_handle.as(), 1); + } + + CHECK(sub.top_simplex_type() == PF); + CHECK(sub.get_all(PF).size() == 2); + CHECK(sub.get_all(PE).size() == 6); + CHECK(sub.get_all(PV).size() == 5); +} + +TEST_CASE("submesh_top_dimension_cofaces", "[mesh][submesh][.]") +{ + REQUIRE(false); // test not implemented + + // logger().set_level(spdlog::level::off); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + const Tuple edge45 = m.edge_tuple_from_vids(4, 5); + + Embedding emb(mesh_in); + std::shared_ptr sub_ptr = emb.add_submesh(); + SubMesh& sub = *sub_ptr; + + CHECK_THROWS(sub.top_simplex_type()); + + sub.add_simplex(edge45, PE); +} + +TEST_CASE("submesh_init_multi", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + + Embedding emb(mesh_in); + CHECK_FALSE(emb.has_child_mesh()); + + std::shared_ptr sub1_ptr = emb.add_submesh(); + SubMesh& sub1 = *sub1_ptr; + + std::shared_ptr sub2_ptr = emb.add_submesh(); + SubMesh& sub2 = *sub2_ptr; + + CHECK(emb.has_child_mesh()); + + sub1.add_simplex(m.face_tuple_from_vids(0, 3, 4), PF); + sub1.add_simplex(m.face_tuple_from_vids(1, 4, 5), PF); + + sub2.add_simplex(m.face_tuple_from_vids(1, 4, 5), PF); + + { + using submesh::utils::write; + CHECK_NOTHROW(write("submesh_init_multi", "vertices", emb, false, false, true, false)); + CHECK_NOTHROW( + write("submesh_init_multi_sub1", "vertices", sub1, false, false, true, false)); + CHECK_NOTHROW( + write("submesh_init_multi_sub2", "vertices", sub2, false, false, true, false)); + } +} + +TEST_CASE("submesh_split", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + + const auto pos = m.get_attribute_handle("vertices", PrimitiveType::Vertex); + + Embedding emb(mesh_in); + CHECK_FALSE(emb.has_child_mesh()); + + std::shared_ptr sub1_ptr = emb.add_submesh(); + sub1_ptr->add_simplex(m.face_tuple_from_vids(1, 4, 5), PF); + + std::shared_ptr sub2_ptr = emb.add_submesh(); + sub2_ptr->add_simplex(m.face_tuple_from_vids(4, 5, 8), PF); + + CHECK(emb.has_child_mesh()); + + + // operations::EdgeSplit split(m); + // emb.set_split_strategies(split); + operations::EdgeSplit split(emb); + split.set_new_attribute_strategy(pos); + + const simplex::Simplex edge45(PE, m.edge_tuple_from_vids(4, 5)); + split(edge45); + + CHECK(sub1_ptr->contains(m.vertex_tuple_from_id(4), PV)); + CHECK(sub1_ptr->contains(m.edge_tuple_from_vids(4, 10), PE)); + CHECK(sub1_ptr->contains(m.edge_tuple_from_vids(5, 10), PE)); + CHECK(sub1_ptr->contains(m.face_tuple_from_vids(4, 10, 1), PF)); + CHECK(sub1_ptr->contains(m.face_tuple_from_vids(5, 10, 1), PF)); + CHECK(sub2_ptr->contains(m.vertex_tuple_from_id(4), PV)); + CHECK(sub2_ptr->contains(m.edge_tuple_from_vids(4, 10), PE)); + CHECK(sub2_ptr->contains(m.edge_tuple_from_vids(5, 10), PE)); + CHECK(sub2_ptr->contains(m.face_tuple_from_vids(4, 10, 8), PF)); + CHECK(sub2_ptr->contains(m.face_tuple_from_vids(5, 10, 8), PF)); + + { + using submesh::utils::write; + CHECK_NOTHROW(write("submesh_split", "vertices", emb, true, true, true, false)); + } +} + +TEST_CASE("submesh_collapse", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + + const auto pos = m.get_attribute_handle("vertices", PrimitiveType::Vertex); + + Embedding emb(mesh_in); + CHECK_FALSE(emb.has_child_mesh()); + + std::shared_ptr sub1_ptr = emb.add_submesh(); + sub1_ptr->add_simplex(m.face_tuple_from_vids(1, 4, 5), PF); + + std::shared_ptr sub2_ptr = emb.add_submesh(); + sub2_ptr->add_simplex(m.face_tuple_from_vids(4, 5, 8), PF); + + CHECK(emb.has_child_mesh()); + + + // operations::EdgeCollapse collapse(m); + // emb.set_collapse_strategies(collapse); + operations::EdgeCollapse collapse(emb); + collapse.set_new_attribute_strategy(pos); + + const simplex::Simplex edge45(PE, m.edge_tuple_from_vids(4, 5)); + collapse(edge45); + + CHECK(sub1_ptr->contains(m.vertex_tuple_from_id(5), PV)); + CHECK(sub1_ptr->contains(m.edge_tuple_from_vids(5, 1), PE)); + CHECK(sub2_ptr->contains(m.vertex_tuple_from_id(5), PV)); + CHECK(sub2_ptr->contains(m.edge_tuple_from_vids(5, 8), PE)); + + { + using submesh::utils::write; + CHECK_NOTHROW(write("submesh_collapse", "vertices", emb, true, true, true, false)); + } +} + +TEST_CASE("submesh_collapse_towards_submesh", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + + const auto pos = m.get_attribute_handle("vertices", PrimitiveType::Vertex); + const auto pos_acc = m.create_const_accessor(pos); + + Embedding emb(mesh_in); + CHECK_FALSE(emb.has_child_mesh()); + + std::shared_ptr sub1_ptr = emb.add_submesh(); + sub1_ptr->add_simplex(m.edge_tuple_from_vids(1, 4), PE); + sub1_ptr->add_simplex(m.edge_tuple_from_vids(4, 8), PE); + + // std::shared_ptr sub2_ptr = emb.add_submesh(); + // sub2_ptr->add_simplex(m.face_tuple_from_vids(4, 5, 8), PF); + Eigen::VectorXd p4 = pos_acc.const_vector_attribute(m.vertex_tuple_from_id(4)); + Eigen::VectorXd p5 = pos_acc.const_vector_attribute(m.vertex_tuple_from_id(5)); + CHECK(p4 != p5); + + CHECK(emb.has_child_mesh()); + + + operations::EdgeCollapse collapse(emb); + // emb.set_collapse_strategies(collapse); + auto clps_strat = std::make_shared>(pos); + clps_strat->set_simplex_predicate(emb.substructure_predicate()); + clps_strat->set_strategy(operations::CollapseBasicStrategy::Default); + + collapse.set_new_attribute_strategy(pos, clps_strat); + + const simplex::Simplex edge45(PE, m.edge_tuple_from_vids(4, 5)); + collapse(edge45); + + CHECK(sub1_ptr->contains(m.vertex_tuple_from_id(5), PV)); + CHECK(sub1_ptr->contains(m.edge_tuple_from_vids(5, 1), PE)); + + p5 = pos_acc.const_vector_attribute(m.vertex_tuple_from_id(5)); + CHECK(p4 == p5); + + { + using submesh::utils::write; + CHECK_NOTHROW( + write("submesh_collapse_towards_submesh", "vertices", emb, true, true, true, false)); + } +} + +TEST_CASE("submesh_from_multimesh", "[mesh][submesh]") +{ + // logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + // add child + { + std::shared_ptr child_ptr = + std::make_shared(tests::single_triangle()); + + const Tuple child_tuple = child_ptr->face_tuple_from_vids(0, 1, 2); + const Tuple parent_tuple = mesh_in->face_tuple_from_vids(8, 9, 5); + std::vector> map_tuples; + map_tuples.emplace_back(std::array{child_tuple, parent_tuple}); + mesh_in->register_child_mesh(child_ptr, map_tuples); + } + // add child + { + std::shared_ptr child_ptr = + std::make_shared(tests::single_triangle()); + + const Tuple child_tuple = child_ptr->face_tuple_from_vids(0, 1, 2); + const Tuple parent_tuple = mesh_in->face_tuple_from_vids(0, 3, 4); + std::vector> map_tuples; + map_tuples.emplace_back(std::array{child_tuple, parent_tuple}); + mesh_in->register_child_mesh(child_ptr, map_tuples); + } + + tests::DEBUG_TriMesh& m = *mesh_in; + + const auto pos = m.get_attribute_handle("vertices", PrimitiveType::Vertex); + const auto pos_acc = m.create_const_accessor(pos); + + auto emb_ptr = submesh::utils::submesh_from_multimesh(mesh_in); + + Embedding& emb = *emb_ptr; + CHECK(emb.has_child_mesh()); + CHECK(emb.get_child_meshes().size() == 2); + { + using submesh::utils::write; + CHECK_NOTHROW(write("submesh_from_multimesh", "vertices", emb, true, true, true, false)); + } +} + +TEST_CASE("submesh_type_casts", "[mesh][submesh]") +{ + logger().set_level(spdlog::level::off); + // logger().set_level(spdlog::level::trace); + + // basic test for implementing + std::shared_ptr mesh_in = + std::make_shared(tests::edge_region_with_position()); + + tests::DEBUG_TriMesh& m = *mesh_in; + CHECK_FALSE(m.has_embedding()); + const Tuple edge45 = m.edge_tuple_from_vids(4, 5); + + std::shared_ptr emb_ptr = std::make_shared(mesh_in); + CHECK(m.has_embedding()); + REQUIRE(emb_ptr->mesh_type() == MeshType::Embedding); + Embedding& emb = static_cast(*emb_ptr); + + std::shared_ptr sub_ptr = emb.add_submesh(); + CHECK(emb.get_child_meshes().size() == 1); + REQUIRE(sub_ptr->mesh_type() == MeshType::SubMesh); + SubMesh& sub = static_cast(*sub_ptr); + + sub.add_simplex(edge45, PE); + CHECK(sub.contains(edge45, PE)); + + // test linking of embedding in mesh + { + Mesh& emb_mesh = emb.mesh(); + REQUIRE(emb_mesh.has_embedding()); + Embedding& mesh_emb = emb_mesh.get_embedding(); + CHECK(&mesh_emb == &emb); + } +} \ No newline at end of file diff --git a/tests/tools/DEBUG_TriMesh.cpp b/tests/tools/DEBUG_TriMesh.cpp index f9dc176a80..9842625e98 100644 --- a/tests/tools/DEBUG_TriMesh.cpp +++ b/tests/tools/DEBUG_TriMesh.cpp @@ -77,6 +77,9 @@ auto DEBUG_TriMesh::edge_tuple_from_vids(const int64_t v1, const int64_t v2) con const attribute::Accessor fv = create_const_accessor(m_fv_handle); auto& fv_base = create_base_accessor(m_fv_handle); for (int64_t fid = 0; fid < capacity(PrimitiveType::Triangle); ++fid) { + if (is_removed(fid)) { + continue; + } Tuple face = face_tuple_from_id(fid); auto fv0 = fv.const_vector_attribute(face); int64_t local_vid1 = -1, local_vid2 = -1; diff --git a/tests/tools/TriMesh_examples.cpp b/tests/tools/TriMesh_examples.cpp index 2b59810063..774a348d2f 100644 --- a/tests/tools/TriMesh_examples.cpp +++ b/tests/tools/TriMesh_examples.cpp @@ -287,6 +287,26 @@ TriMesh edge_region_with_position() return m; } +TriMesh edge_region_with_position_2D() +{ + TriMesh m = edge_region(); + + Eigen::MatrixXd V; + V.resize(10, 2); + V.row(0) << 0.5, 1; + V.row(1) << 1.5, 1; + V.row(2) << 2.5, 1; + V.row(3) << 0, 0; + V.row(4) << 1, 0; + V.row(5) << 2, 0; + V.row(6) << 3, 0; + V.row(7) << 0.5, -1; + V.row(8) << 1.5, -1; + V.row(9) << 2.5, -1; + mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, m); + return m; +} + TriMesh embedded_diamond() { // 0---1 diff --git a/tests/tools/TriMesh_examples.hpp b/tests/tools/TriMesh_examples.hpp index e996639e74..72477d032f 100644 --- a/tests/tools/TriMesh_examples.hpp +++ b/tests/tools/TriMesh_examples.hpp @@ -176,6 +176,8 @@ TriMesh nine_triangles_with_a_hole(); TriMesh ten_triangles_with_position(int dimension); TriMesh edge_region_with_position(); + +TriMesh edge_region_with_position_2D(); // 0---1 // / \ / \ . // 2---3---4