From 9fe4e481da9cf2445da41e05acf51a900068ab22 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 30 Sep 2024 03:38:56 -0400 Subject: [PATCH 01/12] add cdt optimization component --- applications/cdt_sec/CMakeLists.txt | 1 + applications/cdt_sec/cdt_sec_main.cpp | 4 +- components/CMakeLists.txt | 1 + .../cdt_optimization/CMakeLists.txt | 13 + .../cdt_optimization/cdt_optimization.cpp | 304 ++++++++++++++++++ .../cdt_optimization/cdt_optimization.hpp | 18 ++ src/wmtk/invariants/CMakeLists.txt | 3 + .../invariants/Swap2dEdgeLengthInvariant.cpp | 39 +++ .../invariants/Swap2dEdgeLengthInvariant.hpp | 20 ++ 9 files changed, 402 insertions(+), 1 deletion(-) create mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt create mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp create mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp create mode 100644 src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp create mode 100644 src/wmtk/invariants/Swap2dEdgeLengthInvariant.hpp diff --git a/applications/cdt_sec/CMakeLists.txt b/applications/cdt_sec/CMakeLists.txt index bce93b336c..55754bab37 100644 --- a/applications/cdt_sec/CMakeLists.txt +++ b/applications/cdt_sec/CMakeLists.txt @@ -11,6 +11,7 @@ target_link_libraries(cdt_sec_app PRIVATE wmtk::input wmtk::multimesh wmtk::shortestedge_collapse +wmtk::cdt_optimization wmtk::CDT wmtk::output) diff --git a/applications/cdt_sec/cdt_sec_main.cpp b/applications/cdt_sec/cdt_sec_main.cpp index d692195c96..c7c83a5dd3 100644 --- a/applications/cdt_sec/cdt_sec_main.cpp +++ b/applications/cdt_sec/cdt_sec_main.cpp @@ -12,6 +12,7 @@ #include +#include #include #include #include @@ -102,10 +103,11 @@ int main(int argc, char* argv[]) parent_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); - auto mesh_after_sec = wmtk::components::shortestedge_collapse( + auto mesh_after_sec = wmtk::components::cdt_optimization( static_cast(*child_mesh), child_mesh_position_handle, parent_mesh_position_handle, + 5, false, j["length_rel"], false, diff --git a/components/CMakeLists.txt b/components/CMakeLists.txt index a259ebd2bc..7317710528 100644 --- a/components/CMakeLists.txt +++ b/components/CMakeLists.txt @@ -29,6 +29,7 @@ add_subdirectory("tetwild_simplification/wmtk/components/tetwild_simplification" add_subdirectory("multimesh/wmtk/components/multimesh") add_subdirectory("shortestedge_collapse/wmtk/components/shortestedge_collapse") add_subdirectory("CDT/wmtk/components/CDT") +add_subdirectory("cdt_optimization/wmtk/components/cdt_optimization") # add_component("export_cache") # add_component("import_cache") diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt b/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt new file mode 100644 index 0000000000..68e0434734 --- /dev/null +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt @@ -0,0 +1,13 @@ +set(COMPONENT_NAME cdt_optimization) +add_component(${COMPONENT_NAME}) +if(NOT ${WMTK_ENABLE_COMPONENT_${COMPONENT_NAME}}) + return() +endif() + +set(SRC_FILES + cdt_optimization.hpp + cdt_optimization.cpp + ) + + +target_sources(wmtk_${COMPONENT_NAME} PRIVATE ${SRC_FILES}) diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp new file mode 100644 index 0000000000..90ee5e4e5b --- /dev/null +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp @@ -0,0 +1,304 @@ +#include "cdt_optimization.hpp" + +#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 { +std::shared_ptr cdt_optimization( + TriMesh& mesh, + const attribute::MeshAttributeHandle& position_handle, + std::optional& inversion_position_handle, + const int64_t passes, + bool update_other_position, + const double length_rel, + bool lock_boundary, + double envelope_size, + const std::vector& pass_through) +{ + if (mesh.top_simplex_type() != PrimitiveType::Triangle) { + log_and_throw_error( + "shortest edge collapse works only for triangle meshes: {}", + primitive_type_name(mesh.top_simplex_type())); + } + + // TriMesh& mesh = static_cast(position_handle.mesh()); + + auto pass_through_attributes = pass_through; + + bool check_inversion = inversion_position_handle.has_value(); + + ///////////////////////////////////////////// + + auto visited_edge_flag = + mesh.register_attribute("visited_edge", PrimitiveType::Edge, 1, false, char(1)); + + auto update_flag_func = [](Eigen::Ref 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, + position_handle, + update_flag_func); + + ////////////////////////////////// + // 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 = [](Eigen::Ref P) -> Eigen::VectorXd { + assert(P.cols() == 2); + assert(P.rows() == 2 || P.rows() == 3); + return Eigen::VectorXd::Constant(1, (P.col(0) - P.col(1)).norm()); + }; + auto edge_length_update = + std::make_shared>( + edge_length_attribute, + position_handle, + compute_edge_length); + edge_length_update->run_on_all(); + + ////////////////////////////////// + // computing bbox diagonal + Eigen::VectorXd bmin(mesh.top_cell_dimension()); + bmin.setConstant(std::numeric_limits::max()); + Eigen::VectorXd bmax(mesh.top_cell_dimension()); + bmax.setConstant(std::numeric_limits::lowest()); + + auto pt_accessor = mesh.create_const_accessor(position_handle); + + const auto vertices = mesh.get_all(PrimitiveType::Vertex); + for (const auto& v : vertices) { + const auto p = pt_accessor.vector_attribute(v); + for (int64_t d = 0; d < bmax.size(); ++d) { + bmin[d] = std::min(bmin[d], p[d]); + bmax[d] = std::max(bmax[d], p[d]); + } + } + + const double bbdiag = (bmax - bmin).norm(); + + const double length_abs = bbdiag * length_rel; + + wmtk::logger().info( + "bbox max {}, bbox min {}, diag {}, target edge length {}", + bmax, + bmin, + bbdiag, + length_abs); + + auto short_edges_first_priority = [&](const simplex::Simplex& s) { + assert(s.primitive_type() == PrimitiveType::Edge); + return edge_length_accessor.const_scalar_attribute(s.tuple()); + }; + pass_through_attributes.push_back(edge_length_attribute); + auto todo = std::make_shared( + mesh, + edge_length_attribute.as(), + 4. / 5. * length_abs); // MTAO: why is this 4/5? + + //////////////////////////invariants + + auto invariant_link_condition = std::make_shared(mesh); + + auto invariant_interior_edge = std::make_shared(mesh); + auto invariant_interior_vertex = std::make_shared(mesh); + + auto set_all_invariants = [&](auto&& m) { + invariant_interior_edge->add( + std::make_shared(m, PrimitiveType::Edge)); + invariant_interior_vertex->add( + std::make_shared(m, PrimitiveType::Vertex)); + }; + multimesh::MultiMeshVisitor visitor(set_all_invariants); + visitor.execute_from_root(mesh); + + auto invariant_mm_map = std::make_shared(mesh); + + ////////////// positions + std::vector positions; + positions.push_back(position_handle); + if (check_inversion) { + positions.push_back(inversion_position_handle.value()); + } + + ////////////////////////////////////////// + // collapse + ////////////////////////////////////////// + auto collapse = std::make_shared(mesh); + collapse->add_invariant(todo); + collapse->add_invariant(invariant_link_condition); + collapse->add_invariant(invariant_mm_map); + + + if (envelope_size > 0) { + collapse->add_invariant(std::make_shared( + position_handle, + bbdiag * envelope_size, + position_handle)); + } + + + if (check_inversion) { + collapse->add_invariant(std::make_shared>( + inversion_position_handle.value().mesh(), + inversion_position_handle.value().as())); + } + collapse->set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::CollapseBasicStrategy::None); + collapse->add_transfer_strategy(tag_update); + + collapse->set_priority(short_edges_first_priority); + collapse->add_transfer_strategy(edge_length_update); + + if (lock_boundary) { + collapse->add_invariant(invariant_interior_edge); + // set collapse towards boundary + for (auto& p : positions) { + auto tmp = std::make_shared>(p); + tmp->set_strategy(wmtk::operations::CollapseBasicStrategy::CopyOther); + tmp->set_simplex_predicate(wmtk::operations::BasicSimplexPredicate::IsInterior); + collapse->set_new_attribute_strategy(p, tmp); + } + } else { + for (auto& p : positions) { + collapse->set_new_attribute_strategy( + p, + wmtk::operations::CollapseBasicStrategy::CopyOther); + } + } + + + for (const auto& attr : pass_through_attributes) { + collapse->set_new_attribute_strategy(attr); + } + + ////////////////////////////////////////// + // swap + ////////////////////////////////////////// + + auto swap = std::make_shared(mesh); + + swap->add_invariant( + std::make_shared(mesh, position_handle.as())); + swap->add_invariant(invariant_mm_map); + + swap->collapse().add_invariant(invariant_link_condition); + + if (envelope_size > 0) { + swap->add_invariant(std::make_shared( + position_handle, + bbdiag * envelope_size, + position_handle)); + } + + + if (check_inversion) { + swap->add_invariant(std::make_shared>( + inversion_position_handle.value().mesh(), + inversion_position_handle.value().as())); + } + + swap->split().set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + swap->collapse().set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::CollapseBasicStrategy::None); + + swap->add_transfer_strategy(tag_update); + + swap->split().add_transfer_strategy(edge_length_update); + swap->collapse().add_transfer_strategy(edge_length_update); + + for (auto& p : positions) { + swap->split().set_new_attribute_strategy(p); + swap->collapse().set_new_attribute_strategy( + p, + wmtk::operations::CollapseBasicStrategy::CopyOther); + } + + swap->set_priority(short_edges_first_priority); + + swap->add_transfer_strategy(edge_length_update); + + for (const auto& attr : pass_through_attributes) { + swap->split().set_new_attribute_strategy(attr); + swap->collapse().set_new_attribute_strategy(attr); + } + + ////////////////////////////////////////// + // scheduler + ////////////////////////////////////////// + + Scheduler scheduler; + + for (int64_t i = 0; i < passes; ++i) { + auto collapse_stats = + scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); + + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + "collapse", + collapse_stats.number_of_performed_operations(), + collapse_stats.number_of_successful_operations(), + collapse_stats.number_of_failed_operations(), + collapse_stats.collecting_time, + collapse_stats.sorting_time, + collapse_stats.executing_time); + + auto swap_stats = scheduler.run_operation_on_all(*swap, visited_edge_flag.as()); + + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + "swap", + swap_stats.number_of_performed_operations(), + swap_stats.number_of_successful_operations(), + swap_stats.number_of_failed_operations(), + swap_stats.collecting_time, + swap_stats.sorting_time, + swap_stats.executing_time); + + multimesh::consolidate(mesh); + } + + + // // debug code + // for (const auto& e : mesh.get_all(PrimitiveType::Edge)) { + // if (mesh.is_boundary(PrimitiveType::Edge, e)) { + // logger().error("after sec mesh has nonmanifold edges"); + // } + // } + + + return mesh.shared_from_this(); +} +} // namespace wmtk::components diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp new file mode 100644 index 0000000000..f34727db50 --- /dev/null +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp @@ -0,0 +1,18 @@ +#pragma once +#include +#include + +namespace wmtk::components { + +std::shared_ptr cdt_optimization( + TriMesh& mesh, + const attribute::MeshAttributeHandle& position_handle, + std::optional& inversion_position_handle, + const int64_t passes, + bool update_other_position, + const double length_rel, + bool lock_boundary, + double envelope_size, + const std::vector& pass_through); + +} // namespace wmtk::components diff --git a/src/wmtk/invariants/CMakeLists.txt b/src/wmtk/invariants/CMakeLists.txt index 42e91e948b..c59dad4cac 100644 --- a/src/wmtk/invariants/CMakeLists.txt +++ b/src/wmtk/invariants/CMakeLists.txt @@ -92,6 +92,9 @@ set(SRC_FILES EnergyFilterInvariant.hpp EnergyFilterInvariant.cpp + + Swap2dEdgeLengthInvariant.hpp + Swap2dEdgeLengthInvariant.cpp ) target_sources(wildmeshing_toolkit PRIVATE ${SRC_FILES}) diff --git a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp new file mode 100644 index 0000000000..51ece7635e --- /dev/null +++ b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp @@ -0,0 +1,39 @@ +#include "Swap2dEdgeLengthInvariant.hpp" +#include +#include +#include + +namespace wmtk { +Swap2dEdgeLengthInvariant::Swap2dEdgeLengthInvariant( + const Mesh& m, + const attribute::TypedAttributeHandle& coordinate) + : Invariant(m, true, false, false) + , m_coordinate_handle(coordinate) +{} + +bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const +{ + assert(m.top_cell_dimension() == 2); + constexpr static PrimitiveType PV = PrimitiveType::Vertex; + constexpr static PrimitiveType PE = PrimitiveType::Edge; + constexpr static PrimitiveType PF = PrimitiveType::Triangle; + + auto accessor = mesh().create_const_accessor(m_coordinate_handle); + + // get the coords of the vertices + // input face end points + const Tuple v0 = t.tuple(); + const Tuple v1 = mesh().switch_tuple(v0, PV); + // other 2 vertices + const Tuple v2 = mesh().switch_tuples(v0, {PE, PV}); + const Tuple v3 = mesh().switch_tuples(v0, {PF, PE, PV}); + + const double length_old = + (accessor.const_vector_attribute(v0) - accessor.const_vector_attribute(v1)).norm(); + const double length_new = + (accessor.const_vector_attribute(v2) - accessor.const_vector_attribute(v3)).norm(); + + return length_new > length_old; +} + +} // namespace wmtk diff --git a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.hpp b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.hpp new file mode 100644 index 0000000000..7c66908ca8 --- /dev/null +++ b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include +#include "Invariant.hpp" + +namespace wmtk { +class Swap2dEdgeLengthInvariant : public Invariant +{ +public: + Swap2dEdgeLengthInvariant( + const Mesh& m, + const attribute::TypedAttributeHandle& coordinate); + using Invariant::Invariant; + + bool before(const simplex::Simplex& t) const override; + +private: + const attribute::TypedAttributeHandle m_coordinate_handle; +}; +} // namespace wmtk From 322775246a900a28d5affd97c9c83216531acaf9 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 30 Sep 2024 07:25:02 -0400 Subject: [PATCH 02/12] add debug code --- .../cdt_optimization/cdt_optimization.cpp | 53 +++++++++++++++---- .../invariants/Swap2dEdgeLengthInvariant.cpp | 30 ++++++++--- 2 files changed, 68 insertions(+), 15 deletions(-) diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp index 90ee5e4e5b..af7d1b0d0d 100644 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp @@ -3,10 +3,13 @@ #include #include #include +#include +#include #include #include #include #include +#include #include #include #include @@ -22,6 +25,8 @@ #include #include #include +#include +#include "predicates.h" namespace wmtk::components { @@ -202,10 +207,16 @@ std::shared_ptr cdt_optimization( // swap ////////////////////////////////////////// + std::shared_ptr amips = + std::make_shared(mesh, position_handle); + auto function_invariant = + std::make_shared(mesh.top_simplex_type(), amips); + auto swap = std::make_shared(mesh); swap->add_invariant( std::make_shared(mesh, position_handle.as())); + swap->add_invariant(invariant_mm_map); swap->collapse().add_invariant(invariant_link_condition); @@ -219,11 +230,13 @@ std::shared_ptr cdt_optimization( if (check_inversion) { - swap->add_invariant(std::make_shared>( + swap->collapse().add_invariant(std::make_shared>( inversion_position_handle.value().mesh(), inversion_position_handle.value().as())); } + // swap->add_invariant(function_invariant); + swap->split().set_new_attribute_strategy( visited_edge_flag, wmtk::operations::SplitBasicStrategy::None, @@ -234,8 +247,7 @@ std::shared_ptr cdt_optimization( swap->add_transfer_strategy(tag_update); - swap->split().add_transfer_strategy(edge_length_update); - swap->collapse().add_transfer_strategy(edge_length_update); + swap->add_transfer_strategy(edge_length_update); for (auto& p : positions) { swap->split().set_new_attribute_strategy(p); @@ -291,12 +303,35 @@ std::shared_ptr cdt_optimization( } - // // debug code - // for (const auto& e : mesh.get_all(PrimitiveType::Edge)) { - // if (mesh.is_boundary(PrimitiveType::Edge, e)) { - // logger().error("after sec mesh has nonmanifold edges"); - // } - // } + // debug code + for (const auto& e : mesh.get_all(PrimitiveType::Edge)) { + if (mesh.is_boundary(PrimitiveType::Edge, e)) { + logger().error("after sec mesh has nonmanifold edges"); + } + } + + auto inversion_pos_accessor = + inversion_position_handle.value().mesh().create_const_accessor( + inversion_position_handle.value()); + + // debug code + for (const auto& t : + inversion_position_handle.value().mesh().get_all(PrimitiveType::Tetrahedron)) { + const auto vertices = inversion_position_handle.value().mesh().orient_vertices(t); + std::vector pos; + for (int i = 0; i < vertices.size(); ++i) { + pos.push_back(inversion_pos_accessor.const_vector_attribute(vertices[i])); + } + if (orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data()) <= 0) { + wmtk::logger().error( + "Flipped or degenerated tet! orient3d = {}", + orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data())); + } else if (orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data()) <= 1e-5) { + wmtk::logger().error( + "Nearly degenerated tet! orient3d = {}", + orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data())); + } + } return mesh.shared_from_this(); diff --git a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp index 51ece7635e..f3f0374b56 100644 --- a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp +++ b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp @@ -13,7 +13,7 @@ Swap2dEdgeLengthInvariant::Swap2dEdgeLengthInvariant( bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const { - assert(m.top_cell_dimension() == 2); + assert(mesh().top_cell_dimension() == 2); constexpr static PrimitiveType PV = PrimitiveType::Vertex; constexpr static PrimitiveType PE = PrimitiveType::Edge; constexpr static PrimitiveType PF = PrimitiveType::Triangle; @@ -28,12 +28,30 @@ bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const const Tuple v2 = mesh().switch_tuples(v0, {PE, PV}); const Tuple v3 = mesh().switch_tuples(v0, {PF, PE, PV}); - const double length_old = - (accessor.const_vector_attribute(v0) - accessor.const_vector_attribute(v1)).norm(); - const double length_new = - (accessor.const_vector_attribute(v2) - accessor.const_vector_attribute(v3)).norm(); + // use length + // const double length_old = + // (accessor.const_vector_attribute(v0) - accessor.const_vector_attribute(v1)).norm(); + // const double length_new = + // (accessor.const_vector_attribute(v2) - accessor.const_vector_attribute(v3)).norm(); - return length_new > length_old; + // return length_new > length_old; + + // use height + + Eigen::Vector3d p0 = accessor.const_vector_attribute(v0); + Eigen::Vector3d p1 = accessor.const_vector_attribute(v1); + Eigen::Vector3d p2 = accessor.const_vector_attribute(v2); + Eigen::Vector3d p3 = accessor.const_vector_attribute(v3); + + const double h0 = ((p0 - p2).cross(p0 - p3)).norm() / (p2 - p3).norm(); + const double h1 = ((p1 - p2).cross(p1 - p3)).norm() / (p2 - p3).norm(); + const double h2 = ((p2 - p0).cross(p2 - p1)).norm() / (p0 - p1).norm(); + const double h3 = ((p3 - p0).cross(p3 - p1)).norm() / (p0 - p1).norm(); + + const double min_old = std::min(h0, h1); + const double min_new = std::min(h2, h3); + + return min_old < min_new; } } // namespace wmtk From 51079cbe4aee80104106b20234478bacbfe434d9 Mon Sep 17 00:00:00 2001 From: JcDai Date: Tue, 8 Oct 2024 12:11:12 +0800 Subject: [PATCH 03/12] add split draft --- .../cdt_optimization/cdt_optimization.cpp | 49 ++++++++++++++++++- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp index af7d1b0d0d..086139040b 100644 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -120,12 +121,21 @@ std::shared_ptr cdt_optimization( assert(s.primitive_type() == PrimitiveType::Edge); return edge_length_accessor.const_scalar_attribute(s.tuple()); }; + auto long_edges_first_priority = [&](const simplex::Simplex& s) { + assert(s.primitive_type() == PrimitiveType::Edge); + return -edge_length_accessor.const_scalar_attribute(s.tuple()); + }; pass_through_attributes.push_back(edge_length_attribute); - auto todo = std::make_shared( + auto todo_smaller = std::make_shared( mesh, edge_length_attribute.as(), 4. / 5. * length_abs); // MTAO: why is this 4/5? + auto todo_larger= = std::make_shared( + mesh, + edge_length_attribute.as(), + 4.0 / 3.0); + //////////////////////////invariants auto invariant_link_condition = std::make_shared(mesh); @@ -151,11 +161,46 @@ std::shared_ptr cdt_optimization( positions.push_back(inversion_position_handle.value()); } + ////////////////////////////////////////// + // split + ////////////////////////////////////////// + + auto split = std::make_shared(mesh); + split->add_invariant(todo_larger); + if (check_inversion) { + split->add_invariant(std::make_shared>( + inversion_position_handle.value().mesh(), + inversion_position_handle.value().as())); + } + + split->set_priority(long_edges_first_priority); + + split->set_new_attribute_strategy( + visited_edge_flag, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + + split->add_transfer_strategy(tag_update); + split->add_transfer_strategy(edge_length_update); + + for (auto& p : positions) { + split->set_new_attribute_strategy(p); + } + + for (const auto& attr : pass_through_attributes) { + split->set_new_attribute_strategy( + attr, + wmtk::operations::SplitBasicStrategy::None, + wmtk::operations::SplitRibBasicStrategy::None); + } + + + ////////////////////////////////////////// // collapse ////////////////////////////////////////// auto collapse = std::make_shared(mesh); - collapse->add_invariant(todo); + collapse->add_invariant(todo_smaller); collapse->add_invariant(invariant_link_condition); collapse->add_invariant(invariant_mm_map); From d0ad2296dc6ea2a910914066758fb0599c161684 Mon Sep 17 00:00:00 2001 From: JcDai Date: Tue, 8 Oct 2024 00:50:59 -0400 Subject: [PATCH 04/12] add split fix swap --- .../cdt_optimization/cdt_optimization.cpp | 49 +++++++++++++++++-- .../invariants/Swap2dEdgeLengthInvariant.cpp | 8 ++- 2 files changed, 51 insertions(+), 6 deletions(-) diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp index 086139040b..e16d64cd82 100644 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp +++ b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp @@ -131,10 +131,10 @@ std::shared_ptr cdt_optimization( edge_length_attribute.as(), 4. / 5. * length_abs); // MTAO: why is this 4/5? - auto todo_larger= = std::make_shared( + auto todo_larger = std::make_shared( mesh, edge_length_attribute.as(), - 4.0 / 3.0); + 4.0 / 3.0 * length_abs); //////////////////////////invariants @@ -172,7 +172,7 @@ std::shared_ptr cdt_optimization( inversion_position_handle.value().mesh(), inversion_position_handle.value().as())); } - + split->set_priority(long_edges_first_priority); split->set_new_attribute_strategy( @@ -195,7 +195,6 @@ std::shared_ptr cdt_optimization( } - ////////////////////////////////////////// // collapse ////////////////////////////////////////// @@ -316,7 +315,34 @@ std::shared_ptr cdt_optimization( Scheduler scheduler; + auto collapse_stats_pre = + scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); + + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + "collapse", + collapse_stats_pre.number_of_performed_operations(), + collapse_stats_pre.number_of_successful_operations(), + collapse_stats_pre.number_of_failed_operations(), + collapse_stats_pre.collecting_time, + collapse_stats_pre.sorting_time, + collapse_stats_pre.executing_time); + for (int64_t i = 0; i < passes; ++i) { + auto split_stats = scheduler.run_operation_on_all(*split, visited_edge_flag.as()); + + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + "split", + split_stats.number_of_performed_operations(), + split_stats.number_of_successful_operations(), + split_stats.number_of_failed_operations(), + split_stats.collecting_time, + split_stats.sorting_time, + split_stats.executing_time); + auto collapse_stats = scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); @@ -347,6 +373,21 @@ std::shared_ptr cdt_optimization( multimesh::consolidate(mesh); } + auto collapse_stats_after = + scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); + + logger().info( + "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + "collapse", + collapse_stats_after.number_of_performed_operations(), + collapse_stats_after.number_of_successful_operations(), + collapse_stats_after.number_of_failed_operations(), + collapse_stats_after.collecting_time, + collapse_stats_after.sorting_time, + collapse_stats_after.executing_time); + + multimesh::consolidate(mesh); // debug code for (const auto& e : mesh.get_all(PrimitiveType::Edge)) { diff --git a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp index f3f0374b56..c667e09dbc 100644 --- a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp +++ b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp @@ -48,8 +48,12 @@ bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const const double h2 = ((p2 - p0).cross(p2 - p1)).norm() / (p0 - p1).norm(); const double h3 = ((p3 - p0).cross(p3 - p1)).norm() / (p0 - p1).norm(); - const double min_old = std::min(h0, h1); - const double min_new = std::min(h2, h3); + const double min_old = std::min(h2, h3); + const double min_new = std::min(h0, h1); + + if (((p0 - p1).norm() * h2 / 2) < 1e-10 || ((p0 - p1).norm() * h3 / 2) < 1e-10) { + return true; + } return min_old < min_new; } From a0bd3f7265d046ad366e9600e9f1389e98cd4336 Mon Sep 17 00:00:00 2001 From: JcDai Date: Wed, 9 Oct 2024 03:42:25 -0400 Subject: [PATCH 05/12] hack on small tris --- src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp index c667e09dbc..ad74cbb452 100644 --- a/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp +++ b/src/wmtk/invariants/Swap2dEdgeLengthInvariant.cpp @@ -52,7 +52,9 @@ bool Swap2dEdgeLengthInvariant::before(const simplex::Simplex& t) const const double min_new = std::min(h0, h1); if (((p0 - p1).norm() * h2 / 2) < 1e-10 || ((p0 - p1).norm() * h3 / 2) < 1e-10) { - return true; + if (((p2 - p3).norm() * h0 / 2) > 1e-10 && ((p2 - p3).norm() * h1 / 2) > 1e-10) { + return true; + } } return min_old < min_new; From 77c5b6fc300aae15975b984c25c718c64186648d Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 17:07:56 -0500 Subject: [PATCH 06/12] target edge length bug --- .../wmtk/components/wildmeshing/internal/wildmeshing3d.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing3d.cpp b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing3d.cpp index 343c1b4e26..89096a93f6 100644 --- a/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing3d.cpp +++ b/components/wildmeshing/wmtk/components/wildmeshing/internal/wildmeshing3d.cpp @@ -226,7 +226,7 @@ std::vector, std::string>> wildmeshing3d( PrimitiveType::Edge, 1, false, - options.target_edge_length); // defaults to target edge length + target_edge_length); // defaults to target edge length // Target edge length update const double min_edge_length = [&]() -> double { From 388a79e3e76dfb0f1bf1762ce9ca97167a419348 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 17:08:19 -0500 Subject: [PATCH 07/12] cdt opt app --- applications/cdt_sec/CMakeLists.txt | 1 + applications/cdt_sec/cdt_sec_main.cpp | 82 +++++++++++++------ .../CDT/wmtk/components/CDT/internal/CDT.cpp | 26 ++++++ 3 files changed, 86 insertions(+), 23 deletions(-) diff --git a/applications/cdt_sec/CMakeLists.txt b/applications/cdt_sec/CMakeLists.txt index 840cf741e6..bdd8f0ba83 100644 --- a/applications/cdt_sec/CMakeLists.txt +++ b/applications/cdt_sec/CMakeLists.txt @@ -11,6 +11,7 @@ wmtk::multimesh wmtk::cdt_optimization wmtk::shortest_edge_collapse wmtk::CDT +wmtk::wildmeshing wmtk::output) wmtk_register_integration_test(EXEC_NAME cdt_sec_app diff --git a/applications/cdt_sec/cdt_sec_main.cpp b/applications/cdt_sec/cdt_sec_main.cpp index 614d8d4d7c..481aa0ab6d 100644 --- a/applications/cdt_sec/cdt_sec_main.cpp +++ b/applications/cdt_sec/cdt_sec_main.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include "cdt_sec_spec.hpp" @@ -58,7 +59,7 @@ int main(int argc, char* argv[]) auto mesh = wmtk::components::input::input(input_file); wmtk::logger().info("mesh has {} vertices", mesh->get_all(PrimitiveType::Vertex).size()); - auto mesh_after_cdt = wmtk::components::CDT(static_cast(*mesh), true, false); + auto mesh_after_cdt = wmtk::components::CDT(static_cast(*mesh), true, true); attribute::MeshAttributeHandle mesh_after_cdt_position_handle; @@ -96,28 +97,63 @@ int main(int argc, char* argv[]) parent_mesh->get_attribute_handle("is_boundary", PrimitiveType::Triangle); pass_through.push_back(boundary_handle); - auto child_mesh_position_handle = - child_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); - - attribute::MeshAttributeHandle parent_mesh_position_handle = - parent_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); - - { - using namespace components::shortest_edge_collapse; - - ShortestEdgeCollapseOptions options; - options.position_handle = child_mesh_position_handle; - options.other_position_handles.emplace_back(parent_mesh_position_handle); - options.length_rel = j["length_rel"]; - options.envelope_size = j["envelope_size"]; - options.check_inversions = true; - options.pass_through_attributes = pass_through; - - shortest_edge_collapse(*parent_mesh, options); - } - - wmtk::components::output::output(*parent_mesh, output_file, "vertices"); - wmtk::components::output::output(*child_mesh, output_file + "_surface", "vertices"); + // auto child_mesh_position_handle = + // child_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); + + // attribute::MeshAttributeHandle parent_mesh_position_handle = + // parent_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); + + // { + // using namespace components::shortest_edge_collapse; + + // ShortestEdgeCollapseOptions options; + // options.position_handle = child_mesh_position_handle; + // options.other_position_handles.emplace_back(parent_mesh_position_handle); + // options.length_rel = j["length_rel"]; + // options.envelope_size = j["envelope_size"]; + // options.check_inversions = true; + // options.pass_through_attributes = pass_through; + + // shortest_edge_collapse(*parent_mesh, options); + // } + + std::vector enves; + + wmtk::components::EnvelopeOptions e_surface; + e_surface.envelope_name = "surface"; + e_surface.envelope_constrained_mesh = child_mesh; + e_surface.envelope_geometry_mesh = child_mesh; + e_surface.constrained_position_name = "vertices"; + e_surface.geometry_position_name = "vertices"; + e_surface.thickness = j["envelope_size"]; + + enves.push_back(e_surface); + + wmtk::components::WildMeshingOptions wmo; + wmo.input_mesh = parent_mesh; + wmo.input_mesh_position = "vertices"; + wmo.target_edge_length = j["length_rel"]; + wmo.target_max_amips = 50; + wmo.max_passes = 10; + wmo.intermediate_output = false; + wmo.replace_double_coordinate = false; + wmo.scheduler_update_frequency = 0; + wmo.intermediate_output_path = ""; + wmo.intermediate_output_name = j["output"]; + wmo.envelopes = enves; + wmo.pass_through = pass_through; + wmo.skip_split = false; + wmo.skip_collapse = false; + wmo.skip_swap = false; + wmo.skip_smooth = true; + + auto meshes_after_wildmeshing = wildmeshing(wmo); + + wmtk::components::output::output(*meshes_after_wildmeshing[0].first, output_file, "vertices"); + wmtk::components::output::output( + *meshes_after_wildmeshing[1].first, + output_file + "_surface", + "vertices"); const std::string report = j["report"]; if (!report.empty()) { diff --git a/components/CDT/wmtk/components/CDT/internal/CDT.cpp b/components/CDT/wmtk/components/CDT/internal/CDT.cpp index 296ecab1c2..f99484cf3e 100644 --- a/components/CDT/wmtk/components/CDT/internal/CDT.cpp +++ b/components/CDT/wmtk/components/CDT/internal/CDT.cpp @@ -2,9 +2,14 @@ // #include "CDT.hpp" +#include +#include #include #include #include + +#include + #include "cdt_lib.hpp" #include "get_vf.hpp" @@ -79,6 +84,27 @@ std::shared_ptr CDT_internal( } mesh_utils::set_matrix_attribute(V, "vertices", PrimitiveType::Vertex, *tm); + + auto rounding_pt_attribute = + tm->get_attribute_handle_typed("vertices", PrimitiveType::Vertex); + + + auto rounding = std::make_shared(*tm, rounding_pt_attribute); + rounding->add_invariant( + std::make_shared>(*tm, rounding_pt_attribute)); + + Scheduler scheduler; + auto stats = scheduler.run_operation_on_all(*rounding); + + logger().info( + "Executed rounding, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " + "executing: {}", + stats.number_of_performed_operations(), + stats.number_of_successful_operations(), + stats.number_of_failed_operations(), + stats.collecting_time, + stats.sorting_time, + stats.executing_time); } else { MatrixX V_double; V_double.resize(V_final_str.size(), 3); From 232f3ecbcfc65b1dbc87b661a662eb16e030e699 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 17:17:42 -0500 Subject: [PATCH 08/12] rename --- applications/CMakeLists.txt | 2 +- applications/{cdt_sec => cdt_opt}/CMakeLists.txt | 12 ++++++------ .../cdt_sec_main.cpp => cdt_opt/cdt_opt_main.cpp} | 6 +++--- .../cdt_sec_spec.hpp => cdt_opt/cdt_opt_spec.hpp} | 2 +- .../cdt_opt_test_config.json} | 2 +- 5 files changed, 12 insertions(+), 12 deletions(-) rename applications/{cdt_sec => cdt_opt}/CMakeLists.txt (51%) rename applications/{cdt_sec/cdt_sec_main.cpp => cdt_opt/cdt_opt_main.cpp} (97%) rename applications/{cdt_sec/cdt_sec_spec.hpp => cdt_opt/cdt_opt_spec.hpp} (95%) rename applications/{cdt_sec/cdt_sec_test_config.json => cdt_opt/cdt_opt_test_config.json} (87%) diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index a04324ed43..07e8a5df4b 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -30,7 +30,7 @@ add_application(isotropic_remeshing OFF) add_application(tetwild_simplification ON) add_application(triwild ON) add_application(tetwild ON) -add_application(cdt_sec ON) +add_application(cdt_opt ON) add_application(shortest_edge_collapse ON) add_application(insertion ON) add_application(longest_edge_split ON) diff --git a/applications/cdt_sec/CMakeLists.txt b/applications/cdt_opt/CMakeLists.txt similarity index 51% rename from applications/cdt_sec/CMakeLists.txt rename to applications/cdt_opt/CMakeLists.txt index bdd8f0ba83..234c2bed64 100644 --- a/applications/cdt_sec/CMakeLists.txt +++ b/applications/cdt_opt/CMakeLists.txt @@ -1,11 +1,11 @@ -wmtk_add_application(cdt_sec_app - cdt_sec_main.cpp - cdt_sec_spec.hpp +wmtk_add_application(cdt_opt_app + cdt_opt_main.cpp + cdt_opt_spec.hpp ) -target_link_libraries(cdt_sec_app PRIVATE +target_link_libraries(cdt_opt_app PRIVATE wmtk::input wmtk::multimesh wmtk::cdt_optimization @@ -14,7 +14,7 @@ wmtk::CDT wmtk::wildmeshing wmtk::output) -wmtk_register_integration_test(EXEC_NAME cdt_sec_app - CONFIG_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cdt_sec_test_config.json +wmtk_register_integration_test(EXEC_NAME cdt_opt_app + CONFIG_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cdt_opt_test_config.json GIT_REPOSITORY "https://github.com/wildmeshing/data.git" GIT_TAG e39994334263bd3fa8f6eb026d8487b5e9ae7191) diff --git a/applications/cdt_sec/cdt_sec_main.cpp b/applications/cdt_opt/cdt_opt_main.cpp similarity index 97% rename from applications/cdt_sec/cdt_sec_main.cpp rename to applications/cdt_opt/cdt_opt_main.cpp index 481aa0ab6d..b04970c248 100644 --- a/applications/cdt_sec/cdt_sec_main.cpp +++ b/applications/cdt_opt/cdt_opt_main.cpp @@ -20,7 +20,7 @@ #include -#include "cdt_sec_spec.hpp" +#include "cdt_opt_spec.hpp" using namespace wmtk; namespace fs = std::filesystem; @@ -45,12 +45,12 @@ int main(int argc, char* argv[]) j = nlohmann::json::parse(ifs); jse::JSE spec_engine; - bool r = spec_engine.verify_json(j, cdt_sec_spec); + bool r = spec_engine.verify_json(j, cdt_opt_spec); if (!r) { wmtk::logger().error("{}", spec_engine.log2str()); return 1; } else { - j = spec_engine.inject_defaults(j, cdt_sec_spec); + j = spec_engine.inject_defaults(j, cdt_opt_spec); } } diff --git a/applications/cdt_sec/cdt_sec_spec.hpp b/applications/cdt_opt/cdt_opt_spec.hpp similarity index 95% rename from applications/cdt_sec/cdt_sec_spec.hpp rename to applications/cdt_opt/cdt_opt_spec.hpp index b93d1d4589..addca1e501 100644 --- a/applications/cdt_sec/cdt_sec_spec.hpp +++ b/applications/cdt_opt/cdt_opt_spec.hpp @@ -2,7 +2,7 @@ #include namespace { -nlohmann::json cdt_sec_spec = R"( +nlohmann::json cdt_opt_spec = R"( [ { "pointer": "/", diff --git a/applications/cdt_sec/cdt_sec_test_config.json b/applications/cdt_opt/cdt_opt_test_config.json similarity index 87% rename from applications/cdt_sec/cdt_sec_test_config.json rename to applications/cdt_opt/cdt_opt_test_config.json index 0cfade5b55..625cce2164 100644 --- a/applications/cdt_sec/cdt_sec_test_config.json +++ b/applications/cdt_opt/cdt_opt_test_config.json @@ -2,7 +2,7 @@ "test_directory": "unit_test", "tests": [], "release_only_tests": [ - "cdt_sec_out.json" + "cdt_opt_out.json" ], "input_tag": "input", "oracle_tag": "report", From 61543351adb930c886e272b3a03d3c958c90c9f4 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 17:25:28 -0500 Subject: [PATCH 09/12] rename --- applications/cdt_opt/cdt_opt_main.cpp | 24 ++---------------------- 1 file changed, 2 insertions(+), 22 deletions(-) diff --git a/applications/cdt_opt/cdt_opt_main.cpp b/applications/cdt_opt/cdt_opt_main.cpp index b04970c248..5f748248c3 100644 --- a/applications/cdt_opt/cdt_opt_main.cpp +++ b/applications/cdt_opt/cdt_opt_main.cpp @@ -89,34 +89,14 @@ int main(int argc, char* argv[]) } std::string output_file = j["output"]; - wmtk::components::output::output(*parent_mesh, output_file + "_before_sec", "vertices"); - wmtk::components::output::output(*child_mesh, output_file + "_surface_before_sec", "vertices"); + wmtk::components::output::output(*parent_mesh, output_file + "_before_opt", "vertices"); + wmtk::components::output::output(*child_mesh, output_file + "_surface_before_opt", "vertices"); std::vector pass_through; auto boundary_handle = parent_mesh->get_attribute_handle("is_boundary", PrimitiveType::Triangle); pass_through.push_back(boundary_handle); - // auto child_mesh_position_handle = - // child_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); - - // attribute::MeshAttributeHandle parent_mesh_position_handle = - // parent_mesh->get_attribute_handle("vertices", PrimitiveType::Vertex); - - // { - // using namespace components::shortest_edge_collapse; - - // ShortestEdgeCollapseOptions options; - // options.position_handle = child_mesh_position_handle; - // options.other_position_handles.emplace_back(parent_mesh_position_handle); - // options.length_rel = j["length_rel"]; - // options.envelope_size = j["envelope_size"]; - // options.check_inversions = true; - // options.pass_through_attributes = pass_through; - - // shortest_edge_collapse(*parent_mesh, options); - // } - std::vector enves; wmtk::components::EnvelopeOptions e_surface; From 978d1de984408dec13508534e77daede7e3b5af0 Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 17:34:26 -0500 Subject: [PATCH 10/12] update cdt unit test tag --- applications/cdt_opt/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/cdt_opt/CMakeLists.txt b/applications/cdt_opt/CMakeLists.txt index 234c2bed64..8b89c40771 100644 --- a/applications/cdt_opt/CMakeLists.txt +++ b/applications/cdt_opt/CMakeLists.txt @@ -17,4 +17,4 @@ wmtk::output) wmtk_register_integration_test(EXEC_NAME cdt_opt_app CONFIG_FILE ${CMAKE_CURRENT_SOURCE_DIR}/cdt_opt_test_config.json GIT_REPOSITORY "https://github.com/wildmeshing/data.git" - GIT_TAG e39994334263bd3fa8f6eb026d8487b5e9ae7191) + GIT_TAG d1063ed50de45a1bcae9f8b6ae9b8b1d42885abe) From 556dd84203e80071f32a9b8cd6384e63a1756fff Mon Sep 17 00:00:00 2001 From: JcDai Date: Mon, 4 Nov 2024 20:03:55 -0500 Subject: [PATCH 11/12] add options to cdt opt --- applications/cdt_opt/cdt_opt_main.cpp | 4 ++-- applications/cdt_opt/cdt_opt_spec.hpp | 14 +++++++++++++- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/applications/cdt_opt/cdt_opt_main.cpp b/applications/cdt_opt/cdt_opt_main.cpp index 5f748248c3..3dde898991 100644 --- a/applications/cdt_opt/cdt_opt_main.cpp +++ b/applications/cdt_opt/cdt_opt_main.cpp @@ -113,8 +113,8 @@ int main(int argc, char* argv[]) wmo.input_mesh = parent_mesh; wmo.input_mesh_position = "vertices"; wmo.target_edge_length = j["length_rel"]; - wmo.target_max_amips = 50; - wmo.max_passes = 10; + wmo.target_max_amips = j["target_max_amips"]; + wmo.max_passes = j["max_passes"]; wmo.intermediate_output = false; wmo.replace_double_coordinate = false; wmo.scheduler_update_frequency = 0; diff --git a/applications/cdt_opt/cdt_opt_spec.hpp b/applications/cdt_opt/cdt_opt_spec.hpp index addca1e501..62b95733c4 100644 --- a/applications/cdt_opt/cdt_opt_spec.hpp +++ b/applications/cdt_opt/cdt_opt_spec.hpp @@ -12,7 +12,9 @@ nlohmann::json cdt_opt_spec = R"( "root", "report", "envelope_size", - "length_rel" + "length_rel", + "target_max_amips", + "max_passes" ] }, { @@ -29,6 +31,16 @@ nlohmann::json cdt_opt_spec = R"( "type": "float", "default": 0.1 }, +{ + "pointer": "/target_max_amips", + "type": "float", + "default": 50 +}, +{ + "pointer": "/max_passes", + "type": "int", + "default": 10 +}, { "pointer": "/output", "type": "string" From ef853602d5ab861116fb6a35da03daa2a7c67964 Mon Sep 17 00:00:00 2001 From: JcDai Date: Wed, 6 Nov 2024 16:07:00 -0500 Subject: [PATCH 12/12] remove cdt_optimization component --- applications/cdt_opt/CMakeLists.txt | 2 - applications/cdt_opt/cdt_opt_main.cpp | 2 - components/CMakeLists.txt | 1 - .../cdt_optimization/CMakeLists.txt | 13 - .../cdt_optimization/cdt_optimization.cpp | 425 ------------------ .../cdt_optimization/cdt_optimization.hpp | 18 - 6 files changed, 461 deletions(-) delete mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt delete mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp delete mode 100644 components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp diff --git a/applications/cdt_opt/CMakeLists.txt b/applications/cdt_opt/CMakeLists.txt index 8b89c40771..6c4d04a807 100644 --- a/applications/cdt_opt/CMakeLists.txt +++ b/applications/cdt_opt/CMakeLists.txt @@ -8,8 +8,6 @@ wmtk_add_application(cdt_opt_app target_link_libraries(cdt_opt_app PRIVATE wmtk::input wmtk::multimesh -wmtk::cdt_optimization -wmtk::shortest_edge_collapse wmtk::CDT wmtk::wildmeshing wmtk::output) diff --git a/applications/cdt_opt/cdt_opt_main.cpp b/applications/cdt_opt/cdt_opt_main.cpp index 3dde898991..e592c0a5f6 100644 --- a/applications/cdt_opt/cdt_opt_main.cpp +++ b/applications/cdt_opt/cdt_opt_main.cpp @@ -12,11 +12,9 @@ #include -#include #include #include #include -#include #include diff --git a/components/CMakeLists.txt b/components/CMakeLists.txt index 739d12bbc8..c749617482 100644 --- a/components/CMakeLists.txt +++ b/components/CMakeLists.txt @@ -34,7 +34,6 @@ add_subdirectory("wildmeshing/wmtk/components/wildmeshing") add_subdirectory("shortest_edge_collapse/wmtk/components/shortest_edge_collapse") add_subdirectory("longest_edge_split/wmtk/components/longest_edge_split") add_subdirectory("CDT/wmtk/components/CDT") -add_subdirectory("cdt_optimization/wmtk/components/cdt_optimization") add_subdirectory("edge_insertion/wmtk/components/edge_insertion") diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt b/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt deleted file mode 100644 index 68e0434734..0000000000 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -set(COMPONENT_NAME cdt_optimization) -add_component(${COMPONENT_NAME}) -if(NOT ${WMTK_ENABLE_COMPONENT_${COMPONENT_NAME}}) - return() -endif() - -set(SRC_FILES - cdt_optimization.hpp - cdt_optimization.cpp - ) - - -target_sources(wmtk_${COMPONENT_NAME} PRIVATE ${SRC_FILES}) diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp deleted file mode 100644 index e16d64cd82..0000000000 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.cpp +++ /dev/null @@ -1,425 +0,0 @@ -#include "cdt_optimization.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 "predicates.h" - - -namespace wmtk::components { -std::shared_ptr cdt_optimization( - TriMesh& mesh, - const attribute::MeshAttributeHandle& position_handle, - std::optional& inversion_position_handle, - const int64_t passes, - bool update_other_position, - const double length_rel, - bool lock_boundary, - double envelope_size, - const std::vector& pass_through) -{ - if (mesh.top_simplex_type() != PrimitiveType::Triangle) { - log_and_throw_error( - "shortest edge collapse works only for triangle meshes: {}", - primitive_type_name(mesh.top_simplex_type())); - } - - // TriMesh& mesh = static_cast(position_handle.mesh()); - - auto pass_through_attributes = pass_through; - - bool check_inversion = inversion_position_handle.has_value(); - - ///////////////////////////////////////////// - - auto visited_edge_flag = - mesh.register_attribute("visited_edge", PrimitiveType::Edge, 1, false, char(1)); - - auto update_flag_func = [](Eigen::Ref 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, - position_handle, - update_flag_func); - - ////////////////////////////////// - // 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 = [](Eigen::Ref P) -> Eigen::VectorXd { - assert(P.cols() == 2); - assert(P.rows() == 2 || P.rows() == 3); - return Eigen::VectorXd::Constant(1, (P.col(0) - P.col(1)).norm()); - }; - auto edge_length_update = - std::make_shared>( - edge_length_attribute, - position_handle, - compute_edge_length); - edge_length_update->run_on_all(); - - ////////////////////////////////// - // computing bbox diagonal - Eigen::VectorXd bmin(mesh.top_cell_dimension()); - bmin.setConstant(std::numeric_limits::max()); - Eigen::VectorXd bmax(mesh.top_cell_dimension()); - bmax.setConstant(std::numeric_limits::lowest()); - - auto pt_accessor = mesh.create_const_accessor(position_handle); - - const auto vertices = mesh.get_all(PrimitiveType::Vertex); - for (const auto& v : vertices) { - const auto p = pt_accessor.vector_attribute(v); - for (int64_t d = 0; d < bmax.size(); ++d) { - bmin[d] = std::min(bmin[d], p[d]); - bmax[d] = std::max(bmax[d], p[d]); - } - } - - const double bbdiag = (bmax - bmin).norm(); - - const double length_abs = bbdiag * length_rel; - - wmtk::logger().info( - "bbox max {}, bbox min {}, diag {}, target edge length {}", - bmax, - bmin, - bbdiag, - length_abs); - - auto short_edges_first_priority = [&](const simplex::Simplex& s) { - assert(s.primitive_type() == PrimitiveType::Edge); - return edge_length_accessor.const_scalar_attribute(s.tuple()); - }; - auto long_edges_first_priority = [&](const simplex::Simplex& s) { - assert(s.primitive_type() == PrimitiveType::Edge); - return -edge_length_accessor.const_scalar_attribute(s.tuple()); - }; - pass_through_attributes.push_back(edge_length_attribute); - auto todo_smaller = std::make_shared( - mesh, - edge_length_attribute.as(), - 4. / 5. * length_abs); // MTAO: why is this 4/5? - - auto todo_larger = std::make_shared( - mesh, - edge_length_attribute.as(), - 4.0 / 3.0 * length_abs); - - //////////////////////////invariants - - auto invariant_link_condition = std::make_shared(mesh); - - auto invariant_interior_edge = std::make_shared(mesh); - auto invariant_interior_vertex = std::make_shared(mesh); - - auto set_all_invariants = [&](auto&& m) { - invariant_interior_edge->add( - std::make_shared(m, PrimitiveType::Edge)); - invariant_interior_vertex->add( - std::make_shared(m, PrimitiveType::Vertex)); - }; - multimesh::MultiMeshVisitor visitor(set_all_invariants); - visitor.execute_from_root(mesh); - - auto invariant_mm_map = std::make_shared(mesh); - - ////////////// positions - std::vector positions; - positions.push_back(position_handle); - if (check_inversion) { - positions.push_back(inversion_position_handle.value()); - } - - ////////////////////////////////////////// - // split - ////////////////////////////////////////// - - auto split = std::make_shared(mesh); - split->add_invariant(todo_larger); - if (check_inversion) { - split->add_invariant(std::make_shared>( - inversion_position_handle.value().mesh(), - inversion_position_handle.value().as())); - } - - split->set_priority(long_edges_first_priority); - - split->set_new_attribute_strategy( - visited_edge_flag, - wmtk::operations::SplitBasicStrategy::None, - wmtk::operations::SplitRibBasicStrategy::None); - - split->add_transfer_strategy(tag_update); - split->add_transfer_strategy(edge_length_update); - - for (auto& p : positions) { - split->set_new_attribute_strategy(p); - } - - for (const auto& attr : pass_through_attributes) { - split->set_new_attribute_strategy( - attr, - wmtk::operations::SplitBasicStrategy::None, - wmtk::operations::SplitRibBasicStrategy::None); - } - - - ////////////////////////////////////////// - // collapse - ////////////////////////////////////////// - auto collapse = std::make_shared(mesh); - collapse->add_invariant(todo_smaller); - collapse->add_invariant(invariant_link_condition); - collapse->add_invariant(invariant_mm_map); - - - if (envelope_size > 0) { - collapse->add_invariant(std::make_shared( - position_handle, - bbdiag * envelope_size, - position_handle)); - } - - - if (check_inversion) { - collapse->add_invariant(std::make_shared>( - inversion_position_handle.value().mesh(), - inversion_position_handle.value().as())); - } - collapse->set_new_attribute_strategy( - visited_edge_flag, - wmtk::operations::CollapseBasicStrategy::None); - collapse->add_transfer_strategy(tag_update); - - collapse->set_priority(short_edges_first_priority); - collapse->add_transfer_strategy(edge_length_update); - - if (lock_boundary) { - collapse->add_invariant(invariant_interior_edge); - // set collapse towards boundary - for (auto& p : positions) { - auto tmp = std::make_shared>(p); - tmp->set_strategy(wmtk::operations::CollapseBasicStrategy::CopyOther); - tmp->set_simplex_predicate(wmtk::operations::BasicSimplexPredicate::IsInterior); - collapse->set_new_attribute_strategy(p, tmp); - } - } else { - for (auto& p : positions) { - collapse->set_new_attribute_strategy( - p, - wmtk::operations::CollapseBasicStrategy::CopyOther); - } - } - - - for (const auto& attr : pass_through_attributes) { - collapse->set_new_attribute_strategy(attr); - } - - ////////////////////////////////////////// - // swap - ////////////////////////////////////////// - - std::shared_ptr amips = - std::make_shared(mesh, position_handle); - auto function_invariant = - std::make_shared(mesh.top_simplex_type(), amips); - - auto swap = std::make_shared(mesh); - - swap->add_invariant( - std::make_shared(mesh, position_handle.as())); - - swap->add_invariant(invariant_mm_map); - - swap->collapse().add_invariant(invariant_link_condition); - - if (envelope_size > 0) { - swap->add_invariant(std::make_shared( - position_handle, - bbdiag * envelope_size, - position_handle)); - } - - - if (check_inversion) { - swap->collapse().add_invariant(std::make_shared>( - inversion_position_handle.value().mesh(), - inversion_position_handle.value().as())); - } - - // swap->add_invariant(function_invariant); - - swap->split().set_new_attribute_strategy( - visited_edge_flag, - wmtk::operations::SplitBasicStrategy::None, - wmtk::operations::SplitRibBasicStrategy::None); - swap->collapse().set_new_attribute_strategy( - visited_edge_flag, - wmtk::operations::CollapseBasicStrategy::None); - - swap->add_transfer_strategy(tag_update); - - swap->add_transfer_strategy(edge_length_update); - - for (auto& p : positions) { - swap->split().set_new_attribute_strategy(p); - swap->collapse().set_new_attribute_strategy( - p, - wmtk::operations::CollapseBasicStrategy::CopyOther); - } - - swap->set_priority(short_edges_first_priority); - - swap->add_transfer_strategy(edge_length_update); - - for (const auto& attr : pass_through_attributes) { - swap->split().set_new_attribute_strategy(attr); - swap->collapse().set_new_attribute_strategy(attr); - } - - ////////////////////////////////////////// - // scheduler - ////////////////////////////////////////// - - Scheduler scheduler; - - auto collapse_stats_pre = - scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); - - logger().info( - "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " - "executing: {}", - "collapse", - collapse_stats_pre.number_of_performed_operations(), - collapse_stats_pre.number_of_successful_operations(), - collapse_stats_pre.number_of_failed_operations(), - collapse_stats_pre.collecting_time, - collapse_stats_pre.sorting_time, - collapse_stats_pre.executing_time); - - for (int64_t i = 0; i < passes; ++i) { - auto split_stats = scheduler.run_operation_on_all(*split, visited_edge_flag.as()); - - logger().info( - "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " - "executing: {}", - "split", - split_stats.number_of_performed_operations(), - split_stats.number_of_successful_operations(), - split_stats.number_of_failed_operations(), - split_stats.collecting_time, - split_stats.sorting_time, - split_stats.executing_time); - - auto collapse_stats = - scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); - - logger().info( - "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " - "executing: {}", - "collapse", - collapse_stats.number_of_performed_operations(), - collapse_stats.number_of_successful_operations(), - collapse_stats.number_of_failed_operations(), - collapse_stats.collecting_time, - collapse_stats.sorting_time, - collapse_stats.executing_time); - - auto swap_stats = scheduler.run_operation_on_all(*swap, visited_edge_flag.as()); - - logger().info( - "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " - "executing: {}", - "swap", - swap_stats.number_of_performed_operations(), - swap_stats.number_of_successful_operations(), - swap_stats.number_of_failed_operations(), - swap_stats.collecting_time, - swap_stats.sorting_time, - swap_stats.executing_time); - - multimesh::consolidate(mesh); - } - - auto collapse_stats_after = - scheduler.run_operation_on_all(*collapse, visited_edge_flag.as()); - - logger().info( - "Executed {}, {} ops (S/F) {}/{}. Time: collecting: {}, sorting: {}, " - "executing: {}", - "collapse", - collapse_stats_after.number_of_performed_operations(), - collapse_stats_after.number_of_successful_operations(), - collapse_stats_after.number_of_failed_operations(), - collapse_stats_after.collecting_time, - collapse_stats_after.sorting_time, - collapse_stats_after.executing_time); - - multimesh::consolidate(mesh); - - // debug code - for (const auto& e : mesh.get_all(PrimitiveType::Edge)) { - if (mesh.is_boundary(PrimitiveType::Edge, e)) { - logger().error("after sec mesh has nonmanifold edges"); - } - } - - auto inversion_pos_accessor = - inversion_position_handle.value().mesh().create_const_accessor( - inversion_position_handle.value()); - - // debug code - for (const auto& t : - inversion_position_handle.value().mesh().get_all(PrimitiveType::Tetrahedron)) { - const auto vertices = inversion_position_handle.value().mesh().orient_vertices(t); - std::vector pos; - for (int i = 0; i < vertices.size(); ++i) { - pos.push_back(inversion_pos_accessor.const_vector_attribute(vertices[i])); - } - if (orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data()) <= 0) { - wmtk::logger().error( - "Flipped or degenerated tet! orient3d = {}", - orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data())); - } else if (orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data()) <= 1e-5) { - wmtk::logger().error( - "Nearly degenerated tet! orient3d = {}", - orient3d(pos[0].data(), pos[1].data(), pos[2].data(), pos[3].data())); - } - } - - - return mesh.shared_from_this(); -} -} // namespace wmtk::components diff --git a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp b/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp deleted file mode 100644 index f34727db50..0000000000 --- a/components/cdt_optimization/wmtk/components/cdt_optimization/cdt_optimization.hpp +++ /dev/null @@ -1,18 +0,0 @@ -#pragma once -#include -#include - -namespace wmtk::components { - -std::shared_ptr cdt_optimization( - TriMesh& mesh, - const attribute::MeshAttributeHandle& position_handle, - std::optional& inversion_position_handle, - const int64_t passes, - bool update_other_position, - const double length_rel, - bool lock_boundary, - double envelope_size, - const std::vector& pass_through); - -} // namespace wmtk::components