From e0461bb87878f71c74706adfa0d70686e7478117 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 20 May 2026 13:44:24 +0200 Subject: [PATCH 1/3] Remove IPC toolkit. --- CMakeLists.txt | 2 - cmake/recipes/ipc_toolkit.cmake | 20 - .../image_simulation/ImageSimulationMesh.h | 52 +-- .../ImageSimulationMeshTri.cpp | 344 +----------------- .../ImageSimulationMeshTri.hpp | 73 +--- .../components/image_simulation/Parameters.h | 9 +- .../components/image_simulation/Smooth.cpp | 329 ----------------- .../image_simulation/image_simulation.cpp | 4 - .../image_simulation_spec.hpp | 18 +- .../image_simulation_spec.json | 18 +- .../tests/test_image_simulation_energies.cpp | 72 +--- pyproject.toml | 6 +- src/wmtk/optimization/BarrierEnergy.cpp | 191 ---------- src/wmtk/optimization/BarrierEnergy.hpp | 131 ------- tests/test_optimization.cpp | 204 +---------- 15 files changed, 24 insertions(+), 1449 deletions(-) delete mode 100644 cmake/recipes/ipc_toolkit.cmake delete mode 100644 src/wmtk/optimization/BarrierEnergy.cpp delete mode 100644 src/wmtk/optimization/BarrierEnergy.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bd40c119b6..a16ae7c3f2 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -82,7 +82,6 @@ include(paraviewo) include(fenvelope) include(nanoflann) include(polysolve) -include(ipc_toolkit) include(simple_bvh) # Sort projects inside the solution @@ -142,7 +141,6 @@ target_link_libraries(wildmeshing_toolkit PUBLIC FastEnvelope::FastEnvelope simple_bvh::simple_bvh polysolve::polysolve - ipc::toolkit ) if(MSVC) diff --git a/cmake/recipes/ipc_toolkit.cmake b/cmake/recipes/ipc_toolkit.cmake deleted file mode 100644 index 0f81e95b3f..0000000000 --- a/cmake/recipes/ipc_toolkit.cmake +++ /dev/null @@ -1,20 +0,0 @@ -# IPC Toolkit (https://github.com/ipc-sim/ipc-toolkit) -# License: MIT - -if(TARGET ipc::toolkit) - return() -endif() - -message(STATUS "Third-party: creating target 'ipc::toolkit'") - -include(CPM) -CPMAddPackage( - NAME ipc_toolkit - GITHUB_REPOSITORY ipc-sim/ipc-toolkit - GIT_TAG e94c11d174de8f9805fcebc3b0cfe3b9f1306ab5 - OPTIONS - "IPC_TOOLKIT_WITH_CUDA OFF" - "IPC_TOOLKIT_WITH_SIMD OFF" -) - -set_target_properties(ipc_toolkit PROPERTIES FOLDER third_party) \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h index 20da84ea06..b979780503 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMesh.h @@ -144,14 +144,10 @@ class ImageSimulationMesh : public wmtk::TetMesh std::shared_ptr m_order_2_edge_envelope; // todo: add sample envelope option tbb::enumerable_thread_specific> m_solver; - Eigen::SparseMatrix m_surface_mass; // the mass matrix for surface vertices - Eigen::SparseMatrix m_surface_stiffness; // stiffness matrix for surface vertices // scaling factors double m_s_amips = -1; - double m_s_smooth = -1; double m_s_envelope = -1; - double m_s_barrier = -1; // When set, split_edge_after binary-searches vmid onto the zero-crossing of this function. // Negative = stays on v1 side, positive = stays on v2 side. @@ -173,39 +169,12 @@ class ImageSimulationMesh : public wmtk::TetMesh optimization::deactivate_opt_logger(); m_s_amips = 1.; - /** - * TODO - */ - m_s_smooth = 1.; - /** - * TODO - */ m_s_envelope = 1. / (m_params.diag_l * m_params.eps * m_params.eps); - /** - * TODO - */ - m_s_barrier = 1. / (m_params.diag_l4); - - // check weights (ignoring barrier here) - { - double& wa = m_params.w_amips; - double& ws = m_params.w_smooth; - const double sum = wa + ws; - if (sum > 1) { - wa /= sum; - ws /= sum; - logger().warn( - "Weights for AMIPS and smooth sum up to greater than 1. Rescaling to \n " - "w_amips = {}, \n w_smooth = {}", - wa, - ws); - } - double& we = m_params.w_envelope; - we = 1 - (wa + ws); - logger().info("w_envelope = {}", we); - } - init_separation_weight(); + double& wa = m_params.w_amips; + double& we = m_params.w_envelope; + we = 1 - wa; + logger().info("w_envelope = {}", we); } ~ImageSimulationMesh() {} @@ -354,8 +323,6 @@ class ImageSimulationMesh : public wmtk::TetMesh void init_envelope(const MatrixXd& V, const MatrixXi& F); - void init_separation_weight(); - CellTag string_set_to_cell_tag(const std::set& str_set); double get_length2(const Tuple& l) const; @@ -424,21 +391,10 @@ class ImageSimulationMesh : public wmtk::TetMesh double get_quality(const std::array& vs) const; double get_quality(const Tuple& loc) const; - void build_mass_matrix(); - std::vector> get_amips_assembles(const Tuple& t) const; std::shared_ptr get_amips_energy(const Tuple& t) const; - - std::shared_ptr get_smooth_energy(const Tuple& t) const; std::shared_ptr get_envelope_energy(const Tuple& t) const; - std::shared_ptr get_barrier_energy( - const Tuple& t, - const bool use_full_surface = false) const; - - void substructure_region(const Tuple& v_tuple, MatrixXd& V, MatrixXi& E, size_t& vid_local) - const; - /** * @brief Round a vertex position to floating point. diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 22f6916fac..90f09fb60e 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -16,7 +16,6 @@ #include #include #include -#include #include #include #include @@ -350,52 +349,6 @@ void ImageSimulationMeshTri::init_envelope(const MatrixXd& V, const MatrixXi& E) } } -void ImageSimulationMeshTri::init_separation_weight() -{ - if (m_params.separation_factor <= 0) { - m_params.w_separate = 0; - logger().warn( - "Separation factor is {} -> No separation. IPC barrier weight = {}", - m_params.separation_factor, - m_params.w_separate); - return; - } - //m_params.w_separate = m_params.w_envelope * (16. / 5.) * m_params.dhat; - - /** - * This assumes that the barrier term is the clamped log barrier term: - * E_B = b(d) = -(d-\hat{d})^2\ln\left(\frac{d}{\hat{d}}\right) - * - * E_B' = b'(\hat{d}/2) = \hat{d} * (\ln(2) - 0.5) - * - * barrier scaling: s_B = 1/l^4 - * envelope scaling: s_E = 1 / (l * eps^2) - * - * l is the bounding box diagonal. - * - * The envelope energy E_E = \Vert p - p_0 \Vert^2 - * E_E'(\hat{d}/2) = \hat{d} - * - * We compute w_E from - * s_B w_B E_B' = s_E w_E E_E' - * - * w_B = w_E l^3 / (eps^2 * (\ln(2) - 0.5)) - */ - - const double l3 = m_params.diag_l3; - const double eps2 = m_params.eps * m_params.eps; - const double ln2 = std::log(2.0); - - m_params.w_separate = m_params.w_envelope * l3 / (eps2 * (ln2 - 0.5)); - - logger().info("IPC barrier weight = {}, dhat = {}", m_params.w_separate, m_params.dhat); - - if (m_params.w_separate < 0) { - logger().error("Negative separation weight: {}. Setting weight to 0.", m_params.w_separate); - m_params.w_separate = 0; - } -} - CellTag ImageSimulationMeshTri::string_set_to_cell_tag(const std::set& str_set) { CellTag cell_tag; @@ -770,10 +723,6 @@ void ImageSimulationMeshTri::write_vtu_with_energies(const std::string& path) co VectorXd v_sizing_field(vert_capacity()); v_sizing_field.setZero(); - MatrixXd v_energy_grad_barrier(vert_capacity(), 2); - v_energy_grad_barrier.setZero(); - MatrixXd v_energy_grad_smooth(vert_capacity(), 2); - v_energy_grad_smooth.setZero(); MatrixXd v_energy_grad_amips(vert_capacity(), 2); v_energy_grad_amips.setZero(); MatrixXd v_energy_grad_envelope(vert_capacity(), 2); @@ -815,28 +764,16 @@ void ImageSimulationMeshTri::write_vtu_with_energies(const std::string& path) co auto energy_sum = std::make_shared(); auto amips_energy = get_amips_energy(v); - auto smooth_energy = get_smooth_energy(v); auto envelope_energy = get_envelope_energy(v); - // auto barrier_energy = get_barrier_energy(v); - - if (!smooth_energy) { - continue; - } energy_sum->add_energy(amips_energy); - energy_sum->add_energy(smooth_energy); energy_sum->add_energy(envelope_energy); - // energy_sum->add_energy(barrier_energy); VectorXd g; amips_energy->gradient(x, g); v_energy_grad_amips.row(vid) = g; - smooth_energy->gradient(x, g); - v_energy_grad_smooth.row(vid) = g; envelope_energy->gradient(x, g); v_energy_grad_envelope.row(vid) = g; - // barrier_energy->gradient(x, g); - // v_energy_grad_barrier.row(vid) = g; energy_sum->gradient(x, g); v_energy_grad_sum.row(vid) = g; } @@ -859,9 +796,7 @@ void ImageSimulationMeshTri::write_vtu_with_energies(const std::string& path) co surf_writer = std::make_shared(); surf_writer->add_field("sizing_field", v_sizing_field); surf_writer->add_field("amips_grad", v_energy_grad_amips); - surf_writer->add_field("smooth_grad", v_energy_grad_smooth); surf_writer->add_field("envelope_grad", v_energy_grad_envelope); - surf_writer->add_field("barrier_grad", v_energy_grad_barrier); surf_writer->add_field("sum_grad", v_energy_grad_sum); @@ -1528,13 +1463,6 @@ bool ImageSimulationMeshTri::swap_edge_after(const Tuple& t) void ImageSimulationMeshTri::smooth_all_vertices(const size_t n_iters) { - /** - * Mass matrix is only necessary if w_smooth is positive but we still compute the mass matrix - * all the time. It's rather cheap to do and allows for printing all the energy terms. - */ - build_mass_matrix(); - - // actual smoothing for (size_t i = 0; i < n_iters; ++i) { // log_total_surface_energy(); igl::Timer timer; @@ -1638,27 +1566,14 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) if (VA[vid].m_is_on_surface) { assert(!surf_assembles.empty()); - if (surf_assembles.size() != 3) { - return false; - } - auto energy_sum = std::make_shared(); - if (m_params.w_smooth > 0) { - auto smooth_energy = get_smooth_energy(t); - assert(smooth_energy); - energy_sum->add_energy(smooth_energy); - } - if (m_params.w_amips > 0) { energy_sum->add_energy(amips_energy); } if (m_params.w_envelope > 0) { energy_sum->add_energy(get_envelope_energy(t)); } - // if (m_params.w_separate > 0) { - // energy_sum->add_energy(get_barrier_energy(t)); - // } total_energy = energy_sum; } else { @@ -1711,36 +1626,6 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) return true; } -void ImageSimulationMeshTri::build_mass_matrix() -{ - m_surface_mass.resize(vert_capacity()); - m_surface_stiffness.resize(vert_capacity()); - for (const Tuple& t : get_vertices()) { - const size_t vid = t.vid(*this); - if (!m_vertex_attribute.at(vid).m_is_on_surface) { - continue; - } - const auto es = get_surface_edges_for_vertex(vid); - if (es.size() != 2) { - continue; - } - auto& M = m_surface_mass[vid]; - auto& L_w = m_surface_stiffness[vid]; - - std::array pts; - pts[0] = m_vertex_attribute.at(vid).m_pos; - - for (size_t i = 0; i < 2; ++i) { - const auto& vs = es.edges()[i].vertices(); - size_t neighbor_id = vs[0] != vid ? vs[0] : vs[1]; - pts[i + 1] = m_vertex_attribute.at(neighbor_id).m_pos; - } - - optimization::BiharmonicEnergy2D::local_mass_and_stiffness(pts, M, L_w); - // optimization::SmoothingEnergy2D::uniform_mass_and_stiffness(pts, M, L_w); - } -} - std::vector ImageSimulationMeshTri::get_surface_assembles(const Tuple& t) const { const auto& VA = m_vertex_attribute; @@ -1766,37 +1651,11 @@ std::vector ImageSimulationMeshTri::get_surface_assembles(const Tuple& return surface_pts; } -std::shared_ptr ImageSimulationMeshTri::get_smooth_energy( - const Tuple& t) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid = t.vid(*this); - - auto surf_assembles_vec = get_surface_assembles(t); - if (surf_assembles_vec.size() != 3) { - return nullptr; - } - - // mass and stiffness - const auto& M = m_surface_mass[vid]; - const auto& L_w = m_surface_stiffness[vid]; - - const double w = m_s_smooth * m_params.w_smooth; - - std::array pts{ - surf_assembles_vec[0], - surf_assembles_vec[1], - surf_assembles_vec[2]}; - auto smooth_energy = std::make_shared(pts, M, L_w, w); - - return smooth_energy; -} std::shared_ptr ImageSimulationMeshTri::get_envelope_energy( const Tuple& t) const { - const auto& M = m_surface_mass[t.vid(*this)]; - const double w = m_s_envelope * M * m_params.w_envelope; + const double w = m_s_envelope * m_params.w_envelope; auto envelope_energy = std::make_shared( m_envelope_orig, @@ -1805,74 +1664,6 @@ std::shared_ptr ImageSimulationMeshTri::get_envel return envelope_energy; } -std::shared_ptr ImageSimulationMeshTri::get_barrier_energy( - const Tuple& t, - const bool use_full_surface) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid = t.vid(*this); - - if (!VA[vid].m_is_on_surface) { - return nullptr; - } - - const auto& M = m_surface_mass[t.vid(*this)]; - if (M <= 0) { - return nullptr; - } - - // barrier energy - MatrixXd V_barrier; - MatrixXi E_barrier; - size_t vid_barrier; - if (!use_full_surface) { - substructure_region(t, V_barrier, E_barrier, vid_barrier); - } else { - // this is horribly expensive and is only here for testing - - std::vector surf_points; - std::vector global_to_local_vid_map(vert_capacity()); - for (size_t i = 0; i < vert_capacity(); ++i) { - const Tuple v = tuple_from_vertex(i); - if (!v.is_valid(*this)) { - continue; - } - const size_t _vid = v.vid(*this); - if (!m_vertex_attribute.at(_vid).m_is_on_surface) { - continue; - } - surf_points.push_back(m_vertex_attribute.at(_vid).m_pos); - global_to_local_vid_map[_vid] = surf_points.size() - 1; - } - - V_barrier.resize(surf_points.size(), 2); - for (size_t i = 0; i < surf_points.size(); ++i) { - V_barrier.row(i) = surf_points[i]; - } - - const auto surf_edges = get_edges_by_condition([](auto& f) { return f.m_is_surface_fs; }); - E_barrier.resize(surf_edges.size(), 2); - for (size_t i = 0; i < surf_edges.size(); ++i) { - const size_t v0 = global_to_local_vid_map[surf_edges[i][0]]; - const size_t v1 = global_to_local_vid_map[surf_edges[i][1]]; - E_barrier.row(i) = Vector2i(v0, v1); - } - - vid_barrier = global_to_local_vid_map[vid]; - } - - // The inverse mass comes from the weight compuation based on the envelope energy. This contains - // the mass and therefore the barrier energy must contain the inverse. - auto barrier_energy = std::make_shared( - V_barrier, - E_barrier, - vid_barrier, - m_params.dhat, - m_s_barrier * m_params.w_separate / M); - - return barrier_energy; -} - std::vector> ImageSimulationMeshTri::get_amips_assembles(const Tuple& t) const { const size_t vid = t.vid(*this); @@ -1927,13 +1718,10 @@ std::shared_ptr ImageSimulationMeshTri::get_amips void ImageSimulationMeshTri::log_total_surface_energy() { - build_mass_matrix(); - double e_sum = 0; double e_amips = 0; double e_smooth = 0; double e_envelope = 0; - double e_barrier = 0; size_t n_pts = 0; const auto& VA = m_vertex_attribute; @@ -1946,38 +1734,25 @@ void ImageSimulationMeshTri::log_total_surface_energy() const Vector2d old_pos = VA[vid].m_pos; auto amips_energy = get_amips_energy(t); - auto smooth_energy = get_smooth_energy(t); - - if (!smooth_energy) { - continue; - } auto envelope_energy = get_envelope_energy(t); - // auto barrier_energy = get_barrier_energy(t); auto energy_sum = std::make_shared(); energy_sum->add_energy(amips_energy); - energy_sum->add_energy(smooth_energy); energy_sum->add_energy(envelope_energy); - // energy_sum->add_energy(barrier_energy); e_sum += energy_sum->value(old_pos); e_amips += amips_energy->value(old_pos); - e_smooth += smooth_energy->value(old_pos); e_envelope += envelope_energy->value(old_pos); - // e_barrier += barrier_energy->value(old_pos); ++n_pts; } logger().warn( - ">>>>> Energies <<<<<\nSUM = {}\n AMIPS = {}\n Smooth = {}\n Envelope = {}\n Barrier = " - "{}\n #V = {}", + ">>>>> Energies <<<<<\nSUM = {}\n AMIPS = {}\n Envelope = {}\n #V = {}", e_sum, e_amips, - e_smooth, e_envelope, - e_barrier, n_pts); } @@ -2264,119 +2039,4 @@ std::tuple ImageSimulationMeshTri::get_max_avg_energy() return std::make_tuple(max_energy, avg_energy); } -void ImageSimulationMeshTri::substructure_region( - const Tuple& v_tuple, - MatrixXd& V, - MatrixXi& E, - size_t& vid_local) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid_global = v_tuple.vid(*this); - - if (!VA[vid_global].m_is_on_surface) { - log_and_throw_error( - "Cannot compute substructure region for vertex that is not on the surface"); - // We actually could compute this but it doesn't make sense. Remove this if-statement if we - // should ever need to compute regions for vertices not on the surface. - } - - const Vector2d& p0 = VA[vid_global].m_pos; - - double l_max = 0; // longest incident edge length - for (const Tuple& t : get_one_ring_edges_for_vertex(v_tuple)) { - const Vector2d& p1 = VA[t.vid(*this)].m_pos; - const double l_squared = (p1 - p0).squaredNorm(); - l_max = std::max(l_max, l_squared); - } - l_max = std::sqrt(l_max); - - const double r = 1.5 * m_params.dhat + l_max; // region radius - const double r2 = r * r; - - simplex::SimplexCollection candidates; // all edges within the region - - // make a BFS to find edges within distance - std::unordered_set visited; - std::queue q; - q.push(vid_global); - visited.insert(vid_global); - - while (!q.empty()) { - const size_t v_a = q.front(); - q.pop(); - - const Vector2d& p_a = VA[v_a].m_pos; - - if ((p_a - p0).squaredNorm() > r2) { - // add vertices that are connected through edges within the radius - for (const Tuple& t : get_one_ring_edges_for_vertex(v_a)) { - const size_t v_b = t.vid(*this); - const Vector2d& p_b = VA[v_b].m_pos; - const double l = (p_a - p_b).norm(); - const Vector2d p_mid = 0.5 * (p_a + p_b); - const double d = (p_mid - p0).norm() - l; - if (d > r) { - continue; // edge is too far away - } - candidates.add(simplex::Edge(v_a, v_b)); - if (visited.count(v_b) != 0) { - continue; - } - q.push(v_b); - visited.insert(v_b); - } - } else { - // no need to check edge distance as p_a is already within the radius - for (const Tuple& t : get_one_ring_edges_for_vertex(v_a)) { - const size_t v_b = t.vid(*this); - candidates.add(simplex::Edge(v_a, v_b)); - if (visited.count(v_b) != 0) { - continue; - } - q.push(v_b); - visited.insert(v_b); - } - } - } - candidates.sort_and_clean(); - - std::vector edges; // edges on the surface - for (const auto& e : candidates.edges()) { - if (is_edge_on_surface(e.vertices())) { - edges.push_back(e); - } - } - - std::set region_vertices; - for (const auto& e : edges) { - region_vertices.insert(e.vertices()[0]); - region_vertices.insert(e.vertices()[1]); - } - - // build V and E - std::vector points; - std::unordered_map global_to_local_vid_map; - for (const size_t i : region_vertices) { - const Tuple v = tuple_from_vertex(i); - assert(v.is_valid(*this)); - assert(VA[i].m_is_on_surface); - points.push_back(VA[i].m_pos); - global_to_local_vid_map[i] = points.size() - 1; - } - - V.resize(points.size(), 2); - for (size_t i = 0; i < points.size(); ++i) { - V.row(i) = points[i]; - } - - E.resize(edges.size(), 2); - for (size_t i = 0; i < edges.size(); ++i) { - const size_t v0 = global_to_local_vid_map[edges[i].vertices()[0]]; - const size_t v1 = global_to_local_vid_map[edges[i].vertices()[1]]; - E.row(i) = Vector2i(v0, v1); - } - - vid_local = global_to_local_vid_map[vid_global]; -} - } // namespace wmtk::components::image_simulation::tri \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 259c4d9e26..eb3aed6e40 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -112,14 +112,10 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool m_collapse_check_manifold = false; // manifoldness check after collapse tbb::enumerable_thread_specific> m_solver; - std::vector m_surface_mass; // the mass matrix for surface vertices - std::vector m_surface_stiffness; // stiffness matrix for surface vertices // scaling factors double m_s_amips = -1; - double m_s_smooth = -1; double m_s_envelope = -1; - double m_s_barrier = -1; ImageSimulationMeshTri(Parameters& _m_params, double envelope_eps, int _num_threads = 0) : m_params(_m_params) @@ -134,48 +130,16 @@ class ImageSimulationMeshTri : public wmtk::TriMesh m_s_amips = 1.; /** - * The bilaplacian energy for smoothing in 2D scales inverse to the length, because the - * Laplace operator times the positions scales inverse to the length. The mass scales with - * the length but as the bilaplacian contains twice the Laplace operator times positions but - * only once the mass, the entire energy scales inverse to the length. - * - * In 3D, the mass scales to length^2 and therefore the energy is dimensionless. That is not - * the case in 2D. - */ - m_s_smooth = m_params.diag_l; - /** - * The diagonal compensates for the mass (in 2D its just a length, in 3D its an - * area and we need the squared diagonal). * eps makes it such that the energy is relative to the envelope thickness. As it's a * squared energy, we need eps^2. */ - m_s_envelope = 1. / (m_params.diag_l * m_params.eps * m_params.eps); - /** - * The barrier energy computes a double-sided squared distance and therefore needs to be - * scaled by length^4. - */ - m_s_barrier = 1. / (m_params.diag_l4); - - // check weights (ignoring barrier here) - { - double& wa = m_params.w_amips; - double& ws = m_params.w_smooth; - const double sum = wa + ws; - if (sum > 1) { - wa /= sum; - ws /= sum; - logger().warn( - "Weights for AMIPS and smooth sum up to greater than 1. Rescaling to \n " - "w_amips = {}, \n w_smooth = {}", - wa, - ws); - } - double& we = m_params.w_envelope; - we = 1 - (wa + ws); - logger().info("w_envelope = {}", we); - } + m_s_envelope = 1. / (m_params.eps * m_params.eps); + - init_separation_weight(); + double& wa = m_params.w_amips; + double& we = m_params.w_envelope; + we = 1 - wa; + logger().info("w_envelope = {}", we); } ~ImageSimulationMeshTri() {} @@ -214,8 +178,6 @@ class ImageSimulationMeshTri : public wmtk::TriMesh void init_envelope(const MatrixXd& V, const MatrixXi& F); - void init_separation_weight(); - CellTag string_set_to_cell_tag(const std::set& str_set); bool adjust_sizing_field_serial(double max_energy); @@ -252,24 +214,13 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool smooth_before(const Tuple& t) override; bool smooth_after(const Tuple& t) override; - void build_mass_matrix(); /** * @brief A vector containing the vertex position and all positions of the surface neighbors. * * Returns an empty vector if vertex is not on the surface. */ std::vector get_surface_assembles(const Tuple& t) const; - std::shared_ptr get_smooth_energy(const Tuple& t) const; std::shared_ptr get_envelope_energy(const Tuple& t) const; - /** - * @brief Get the energy for minimal separation. - * - * By default, this method constructs the energy only using a local neighborhood. Using the - * entire surface is extremely slow. - */ - std::shared_ptr get_barrier_energy( - const Tuple& t, - const bool use_full_surface = false) const; std::vector> get_amips_assembles(const Tuple& t) const; std::shared_ptr get_amips_energy(const Tuple& t) const; @@ -365,18 +316,6 @@ class ImageSimulationMeshTri : public wmtk::TriMesh void tag_priority(const std::vector& tags_order); - /** - * @brief Find the substructure region that might be affected by minimal separation. - * - * The region considers all edges within 1.5 * (dhat + l_max), where l_max is the maximum - * incident edge length to the vertex and dhat is given by the parameters. - * - * @param v_tuple The vertex the region should be found for. - * @param V Nx2 vertex positions in the region - * @param E Nx2 edges in the region - */ - void substructure_region(const Tuple& v_tuple, MatrixXd& V, MatrixXi& E, size_t& vid) const; - bool vertex_is_on_surface(const size_t vid) const override { return m_vertex_attribute.at(vid).m_is_on_surface || diff --git a/components/image_simulation/wmtk/components/image_simulation/Parameters.h b/components/image_simulation/wmtk/components/image_simulation/Parameters.h index 310c4d3706..77afedf29b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/Parameters.h +++ b/components/image_simulation/wmtk/components/image_simulation/Parameters.h @@ -37,13 +37,8 @@ struct Parameters bool smooth_without_envelope = false; // weighting terms for the optimization - double w_amips = 1; - double w_smooth = 0; + double w_amips = 1e-4; double w_envelope = 0; - double w_separate = 0; - - double dhat = -1; - double separation_factor = 1; std::string operation = "remeshing"; @@ -76,8 +71,6 @@ struct Parameters } l_min = 0.5 * eps; - - dhat = separation_factor * l_min; } void init( const std::vector& vertices, diff --git a/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp b/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp index 14e064c56a..7038f7e263 100644 --- a/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/Smooth.cpp @@ -7,10 +7,8 @@ #include #include #include -#include #include #include -#include #include #include #include @@ -97,30 +95,12 @@ bool ImageSimulationMesh::smooth_after(const Tuple& t) // second solve (with proper weights) energy_sum = std::make_shared(); - if (m_params.w_smooth > 0) { - auto smooth_energy = get_smooth_energy(t); - if (smooth_energy) { - energy_sum->add_energy(smooth_energy); - } - } if (m_params.w_amips > 0) { energy_sum->add_energy(amips_energy); } if (m_params.w_envelope > 0) { energy_sum->add_energy(envelope_energy); } - // if (m_params.w_separate > 0) { - // auto barrier_energy = get_barrier_energy(t); - // if (barrier_energy) { - // double val = barrier_energy->value(old_pos); - // // only consider energy if it is non-zero at rest - // if (val > 0) { - // energy_sum->add_energy(barrier_energy); - // } else { - // logger().trace("Ignore barrier energy for zero rest state."); - // } - // } - // } total_energy = energy_sum; solve(); } else { @@ -172,8 +152,6 @@ bool ImageSimulationMesh::smooth_after(const Tuple& t) void ImageSimulationMesh::smooth_all_vertices(const size_t n_iters = 1) { - build_mass_matrix(); - for (size_t i = 0; i < n_iters; ++i) { igl::Timer timer; double time; @@ -238,33 +216,6 @@ void ImageSimulationMesh::smooth_all_vertices(const size_t n_iters = 1) } } -void ImageSimulationMesh::build_mass_matrix() -{ - MatrixXd V; - MatrixXi F; - - V.resize(vert_capacity(), 3); - for (size_t i = 0; i < V.rows(); ++i) { - V.row(i) = m_vertex_attribute.at(i).m_posf; - } - - const auto faces = get_faces_by_condition([](auto& f) { return f.m_is_surface_fs; }); - F.resize(faces.size(), 3); - for (size_t i = 0; i < F.rows(); ++i) { - F.row(i) = Vector3i(faces[i][0], faces[i][1], faces[i][2]); - } - - if (F.rows() == 0) { - logger().warn("No surface faces found. Mass matrix cannot be built."); - return; - } - optimization::BiharmonicEnergy3D::global_mass_and_stiffness( - V, - F, - m_surface_mass, - m_surface_stiffness); -} - std::vector> ImageSimulationMesh::get_amips_assembles(const Tuple& t) const { const size_t vid = t.vid(*this); @@ -303,47 +254,10 @@ std::shared_ptr ImageSimulationMesh::get_amips_en return amips_energy; } -std::shared_ptr ImageSimulationMesh::get_smooth_energy( - const Tuple& t) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid = t.vid(*this); - - std::vector adj; - optimization::BiharmonicEnergy3D::adjacency_from_stiffness(vid, m_surface_stiffness, adj); - - if (adj.empty()) { - return nullptr; - } - - MatrixXd pts; - pts.resize(adj.size() + 1, 3); - pts.row(0) = VA[vid].m_posf; - for (size_t j = 0; j < adj.size(); ++j) { - pts.row(j + 1) = VA[adj[j]].m_posf; - } - - double M; - VectorXd L_w; - optimization::BiharmonicEnergy3D::extract_local_mass_and_stiffness( - vid, - m_surface_mass, - m_surface_stiffness, - M, - L_w); - - const double w = m_s_smooth * m_params.w_smooth; - - auto smooth_energy = std::make_shared(pts, M, L_w, w); - - return smooth_energy; -} - std::shared_ptr ImageSimulationMesh::get_envelope_energy( const Tuple& t) const { size_t vid = t.vid(*this); - // const auto& M = m_surface_mass.coeff(vid, vid); const double w = m_s_envelope * m_params.w_envelope; auto env = m_envelope_orig; @@ -362,247 +276,4 @@ std::shared_ptr ImageSimulationMesh::get_envelope return envelope_energy; } -std::shared_ptr ImageSimulationMesh::get_barrier_energy( - const Tuple& t, - const bool use_full_surface) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid = t.vid(*this); - - if (!VA[vid].m_is_on_surface) { - return nullptr; - } - - // const auto& M = m_surface_mass.coeff(vid, vid); - // if (M <= 0) { - // return nullptr; - // } - - // barrier energy - MatrixXd V_barrier; - MatrixXi F_barrier; - size_t vid_barrier; - if (!use_full_surface) { - substructure_region(t, V_barrier, F_barrier, vid_barrier); - } else { - // this is horribly expensive and is only here for testing - - std::vector surf_points; - std::vector global_to_local_vid_map(vert_capacity()); - for (size_t i = 0; i < vert_capacity(); ++i) { - const Tuple v = tuple_from_vertex(i); - if (!v.is_valid(*this)) { - continue; - } - const size_t _vid = v.vid(*this); - if (!m_vertex_attribute.at(_vid).m_is_on_surface) { - continue; - } - surf_points.push_back(m_vertex_attribute.at(_vid).m_posf); - global_to_local_vid_map[_vid] = surf_points.size() - 1; - } - - V_barrier.resize(surf_points.size(), 3); - for (size_t i = 0; i < surf_points.size(); ++i) { - V_barrier.row(i) = surf_points[i]; - } - - const auto faces = get_faces_by_condition([](auto& f) { return f.m_is_surface_fs; }); - F_barrier.resize(faces.size(), 3); - for (size_t i = 0; i < faces.size(); ++i) { - const size_t v0 = global_to_local_vid_map[faces[i][0]]; - const size_t v1 = global_to_local_vid_map[faces[i][1]]; - const size_t v2 = global_to_local_vid_map[faces[i][2]]; - F_barrier.row(i) = Vector3i(v0, v1, v2); - } - - vid_barrier = global_to_local_vid_map[vid]; - } - - if (F_barrier.rows() == 0) { - return nullptr; - } - - MatrixXi E_barrier; - igl::edges(F_barrier, E_barrier); - - auto barrier_energy = std::make_shared( - V_barrier, - E_barrier, - F_barrier, - vid_barrier, - m_params.dhat, - m_s_barrier * m_params.w_separate); - - return barrier_energy; -} - -void ImageSimulationMesh::substructure_region( - const Tuple& v_tuple, - MatrixXd& V, - MatrixXi& F, - size_t& vid_local) const -{ - const auto& VA = m_vertex_attribute; - const size_t vid_global = v_tuple.vid(*this); - - if (!VA[vid_global].m_is_on_surface) { - log_and_throw_error( - "Cannot compute substructure region for vertex that is not on the surface"); - // We actually could compute this but it doesn't make sense. Remove this if-statement if we - // should ever need to compute regions for vertices not on the surface. - } - - const Vector3d& p0 = VA[vid_global].m_posf; - - double l_max = 0; // longest incident edge length - for (const size_t& v1 : get_one_ring_vids_for_vertex(vid_global)) { - const Vector3d& p1 = VA[v1].m_posf; - const double l_squared = (p1 - p0).squaredNorm(); - l_max = std::max(l_max, l_squared); - } - l_max = std::sqrt(l_max); - - const double r = 1.5 * m_params.dhat + l_max; // region radius - const double r2 = r * r; - - simplex::SimplexCollection candidates; // all faces within the region - - // make a BFS to find faces within distance - std::unordered_set visited; - std::queue q; - q.push(vid_global); - visited.insert(vid_global); - - while (!q.empty()) { - const size_t v_a = q.front(); - q.pop(); - - const Vector3d& p_a = VA[v_a].m_posf; - const double d2_a = (p0 - p_a).squaredNorm(); - - for (const size_t& tid : get_one_ring_tids_for_vertex(vid_global)) { - // get all faces from tet - simplex::Tet tet = simplex_from_tet(tid); - const auto fs = tet.faces(); - - for (size_t i = 0; i < 4; ++i) { - const Vector3d& t1 = VA[fs[i].vertices()[0]].m_posf; - const Vector3d& t2 = VA[fs[i].vertices()[1]].m_posf; - const Vector3d& t3 = VA[fs[i].vertices()[2]].m_posf; - - // if p_a is within the region, skip distance test for triangles - if (d2_a <= r2) { - const double d2 = ipc::point_triangle_distance(p0, t1, t2, t3); - if (d2 > r2) { - // triangle is too far away - continue; - } - } - - candidates.add(fs[i]); - // add vertices to queue - for (size_t j = 0; j < 3; ++j) { - const size_t v_b = fs[i].vertices()[j]; - if (visited.count(v_b) == 0) { - q.push(v_b); - visited.insert(v_b); - } - } - } - } - } - candidates.sort_and_clean(); - - std::vector faces; // faces on the surface - for (const auto& f : candidates.faces()) { - const auto [tup, fid] = tuple_from_face(f.vertices()); - if (face_is_on_surface(fid)) { - faces.push_back(f); - } - } - - std::unordered_set region_vertices; - for (const auto& f : faces) { - region_vertices.insert(f.vertices()[0]); - region_vertices.insert(f.vertices()[1]); - region_vertices.insert(f.vertices()[2]); - } - - // build V and F - std::vector points; - std::unordered_map global_to_local_vid_map; - for (const size_t i : region_vertices) { - const Tuple v = tuple_from_vertex(i); - assert(v.is_valid(*this)); - assert(VA[i].m_is_on_surface); - points.push_back(VA[i].m_posf); - global_to_local_vid_map[i] = points.size() - 1; - } - - V.resize(points.size(), 3); - for (size_t i = 0; i < points.size(); ++i) { - V.row(i) = points[i]; - } - - F.resize(faces.size(), 3); - for (size_t i = 0; i < faces.size(); ++i) { - const size_t v0 = global_to_local_vid_map[faces[i].vertices()[0]]; - const size_t v1 = global_to_local_vid_map[faces[i].vertices()[1]]; - const size_t v2 = global_to_local_vid_map[faces[i].vertices()[2]]; - F.row(i) = Vector3i(v0, v1, v2); - } - - vid_local = global_to_local_vid_map[vid_global]; -} - -void ImageSimulationMesh::init_separation_weight() -{ - if (m_params.preserve_topology) { - // no separation for not preserving topology - m_params.w_separate = 0; - } - if (m_params.separation_factor <= 0) { - m_params.w_separate = 0; - logger().warn( - "Separation factor is {} -> No separation. IPC barrier weight = {}", - m_params.separation_factor, - m_params.w_separate); - return; - } - - /** - * This assumes that the barrier term is the clamped log barrier term: - * E_B = b(d) = -(d-\hat{d})^2\ln\left(\frac{d}{\hat{d}}\right) - * - * E_B' = b'(\hat{d}/2) = \hat{d} * (\ln(2) - 0.5) - * - * barrier scaling: s_B = 1/l^4 - * envelope scaling: s_E = 1 / (l * eps^2) - * - * l is the bounding box diagonal. - * - * The envelope energy E_E = \Vert p - p_0 \Vert^2 - * E_E'(\hat{d}/2) = \hat{d} - * - * We compute w_E from - * s_B w_B E_B' = s_E w_E E_E' - * - * w_B = w_E l^3 / (eps^2 * (\ln(2) - 0.5)) - */ - - const double l3 = m_params.diag_l3; - const double eps2 = m_params.eps * m_params.eps; - const double ln2 = std::log(2.0); - - m_params.w_separate = m_params.w_envelope * l3 / (eps2 * (ln2 - 0.5)); - - logger().info("IPC barrier weight = {}, dhat = {}", m_params.w_separate, m_params.dhat); - - if (m_params.w_separate < 0) { - logger().error("Negative separation weight: {}. Setting weight to 0.", m_params.w_separate); - m_params.w_separate = 0; - } -} - } // namespace wmtk::components::image_simulation \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index b656892dc8..ebae52a487 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -50,8 +50,6 @@ void run_3D(const nlohmann::json& json_params, const InputData& input_data) params.smooth_without_envelope = json_params["smooth_without_envelope"]; params.w_amips = json_params["w_amips"]; - params.w_smooth = json_params["w_smooth"]; - params.separation_factor = json_params["separation_factor"]; const bool write_vtu = json_params["write_vtu"]; @@ -228,8 +226,6 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) params.smooth_without_envelope = json_params["smooth_without_envelope"]; params.w_amips = json_params["w_amips"]; - params.w_smooth = json_params["w_smooth"]; - params.separation_factor = json_params["separation_factor"]; const bool write_vtu = json_params["write_vtu"]; diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp index 7ae7150075..9b2f0cffd3 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp @@ -28,8 +28,6 @@ nlohmann::json image_simulation_spec = R"( "stop_at_float", "preserve_topology", "w_amips", - "w_smooth", - "separation_factor", "smooth_without_envelope", "write_vtu", "write_envelope", @@ -220,19 +218,9 @@ nlohmann::json image_simulation_spec = R"( "pointer": "/w_amips", "type": "float", "default": 1e-4, - "doc": "AMIPS energy. Sum of energy weights must be smaller than 1!" - }, - { - "pointer": "/w_smooth", - "type": "float", - "default": 0, - "doc": "Surface smoothing energy. Sum of energy weights must be smaller than 1!" - }, - { - "pointer": "/separation_factor", - "type": "float", - "default": 1, - "doc": "Usually, surfaces are pushed apart to create a gap in the size of the smallest allowed edge length. The separation factor can increase or decrease that gap. If set to 0 or negative, no separation is performed." + "min": 0, + "max": 1, + "doc": "AMIPS energy" }, { "pointer": "/smooth_without_envelope", diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index 5c19e82b82..c51f58405b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -23,8 +23,6 @@ "stop_at_float", "preserve_topology", "w_amips", - "w_smooth", - "separation_factor", "smooth_without_envelope", "write_vtu", "write_envelope", @@ -215,19 +213,9 @@ "pointer": "/w_amips", "type": "float", "default": 1e-4, - "doc": "AMIPS energy. Sum of energy weights must be smaller than 1!" - }, - { - "pointer": "/w_smooth", - "type": "float", - "default": 0, - "doc": "Surface smoothing energy. Sum of energy weights must be smaller than 1!" - }, - { - "pointer": "/separation_factor", - "type": "float", - "default": 1, - "doc": "Usually, surfaces are pushed apart to create a gap in the size of the smallest allowed edge length. The separation factor can increase or decrease that gap. If set to 0 or negative, no separation is performed." + "min": 0, + "max": 1, + "doc": "AMIPS energy" }, { "pointer": "/smooth_without_envelope", diff --git a/components/image_simulation/wmtk/components/image_simulation/tests/test_image_simulation_energies.cpp b/components/image_simulation/wmtk/components/image_simulation/tests/test_image_simulation_energies.cpp index 435106599e..764f59ace6 100644 --- a/components/image_simulation/wmtk/components/image_simulation/tests/test_image_simulation_energies.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/tests/test_image_simulation_energies.cpp @@ -58,18 +58,14 @@ TEST_CASE("scale-invariance", test_groups + test_release_only) struct Energies { - double amips = -1, smooth = -1, envelope = -1, barrier = -1; + double amips = -1, envelope = -1; }; auto collect_energies = [](double scale) { Parameters params; - params.w_amips = 1e-8; - params.w_smooth = 1e-4; auto mesh_ptr = init_optimization_tests(params, scale); tri::ImageSimulationMeshTri& mesh = *mesh_ptr; - mesh.build_mass_matrix(); - const auto& VA = mesh.m_vertex_attribute; std::vector energies; @@ -83,15 +79,8 @@ TEST_CASE("scale-invariance", test_groups + test_release_only) continue; } - auto smooth_energy = mesh.get_smooth_energy(t); - - if (!smooth_energy) { - continue; - } - auto amips_energy = mesh.get_amips_energy(t); auto envelope_energy = mesh.get_envelope_energy(t); - auto barrier_energy = mesh.get_barrier_energy(t); // find a position to test the envelope energy Tuple n; @@ -108,9 +97,7 @@ TEST_CASE("scale-invariance", test_groups + test_release_only) Energies e; e.amips = amips_energy->value(x); - e.smooth = smooth_energy->value(y); e.envelope = envelope_energy->value(y); - e.barrier = barrier_energy->value(x); energies.push_back(e); } @@ -126,64 +113,7 @@ TEST_CASE("scale-invariance", test_groups + test_release_only) // compare all energies // logger().info("amips"); compare(a[i].amips, b[i].amips); - // logger().info("smooth"); - compare(a[i].smooth, b[i].smooth); // logger().info("env"); compare(a[i].envelope, b[i].envelope); - // logger().info("barrier"); - compare(a[i].barrier, b[i].barrier); - } -} - -TEST_CASE("barrier-bfs", test_groups + test_release_only) -{ - logger().set_level(spdlog::level::off); - - using Tuple = TriMesh::Tuple; - Parameters params; - auto mesh_ptr = init_optimization_tests(params); - tri::ImageSimulationMeshTri& mesh = *mesh_ptr; - - mesh.build_mass_matrix(); - - const auto& VA = mesh.m_vertex_attribute; - - for (const Tuple& t : mesh.get_vertices()) { - const size_t vid = t.vid(mesh); - if (vid < 8000) { - continue; - } - if (!VA[vid].m_is_on_surface) { - continue; - } - - mesh.m_s_barrier = 1; // remove weight - - auto loc = mesh.get_barrier_energy(t, false); - auto glob = mesh.get_barrier_energy(t, true); - - if (!loc) { - continue; - } - - const Vector2d& x = VA[vid].m_pos; - - double v_loc = loc->value(x); - VectorXd g_loc; - loc->gradient(x, g_loc); - MatrixXd h_loc; - loc->hessian(x, h_loc); - - double v_glob = glob->value(x); - VectorXd g_glob; - glob->gradient(x, g_glob); - MatrixXd h_glob; - glob->hessian(x, h_glob); - - // logger() - // .info("v = {}: g_loc = {}, g_glob = {}", vid, g_loc.transpose(), g_glob.transpose()); - - CHECK((g_loc - g_glob).squaredNorm() < 1e-10); - CHECK((h_loc - h_glob).squaredNorm() < 1e-10); } } \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 7292b0c2a4..f311252077 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,8 +27,8 @@ build-dir = "build/RelWithDebInfo" cmake.args = [ "-DWMTK_PYBIND=ON", "-DCMAKE_INSTALL_DEFAULT_COMPONENT_NAME=python_bindings", - "-DBUILD_TESTING=OFF", - "-DWMTK_APP_HARMTET=OFF", - "-DWMTK_APP_INTERIOR_TET_OPT=OFF", + # "-DBUILD_TESTING=OFF", + # "-DWMTK_APP_HARMTET=OFF", + # "-DWMTK_APP_INTERIOR_TET_OPT=OFF", ] install.components = ["python_bindings"] diff --git a/src/wmtk/optimization/BarrierEnergy.cpp b/src/wmtk/optimization/BarrierEnergy.cpp deleted file mode 100644 index f8e4852e4b..0000000000 --- a/src/wmtk/optimization/BarrierEnergy.cpp +++ /dev/null @@ -1,191 +0,0 @@ -#include "BarrierEnergy.hpp" - - -namespace wmtk::optimization { - -BarrierEnergy2D::BarrierEnergy2D( - const MatrixXd& V, - const MatrixXi& E, - const size_t vid, - const double dhat, - const double weight) - : m_collision_mesh(V, E) - , m_V(V) -#ifndef WMTK_SMOOTH_BARRIER - , m_B(dhat) -#endif - , m_x0(V.row(vid)) - , m_vid(vid) - , m_weight(weight) -#ifdef WMTK_SMOOTH_BARRIER - , m_params(dhat, 1.0, 0.0, 1.0, 0.0, 1) - , m_smooth_B(m_params) -#endif -{ - assert(V.cols() == 2); - assert(E.cols() == 2); - -#ifdef WMTK_SMOOTH_BARRIER - m_smooth_collisions.build(m_collision_mesh, m_V, m_params); -#else - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -#endif -} - -BarrierEnergy2D::TVector BarrierEnergy2D::initial_position() const -{ - return m_x0; -} - -void BarrierEnergy2D::replace_vid(const size_t vid) -{ - //// reset m_V - // m_V.row(m_vid) = m_x0; - - // set new vid - m_vid = vid; - m_x0 = m_V.row(m_vid); - -#ifdef WMTK_SMOOTH_BARRIER - m_smooth_collisions.build(m_collision_mesh, m_V, m_params); -#else - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -#endif -} - - -double BarrierEnergy2D::value(const TVector& x) -{ - assert(x.size() == 2); - update_collisions(x); -#ifdef WMTK_SMOOTH_BARRIER - double potential = m_smooth_B(m_smooth_collisions, m_collision_mesh, m_V); -#else - double potential = m_B(m_collisions, m_collision_mesh, m_V); -#endif - return m_weight * potential; -} - -void BarrierEnergy2D::gradient(const TVector& x, TVector& gradv) -{ - assert(x.size() == 2); - update_collisions(x); -#ifdef WMTK_SMOOTH_BARRIER - VectorXd grad = m_smooth_B.gradient(m_smooth_collisions, m_collision_mesh, m_V); -#else - VectorXd grad = m_B.gradient(m_collisions, m_collision_mesh, m_V); -#endif - gradv = m_weight * grad.segment<2>(m_vid * 2); -} - -void BarrierEnergy2D::hessian(const TVector& x, MatrixXd& hessian) -{ - assert(x.size() == 2); - update_collisions(x); -#ifdef WMTK_SMOOTH_BARRIER - Eigen::SparseMatrix hess = - m_smooth_B.hessian(m_smooth_collisions, m_collision_mesh, m_V); -#else - Eigen::SparseMatrix hess = m_B.hessian(m_collisions, m_collision_mesh, m_V); -#endif - - Matrix2d v_hess; - v_hess.setZero(); - v_hess(0, 0) = hess.coeff(2 * m_vid, 2 * m_vid); - v_hess(1, 1) = hess.coeff(2 * m_vid + 1, 2 * m_vid + 1); - hessian = m_weight * v_hess; -} - -void BarrierEnergy2D::update_collisions(const TVector& x) -{ - assert(m_V.cols() == x.size()); - const VectorXd& curr = m_V.row(m_vid); - if (curr == x) { - return; - } - - m_V.row(m_vid) = x; -#ifdef WMTK_SMOOTH_BARRIER - m_smooth_collisions.build(m_collision_mesh, m_V, m_params); -#else - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -#endif -} - - -BarrierEnergy3D::BarrierEnergy3D( - const MatrixXd& V, - const MatrixXi& E, - const MatrixXi& F, - const size_t vid, - const double dhat, - const double weight) - : m_collision_mesh(V, E, F) - , m_V(V) - , m_B(dhat) - , m_x0(V.row(vid)) - , m_vid(vid) - , m_weight(weight) -{ - assert(V.cols() == 3); - assert(F.cols() == 3); - - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -} - -BarrierEnergy3D::TVector BarrierEnergy3D::initial_position() const -{ - return m_x0; -} - -void BarrierEnergy3D::replace_vid(const size_t vid) -{ - // set new vid - m_vid = vid; - m_x0 = m_V.row(m_vid); - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -} - -double BarrierEnergy3D::value(const TVector& x) -{ - assert(x.size() == 3); - update_collisions(x); - double potential = m_B(m_collisions, m_collision_mesh, m_V); - return m_weight * potential; -} - -void BarrierEnergy3D::gradient(const TVector& x, TVector& gradv) -{ - assert(x.size() == 3); - update_collisions(x); - VectorXd grad = m_B.gradient(m_collisions, m_collision_mesh, m_V); - gradv = m_weight * grad.segment<3>(m_vid * 3); -} - -void BarrierEnergy3D::hessian(const TVector& x, MatrixXd& hessian) -{ - assert(x.size() == 3); - update_collisions(x); - Eigen::SparseMatrix hess = m_B.hessian(m_collisions, m_collision_mesh, m_V); - - Matrix3d v_hess; - v_hess.setZero(); - v_hess(0, 0) = hess.coeff(3 * m_vid + 0, 3 * m_vid + 0); - v_hess(1, 1) = hess.coeff(3 * m_vid + 1, 3 * m_vid + 1); - v_hess(2, 2) = hess.coeff(3 * m_vid + 2, 3 * m_vid + 2); - hessian = m_weight * v_hess; -} - -void BarrierEnergy3D::update_collisions(const TVector& x) -{ - assert(m_V.cols() == x.size()); - const VectorXd& curr = m_V.row(m_vid); - if (curr == x) { - return; - } - - m_V.row(m_vid) = x; - m_collisions.build(m_collision_mesh, m_V, m_B.dhat()); -} - -} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/BarrierEnergy.hpp b/src/wmtk/optimization/BarrierEnergy.hpp deleted file mode 100644 index dab31fae33..0000000000 --- a/src/wmtk/optimization/BarrierEnergy.hpp +++ /dev/null @@ -1,131 +0,0 @@ -#pragma once - -// #define WMTK_SMOOTH_BARRIER - -#include -#include -#include -#ifdef WMTK_SMOOTH_BARRIER -#include -#endif -#include -#include - - -namespace wmtk::optimization { - -class BarrierEnergy2D : public polysolve::nonlinear::Problem -{ -public: - using typename polysolve::nonlinear::Problem::Scalar; - using typename polysolve::nonlinear::Problem::THessian; - using typename polysolve::nonlinear::Problem::TVector; - - /** - * @brief Barrier energy for a polyline in 2D - * - * The energy is defined over an entire polyline described by V and E. But the optimization only - * considers one vertex with ID `vid`. The `vid` can be replaced to optimize another vertex - * position later on. - */ - BarrierEnergy2D( - const MatrixXd& V, - const MatrixXi& E, - const size_t vid, - const double dhat, - const double weight = 1); - - TVector initial_position() const; - - void replace_vid(const size_t vid); - - MatrixXd& V() { return m_V; } - - double value(const TVector& x) override; - void gradient(const TVector& x, TVector& gradv) override; - void hessian(const TVector& x, THessian& hessian) override - { - log_and_throw_error("Sparse functions do not exist, use dense solver"); - } - void hessian(const TVector& x, MatrixXd& hessian) override; - - void solution_changed(const TVector& new_x) override {} - - bool is_step_valid(const TVector& x0, const TVector& x1) override { return true; } - - void update_collisions(const TVector& x); - -private: - ipc::CollisionMesh m_collision_mesh; - MatrixXd m_V; -#ifdef WMTK_SMOOTH_BARRIER - ipc::SmoothContactParameters m_params; - ipc::SmoothCollisions m_smooth_collisions; - ipc::SmoothContactPotential m_smooth_B; -#else - ipc::NormalCollisions m_collisions; - ipc::BarrierPotential m_B; -#endif - - - Vector2d m_x0; - size_t m_vid; - - double m_weight; -}; - -class BarrierEnergy3D : public polysolve::nonlinear::Problem -{ -public: - using typename polysolve::nonlinear::Problem::Scalar; - using typename polysolve::nonlinear::Problem::THessian; - using typename polysolve::nonlinear::Problem::TVector; - - /** - * @brief Barrier energy for a polyline in 2D - * - * The energy is defined over an entire polyline described by V and E. But the optimization only - * considers one vertex with ID `vid`. The `vid` can be replaced to optimize another vertex - * position later on. - */ - BarrierEnergy3D( - const MatrixXd& V, - const MatrixXi& E, - const MatrixXi& F, - const size_t vid, - const double dhat, - const double weight = 1); - - TVector initial_position() const; - - void replace_vid(const size_t vid); - - MatrixXd& V() { return m_V; } - - double value(const TVector& x) override; - void gradient(const TVector& x, TVector& gradv) override; - void hessian(const TVector& x, THessian& hessian) override - { - log_and_throw_error("Sparse functions do not exist, use dense solver"); - } - void hessian(const TVector& x, MatrixXd& hessian) override; - - void solution_changed(const TVector& new_x) override {} - - bool is_step_valid(const TVector& x0, const TVector& x1) override { return true; } - - void update_collisions(const TVector& x); - -private: - ipc::CollisionMesh m_collision_mesh; - MatrixXd m_V; - ipc::NormalCollisions m_collisions; - ipc::BarrierPotential m_B; - - Vector3d m_x0; - size_t m_vid; - - double m_weight; -}; - -} // namespace wmtk::optimization \ No newline at end of file diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp index c867cd0e28..ebf43c40e3 100644 --- a/tests/test_optimization.cpp +++ b/tests/test_optimization.cpp @@ -3,13 +3,9 @@ #include #include #include -#include -#include -#include #include #include #include -#include #include #include #include @@ -390,202 +386,4 @@ TEST_CASE("biharmonic_energy_bunny_3d", "[energies][.]") // igl::writeOFF(fmt::format("debug_{}.off", i), V, F); } -} - -TEST_CASE("ipc_2d", "[energies][ipc]") -{ - // two edges - MatrixXd V; - V.resize(4, 2); - // e0 - V.row(0) = Vector2d(0, 0); - V.row(1) = Vector2d(10, 0); - // e1 - V.row(2) = Vector2d(0, 1e-3); - V.row(3) = Vector2d(10, 1e-3); - - MatrixXi E; - E.resize(1, 2); - E.row(0) = Vector2i(0, 1); - // E.row(1) = Vector2i(2, 3); - - ipc::CollisionMesh collision_mesh(V, E); - - const double dhat = 1; - - ipc::Candidates candidates; // gather all interesting pairs in here - candidates.ev_candidates.emplace_back(0, 2); - candidates.ev_candidates.emplace_back(0, 3); - - ipc::NormalCollisions collisions; - // collisions.build(collision_mesh, V, dhat); - collisions.build(candidates, collision_mesh, V, dhat); - - const ipc::BarrierPotential B(dhat); - double barrier_potential = B(collisions, collision_mesh, V); - // logger().info("Potential = {}", barrier_potential); - - VectorXd barrier_potential_grad = B.gradient(collisions, collision_mesh, V); - // logger().info("Grad Potential = \n{}", barrier_potential_grad); - - size_t i = 2; - Vector2d v2_grad(barrier_potential_grad[2 * i + 0], barrier_potential_grad[2 * i + 1]); - // logger().info("Grad Potential v2 = {}", v2_grad.transpose()); - - CHECK(v2_grad[1] < 0); - CHECK(v2_grad[0] == 0); - - Eigen::SparseMatrix barrier_potential_hess = B.hessian(collisions, collision_mesh, V); - Matrix2d v2_hess; - v2_hess.setZero(); - v2_hess(0, 0) = barrier_potential_hess.coeff(2 * i, 2 * i); - v2_hess(1, 1) = barrier_potential_hess.coeff(2 * i + 1, 2 * i + 1); - // logger().info("Hess v2 = \n{}", v2_hess); - - CHECK(v2_hess(0, 0) == 0); - CHECK(v2_hess(0, 1) == 0); - CHECK(v2_hess(1, 0) == 0); - CHECK(v2_hess(1, 1) > 0); - // logger().info("Hess full = \n{}", MatrixXd(barrier_potential_hess)); -} - -TEST_CASE("barrier_energy_2d", "[energies][ipc]") -{ - // 4 edges, two top two bottom - MatrixXd V; - V.resize(6, 2); - // bottom two edges - V.row(0) = Vector2d(0, 0); - V.row(1) = Vector2d(10, 0); - V.row(2) = Vector2d(20, 0); - // top two edges - V.row(3) = Vector2d(0, 1e-3); - V.row(4) = Vector2d(10, 1e-3); - V.row(5) = Vector2d(20, 1e-3); - - MatrixXi E; - E.resize(4, 2); - E.row(0) = Vector2i(0, 1); - E.row(1) = Vector2i(1, 2); - E.row(2) = Vector2i(3, 4); - E.row(3) = Vector2i(4, 5); - - - optimization::BarrierEnergy2D energy(V, E, 4, 1); - const Vector2d p0 = energy.initial_position(); - CHECK(energy.is_step_valid(p0, p0)); - - { - VectorXd g; - energy.gradient(p0, g); - REQUIRE(g.size() == 2); - CHECK(g[1] < 0); - - MatrixXd h; - CHECK_NOTHROW(energy.hessian(p0, h)); - } - auto x = energy.initial_position(); - const double e_before = energy.value(x); - CHECK(e_before > 0); - { - auto linear_solver_params = optimization::basic_linear_solver_params; - auto nonlinear_solver_params = optimization::basic_nonlinear_solver_params; - nonlinear_solver_params["max_iterations"] = 100; - - auto solver = polysolve::nonlinear::Solver::create( - nonlinear_solver_params, - linear_solver_params, - 1, - opt_logger()); - optimization::deactivate_opt_logger(); - - CHECK_NOTHROW(solver->minimize(energy, x)); - } - const double e_after = energy.value(x); - CHECK(e_after < e_before); - { - VectorXd g; - energy.gradient(x, g); - CHECK(g[1] < 1e-6); - } -} - -TEST_CASE("barrier_plus_biharmonic_energy_2d", "[energies][ipc]") -{ - // 4 edges, two top two bottom - MatrixXd V; - V.resize(6, 2); - // bottom two edges - V.row(0) = Vector2d(0, 0); - V.row(1) = Vector2d(10, 0); - V.row(2) = Vector2d(20, 0); - // top two edges - V.row(3) = Vector2d(0, 1e-3); - V.row(4) = Vector2d(10, 1e-3); - V.row(5) = Vector2d(20, 1e-3); - std::array top_vs; - top_vs[0] = V.row(4); // optimized vertex - top_vs[1] = V.row(3); // neighbor - top_vs[2] = V.row(5); // neighbor - - MatrixXi E; - E.resize(4, 2); - E.row(0) = Vector2i(0, 1); - E.row(1) = Vector2i(1, 2); - E.row(2) = Vector2i(3, 4); - E.row(3) = Vector2i(4, 5); - - optimization::EnergySum energy; - - std::shared_ptr barrier_energy = - std::make_shared(V, E, 4, 1); - - double M; - Vector3d L_w; - optimization::BiharmonicEnergy2D::local_mass_and_stiffness(top_vs, M, L_w); - std::shared_ptr biharmonic_energy = - std::make_shared(top_vs, M, L_w); - - energy.add_energy(barrier_energy, 1); - energy.add_energy(biharmonic_energy, 1e3); - - const Vector2d p0 = barrier_energy->initial_position(); - CHECK(barrier_energy->is_step_valid(p0, p0)); - - - { - VectorXd g; - energy.gradient(p0, g); - REQUIRE(g.size() == 2); - CHECK(g[1] < 0); - - MatrixXd h; - CHECK_NOTHROW(energy.hessian(p0, h)); - } - auto x = barrier_energy->initial_position(); - const double e_before = energy.value(x); - CHECK(e_before > 0); - { - auto linear_solver_params = optimization::basic_linear_solver_params; - auto nonlinear_solver_params = optimization::basic_nonlinear_solver_params; - nonlinear_solver_params["max_iterations"] = 100; - - auto solver = polysolve::nonlinear::Solver::create( - nonlinear_solver_params, - linear_solver_params, - 1, - opt_logger()); - optimization::deactivate_opt_logger(); - - CHECK_NOTHROW(solver->minimize(energy, x)); - } - const double e_after = energy.value(x); - logger().info("x = {}", x); - CHECK(e_after < e_before); - - // vertex should be somewhere between 1e-3 and 1 - CHECK(x[1] > 1e-3); - CHECK(x[1] < 1); - // vertex should not move to the side - CHECK(std::abs(x[0] - 10) < 1e-4); -} +} \ No newline at end of file From e6c01a697114396adc7314b9a8b27b9f59184efa Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 20 May 2026 17:08:07 +0200 Subject: [PATCH 2/3] Fix for the substructure link condition in 3D. --- src/wmtk/TetMeshSubstructure.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/wmtk/TetMeshSubstructure.cpp b/src/wmtk/TetMeshSubstructure.cpp index 2f3468ceb3..6755644d38 100644 --- a/src/wmtk/TetMeshSubstructure.cpp +++ b/src/wmtk/TetMeshSubstructure.cpp @@ -239,8 +239,10 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const SimplexCollection link_u_0; SimplexCollection link_u_1; + SimplexCollection link_u_2; SimplexCollection link_v_0; SimplexCollection link_v_1; + SimplexCollection link_v_2; SimplexCollection link_e_0; SimplexCollection link_e_1; @@ -300,9 +302,12 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const Edge ew = fw.opposite_edge(u); link_u_0.add_with_faces(ew); link_u_1.add_with_faces(ew); + const Vertex vw = ew.opposite_vertex(u); + link_u_2.add(vw); } link_u_0.sort_and_clean(); link_u_1.sort_and_clean(); + link_u_2.sort_and_clean(); } // vertex v links { @@ -352,10 +357,19 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const Edge ew = fw.opposite_edge(v); link_v_0.add_with_faces(ew); link_v_1.add_with_faces(ew); + const Vertex vw = ew.opposite_vertex(v); + link_v_2.add(vw); } link_v_0.sort_and_clean(); link_v_1.sort_and_clean(); + link_v_2.sort_and_clean(); } + + const auto link_uv_2 = SimplexCollection::get_intersection(link_u_2, link_v_2); + if (!link_uv_2.empty()) { + return false; + } + // edge links { const Edge e(u_id, v_id); From 88630502cc49f6505765e0c9366ab9e85a7f0a9d Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 20 May 2026 17:59:21 +0200 Subject: [PATCH 3/3] Fix bug in link condition. --- src/wmtk/TetMeshSubstructure.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/wmtk/TetMeshSubstructure.cpp b/src/wmtk/TetMeshSubstructure.cpp index 6755644d38..394cce499a 100644 --- a/src/wmtk/TetMeshSubstructure.cpp +++ b/src/wmtk/TetMeshSubstructure.cpp @@ -302,7 +302,7 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const Edge ew = fw.opposite_edge(u); link_u_0.add_with_faces(ew); link_u_1.add_with_faces(ew); - const Vertex vw = ew.opposite_vertex(u); + const Vertex vw = e.opposite_vertex(u); link_u_2.add(vw); } link_u_0.sort_and_clean(); @@ -357,7 +357,7 @@ bool TetMesh::substructure_link_condition(const Tuple& e_tuple) const const Edge ew = fw.opposite_edge(v); link_v_0.add_with_faces(ew); link_v_1.add_with_faces(ew); - const Vertex vw = ew.opposite_vertex(v); + const Vertex vw = e.opposite_vertex(v); link_v_2.add(vw); } link_v_0.sort_and_clean();