From 175cd00ec3b8130d2300ec1d8adc21c282c30dea Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Tue, 10 Mar 2026 12:41:35 +0100 Subject: [PATCH 01/19] Add polysolve and some 2D energies. In the future, this should replace the current Newton solver in the toolkit. --- CMakeLists.txt | 2 + cmake/recipes/polysolve.cmake | 33 ++++ .../ImageSimulationMeshTri.cpp | 68 +++++--- .../ImageSimulationMeshTri.hpp | 6 + src/wmtk/Types.hpp | 2 + src/wmtk/optimization/AMIPSEnergy.cpp | 76 +++++++++ src/wmtk/optimization/AMIPSEnergy.hpp | 42 +++++ src/wmtk/optimization/DirichletEnergy.cpp | 48 ++++++ src/wmtk/optimization/DirichletEnergy.hpp | 40 +++++ src/wmtk/optimization/EnergySum.cpp | 59 +++++++ src/wmtk/optimization/EnergySum.hpp | 44 +++++ src/wmtk/optimization/solver.hpp | 26 +++ src/wmtk/utils/Rational.hpp | 2 +- src/wmtk/utils/orient.cpp | 75 +++++++++ src/wmtk/utils/orient.hpp | 15 ++ tests/CMakeLists.txt | 2 + tests/test_optimization.cpp | 153 ++++++++++++++++++ tests/test_orient.cpp | 92 +++++++++++ 18 files changed, 766 insertions(+), 19 deletions(-) create mode 100644 cmake/recipes/polysolve.cmake create mode 100644 src/wmtk/optimization/AMIPSEnergy.cpp create mode 100644 src/wmtk/optimization/AMIPSEnergy.hpp create mode 100644 src/wmtk/optimization/DirichletEnergy.cpp create mode 100644 src/wmtk/optimization/DirichletEnergy.hpp create mode 100644 src/wmtk/optimization/EnergySum.cpp create mode 100644 src/wmtk/optimization/EnergySum.hpp create mode 100644 src/wmtk/optimization/solver.hpp create mode 100644 src/wmtk/utils/orient.cpp create mode 100644 src/wmtk/utils/orient.hpp create mode 100644 tests/test_optimization.cpp create mode 100644 tests/test_orient.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b64f2fed9b..12d7418cf3 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -70,6 +70,7 @@ include(spdlog) include(paraviewo) include(fenvelope) include(nanoflann) +include(polysolve) include(simple_bvh) # Sort projects inside the solution @@ -128,6 +129,7 @@ target_link_libraries(wildmeshing_toolkit PUBLIC paraviewo::paraviewo FastEnvelope::FastEnvelope simple_bvh::simple_bvh + polysolve::polysolve ) if(MSVC) diff --git a/cmake/recipes/polysolve.cmake b/cmake/recipes/polysolve.cmake new file mode 100644 index 0000000000..e54b8c781f --- /dev/null +++ b/cmake/recipes/polysolve.cmake @@ -0,0 +1,33 @@ +# PolySolve (https://github.com/polyfem/polysolve) +# License: MIT + +if(TARGET polysolve::polysolve) + return() +endif() + +message(STATUS "Third-party: creating target 'polysolve::polysolve'") + +# # TODO: this requires a conflicting version of Eigen. Reenable when Eigen 3.4+ is available. +# set(POLYSOLVE_WITH_ACCELERATE OFF CACHE BOOL "Enable Apple Accelerate" FORCE) + +include(CPM) +CPMAddPackage( + NAME polysolve + GITHUB_REPOSITORY polyfem/polysolve + GIT_TAG 57dbeadcb65932a38508e3a5a07e98732a1c50bc + OPTIONS + "POLYSOLVE_WITH_ACCELERATE OFF" + "POLYSOLVE_WITH_CHOLMOD OFF" + "POLYSOLVE_WITH_UMFPACK OFF" + "POLYSOLVE_WITH_SUPERLU OFF" + "POLYSOLVE_WITH_SPQR OFF" + "POLYSOLVE_WITH_MKL OFF" + "POLYSOLVE_WITH_CUSOLVER OFF" + "POLYSOLVE_WITH_PARDISO OFF" + "POLYSOLVE_WITH_HYPRE OFF" + "POLYSOLVE_WITH_AMGCL OFF" + "POLYSOLVE_WITH_SPECTRA OFF" +) + +set_target_properties(polysolve PROPERTIES FOLDER third_party) +set_target_properties(polysolve_linear PROPERTIES FOLDER third_party) \ No newline at end of file diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 498167c05a..a90b4b9f2a 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -11,6 +11,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -1080,6 +1083,8 @@ bool ImageSimulationMeshTri::swap_edge_after(const Tuple& t) void ImageSimulationMeshTri::smooth_all_vertices() { + assert(m_solver); + igl::Timer timer; timer.start(); std::vector> collect_all_ops; @@ -1164,8 +1169,20 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) const Vector2d old_pos = VA[vid].m_pos; - m_vertex_attribute[vid].m_pos = - newton_method_from_stack(assembles, AMIPS2D_energy, AMIPS2D_jacobian, AMIPS2D_hessian); + // call to polysolve + { + optimization::AMIPSEnergy2D energy(assembles); + auto x = energy.initial_position(); + try { + m_solver->minimize(energy, x); + } catch (const std::exception& e) { + // polysolve might throw errors that we want to ignore (e.g., line search failed) + } + m_vertex_attribute[vid].m_pos = x; + } + + // m_vertex_attribute[vid].m_pos = + // newton_method_from_stack(assembles, AMIPS2D_energy, AMIPS2D_jacobian, AMIPS2D_hessian); wmtk::logger().trace( "old pos {} -> new pos {}", @@ -1201,12 +1218,27 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) } } } + // project to surface + //{ + // assert(m_envelope->initialized()); + // Vector2d project; + // m_envelope->nearest_point(VA[vid].m_pos, project); + // + // m_vertex_attribute[vid].m_pos = project; + //} + + // call to polysolve { - assert(m_envelope->initialized()); - Vector2d project; - m_envelope->nearest_point(VA[vid].m_pos, project); + optimization::DirichletEnergy2D energy(surface_assemble); + auto x = energy.initial_position(); + try { + m_solver->minimize(energy, x); + } catch (const std::exception& e) { + // polysolve might throw errors that we want to ignore (e.g., line search failed) + } + m_vertex_attribute[vid].m_pos = x; - m_vertex_attribute[vid].m_pos = project; + return true; // DEBUG CODE!!!! } for (auto& n : surface_assemble) { @@ -1216,18 +1248,18 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) } } - // check surface containment - for (const auto& n : surface_assemble) { - std::array edge; - for (int k = 0; k < 2; k++) { - for (int kk = 0; kk < 2; kk++) { - edge[k][kk] = n[k * 2 + kk]; - } - } - if (m_envelope->is_outside(edge)) { - return false; - } - } + //// check surface containment + // for (const auto& n : surface_assemble) { + // std::array edge; + // for (int k = 0; k < 2; k++) { + // for (int kk = 0; kk < 2; kk++) { + // edge[k][kk] = n[k * 2 + kk]; + // } + // } + // if (m_envelope->is_outside(edge)) { + // return false; + // } + // } // quality auto max_after_quality = 0.; diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 0a1caf4bc5..47b186c088 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -5,6 +5,7 @@ #include #include #include +#include // clang-format off #include @@ -102,6 +103,8 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool m_collapse_check_topology = false; // sanity check bool m_collapse_check_manifold = false; // manifoldness check after collapse + std::unique_ptr m_solver; + ImageSimulationMeshTri(Parameters& _m_params, double envelope_eps, int _num_threads = 0) : m_params(_m_params) , m_envelope_eps(envelope_eps) @@ -110,6 +113,9 @@ class ImageSimulationMeshTri : public wmtk::TriMesh p_vertex_attrs = &m_vertex_attribute; p_edge_attrs = &m_edge_attribute; p_face_attrs = &m_face_attribute; + + m_solver = optimization::create_basic_solver(); + optimization::deactivate_opt_logger(); } ~ImageSimulationMeshTri() {} diff --git a/src/wmtk/Types.hpp b/src/wmtk/Types.hpp index 7b306b1848..66b67757ec 100644 --- a/src/wmtk/Types.hpp +++ b/src/wmtk/Types.hpp @@ -13,6 +13,8 @@ using MatrixX = Eigen::Matrix; using MatrixXd = MatrixX; using MatrixXi = MatrixX; using MatrixXr = MatrixX; +using Matrix2d = Eigen::Matrix2d; +using Matrix3d = Eigen::Matrix3d; using Matrix4d = Eigen::Matrix4d; template diff --git a/src/wmtk/optimization/AMIPSEnergy.cpp b/src/wmtk/optimization/AMIPSEnergy.cpp new file mode 100644 index 0000000000..0d993545b1 --- /dev/null +++ b/src/wmtk/optimization/AMIPSEnergy.cpp @@ -0,0 +1,76 @@ +#include "AMIPSEnergy.hpp" + +#include +#include + +namespace wmtk::optimization { + +AMIPSEnergy2D::AMIPSEnergy2D(std::vector>& cells) + : m_cells(cells) +{} + +AMIPSEnergy2D::TVector AMIPSEnergy2D::initial_position() const +{ + Vector2d tmp(m_cells[0][0], m_cells[0][1]); + return tmp; +} + +double AMIPSEnergy2D::value(const TVector& x) +{ + double res = 0; + for (auto c : m_cells) { + assert(x.size() == 2); + c[0] = x[0]; + c[1] = x[1]; + res += wmtk::AMIPS2D_energy(c); + } + return res; +} + +void AMIPSEnergy2D::gradient(const TVector& x, TVector& gradv) +{ + gradv.resize(2); + gradv.setZero(); + + Vector2d tmp; + for (auto c : m_cells) { + assert(x.size() == 2); + c[0] = x[0]; + c[1] = x[1]; + wmtk::AMIPS2D_jacobian(c, tmp); + gradv += tmp; + } +} + +void AMIPSEnergy2D::hessian(const TVector& x, MatrixXd& hessian) +{ + hessian.resize(2, 2); + hessian.setZero(); + + Matrix2d tmp; + for (auto c : m_cells) { + assert(x.size() == 2); + c[0] = x[0]; + c[1] = x[1]; + wmtk::AMIPS2D_hessian(c, tmp); + hessian += tmp; + } +} + +void AMIPSEnergy2D::solution_changed(const TVector& new_x) {} + +bool AMIPSEnergy2D::is_step_valid(const TVector& x0, const TVector& x1) +{ + const Vector2d p0 = x1; + for (const auto& c : m_cells) { + const Vector2d p1(c[2], c[3]); + const Vector2d p2(c[4], c[5]); + + if (utils::orient2d(p0, p1, p2) <= 0) { + return false; + } + } + return true; +} + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/AMIPSEnergy.hpp b/src/wmtk/optimization/AMIPSEnergy.hpp new file mode 100644 index 0000000000..ac52b6cdba --- /dev/null +++ b/src/wmtk/optimization/AMIPSEnergy.hpp @@ -0,0 +1,42 @@ +#pragma once + +#include +#include + +namespace wmtk::optimization { + +class AMIPSEnergy2D : public polysolve::nonlinear::Problem +{ +public: + using typename polysolve::nonlinear::Problem::Scalar; + using typename polysolve::nonlinear::Problem::THessian; + using typename polysolve::nonlinear::Problem::TVector; + + /** + * @brief Sum of AMIPS energies for 2D triangles + * + * Each triangle must be provided as an array of 6 values: {x0, y0, x1, y1, x2, y2}. + * The first two entries (x0, y0) must be the same for all triangles and will be replaced with + * `x` during optimization. + */ + AMIPSEnergy2D(std::vector>& cells); + + TVector initial_position() const; + + 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; + +private: + std::vector> m_cells; +}; + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/DirichletEnergy.cpp b/src/wmtk/optimization/DirichletEnergy.cpp new file mode 100644 index 0000000000..98dbb35324 --- /dev/null +++ b/src/wmtk/optimization/DirichletEnergy.cpp @@ -0,0 +1,48 @@ +#include "DirichletEnergy.hpp" + +namespace wmtk::optimization { + +DirichletEnergy2D::DirichletEnergy2D(std::vector>& cells) + : m_cells(cells) +{} + +DirichletEnergy2D::TVector DirichletEnergy2D::initial_position() const +{ + Vector2d tmp(m_cells[0][0], m_cells[0][1]); + return tmp; +} + +double DirichletEnergy2D::value(const TVector& x) +{ + assert(x.size() == 2); + + double res = 0; + for (const auto& c : m_cells) { + Vector2d y(c[2], c[3]); + res += (x - y).squaredNorm(); + } + res *= 0.5; + return res; +} + +void DirichletEnergy2D::gradient(const TVector& x, TVector& gradv) +{ + assert(x.size() == 2); + + Vector2d tmp = 2 * x; + for (const auto& c : m_cells) { + Vector2d y(c[2], c[3]); + tmp -= y; + } + + gradv = tmp; +} + +void DirichletEnergy2D::hessian(const TVector& x, MatrixXd& hessian) +{ + hessian = 2 * Matrix2d::Identity(); +} + +void DirichletEnergy2D::solution_changed(const TVector& new_x) {} + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/DirichletEnergy.hpp b/src/wmtk/optimization/DirichletEnergy.hpp new file mode 100644 index 0000000000..0a09d97f97 --- /dev/null +++ b/src/wmtk/optimization/DirichletEnergy.hpp @@ -0,0 +1,40 @@ +#pragma once + +#include +#include + +namespace wmtk::optimization { + +class DirichletEnergy2D : public polysolve::nonlinear::Problem +{ +public: + using typename polysolve::nonlinear::Problem::Scalar; + using typename polysolve::nonlinear::Problem::THessian; + using typename polysolve::nonlinear::Problem::TVector; + + /** + * @brief The Dirichlet energy of a vertex position on a polyline. + * + * Each edge must be provided as an array of 4 values: {x0, y0, x1, y1}. + * The first two entries (x0, y0) must be the same for all edges and will be replaced with + * `x` during optimization. + */ + DirichletEnergy2D(std::vector>& cells); + + TVector initial_position() const; + + 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; + +private: + std::vector> m_cells; +}; + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/EnergySum.cpp b/src/wmtk/optimization/EnergySum.cpp new file mode 100644 index 0000000000..8949e9dd1d --- /dev/null +++ b/src/wmtk/optimization/EnergySum.cpp @@ -0,0 +1,59 @@ +#include "EnergySum.hpp" + +#include + +namespace wmtk::optimization { + +void EnergySum::add_energy(const std::shared_ptr& energy, const double weight) +{ + m_energies.push_back(energy); + m_weights.push_back(weight); +} + +double EnergySum::value(const TVector& x) +{ + double res = 0; + for (size_t i = 0; i < m_energies.size(); ++i) { + res += m_weights[i] * m_energies[i]->value(x); + } + return res; +} + +void EnergySum::gradient(const TVector& x, TVector& gradv) +{ + m_energies[0]->gradient(x, gradv); + gradv *= m_weights[0]; + + for (size_t i = 1; i < m_energies.size(); ++i) { + TVector g; + m_energies[i]->gradient(x, g); + gradv += m_weights[i] * g; + } +} + +void EnergySum::hessian(const TVector& x, MatrixXd& hessian) +{ + m_energies[0]->hessian(x, hessian); + hessian *= m_weights[0]; + + for (size_t i = 1; i < m_energies.size(); ++i) { + MatrixXd h; + m_energies[i]->hessian(x, h); + hessian += m_weights[i] * h; + } +} + +void EnergySum::solution_changed(const TVector& new_x) {} + +bool EnergySum::is_step_valid(const TVector& x0, const TVector& x1) +{ + for (const auto& e : m_energies) { + if (!e->is_step_valid(x0, x1)) { + return false; + } + } + + return true; +} + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/EnergySum.hpp b/src/wmtk/optimization/EnergySum.hpp new file mode 100644 index 0000000000..3d55c80721 --- /dev/null +++ b/src/wmtk/optimization/EnergySum.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include +#include + +namespace wmtk::optimization { + +class EnergySum : public polysolve::nonlinear::Problem +{ +public: + using Problem = polysolve::nonlinear::Problem; + + using typename Problem::Scalar; + using typename Problem::THessian; + using typename Problem::TVector; + + /** + * @brief A weighted sum of multiple energies. + */ + EnergySum() = default; + + /** + * @brief Add an energy term to the sum of energies. + */ + void add_energy(const std::shared_ptr& energy, const double weight = 1); + + 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; + +private: + std::vector> m_energies; + std::vector m_weights; +}; + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/solver.hpp b/src/wmtk/optimization/solver.hpp new file mode 100644 index 0000000000..841a868603 --- /dev/null +++ b/src/wmtk/optimization/solver.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include + +namespace wmtk::optimization { + +const polysolve::json basic_linear_solver_params = R"({"solver": "Eigen::LDLT"})"_json; +const polysolve::json basic_nonlinear_solver_params = + R"({"solver": "DenseNewton", "max_iterations": 10, "allow_out_of_iterations": true})"_json; + +inline std::unique_ptr create_basic_solver() +{ + return polysolve::nonlinear::Solver::create( + basic_nonlinear_solver_params, + basic_linear_solver_params, + 1, + opt_logger()); +} + +inline void deactivate_opt_logger() +{ + opt_logger().set_level(spdlog::level::off); +} + +} // namespace wmtk::optimization diff --git a/src/wmtk/utils/Rational.hpp b/src/wmtk/utils/Rational.hpp index e0285f4b25..2304f48860 100644 --- a/src/wmtk/utils/Rational.hpp +++ b/src/wmtk/utils/Rational.hpp @@ -11,7 +11,7 @@ class Rational public: mpq_t value; void canonicalize() { mpq_canonicalize(value); } - int get_sign() { return mpq_sgn(value); } + int get_sign() const { return mpq_sgn(value); } template void init(const T& v) { diff --git a/src/wmtk/utils/orient.cpp b/src/wmtk/utils/orient.cpp new file mode 100644 index 0000000000..83138598d7 --- /dev/null +++ b/src/wmtk/utils/orient.cpp @@ -0,0 +1,75 @@ +#include "orient.hpp" + +#include +#include + +namespace wmtk::utils { + +bool orient3d(const Vector3r& p0, const Vector3r& p1, const Vector3r& p2, const Vector3r& p3) +{ + const Vector3r a(p1 - p0); + const Vector3r b(p2 - p0); + // Vector3r n = a.cross(b); + // Vector3r d = p3 - p0; + // auto res = n.dot(d); + + // cross product + const Rational nx = a[1] * b[2] - a[2] * b[1]; + const Rational ny = a[2] * b[0] - a[0] * b[2]; + const Rational nz = a[0] * b[1] - a[1] * b[0]; + + const Vector3r d = p3 - p0; + + // dot product: n · d + Rational res = nx * d[0] + ny * d[1] + nz * d[2]; + + + if (res > 0) { + // predicates returns pos value: non-inverted + return true; + } else { + return false; + } +} + +bool orient3d(const Vector3d& p0, const Vector3d& p1, const Vector3d& p2, const Vector3d& p3) +{ + igl::predicates::exactinit(); + auto res = igl::predicates::orient3d(p0, p1, p2, p3); + int result; + if (res == igl::predicates::Orientation::POSITIVE) { + result = 1; + } else if (res == igl::predicates::Orientation::NEGATIVE) { + result = -1; + } else { + result = 0; + } + + if (result < 0) { + // neg result == pos tet (tet origin from geogram delaunay) + return true; + } + return false; +} + +bool orient2d(const Vector2r& p0, const Vector2r& p1, const Vector2r& p2) +{ + Eigen::Matrix2 M; + M.row(0) = p1 - p0; + M.row(1) = p2 - p0; + const Rational det = M(0, 0) * M(1, 1) - M(0, 1) * M(1, 0); + // logger().info("{}", det.to_double()); + return det.get_sign() == 1; +} + +bool orient2d(const Vector2d& p0, const Vector2d& p1, const Vector2d& p2) +{ + igl::predicates::exactinit(); + auto res = igl::predicates::orient2d(p0, p1, p2); + if (res == igl::predicates::Orientation::POSITIVE) { + return true; + } + return false; +} + +} // namespace wmtk::utils \ No newline at end of file diff --git a/src/wmtk/utils/orient.hpp b/src/wmtk/utils/orient.hpp new file mode 100644 index 0000000000..61c8ec7fa6 --- /dev/null +++ b/src/wmtk/utils/orient.hpp @@ -0,0 +1,15 @@ +#pragma once + +#include + +namespace wmtk::utils { + +bool orient3d(const Vector3r& p0, const Vector3r& p1, const Vector3r& p2, const Vector3r& p3); + +bool orient3d(const Vector3d& p0, const Vector3d& p1, const Vector3d& p2, const Vector3d& p3); + +bool orient2d(const Vector2r& p0, const Vector2r& p1, const Vector2r& p2); + +bool orient2d(const Vector2d& p0, const Vector2d& p1, const Vector2d& p2); + +} // namespace wmtk::utils \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e3f813a294..f744f73408 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -16,6 +16,8 @@ set(TEST_SOURCES test_raw_simplex_collection.cpp test_simplex_navigation.cpp test_tuple.cpp + test_orient.cpp + test_optimization.cpp ) add_executable(wmtk_tests ${TEST_SOURCES}) diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp new file mode 100644 index 0000000000..ca70df5a9c --- /dev/null +++ b/tests/test_optimization.cpp @@ -0,0 +1,153 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace wmtk; + +TEST_CASE("amips_energy_2d", "[energies]") +{ + const Vector2d p0(0.5, 1e-6); + std::vector> cells; + cells.push_back({{p0[0], p0[1], 0, 0, 1, 0}}); + + optimization::AMIPSEnergy2D energy(cells); + CHECK(energy.is_step_valid(p0, p0)); + CHECK_FALSE(energy.is_step_valid(p0, -p0)); + + { + VectorXd g; + energy.gradient(p0, g); + CHECK(g[1] < 0); + + MatrixXd h; + CHECK_NOTHROW(energy.hessian(p0, h)); + } + auto x = energy.initial_position(); + const double e_before = energy.value(x); + { + auto m_solver = optimization::create_basic_solver(); + optimization::deactivate_opt_logger(); + + CHECK_NOTHROW(m_solver->minimize(energy, x)); + } + const double e_after = energy.value(x); + CHECK(e_after < e_before); + CHECK(e_after < 2 + 1e-6); // should have perfect quality, i.e. AMIPS=2 +} + +TEST_CASE("dirichlet_energy_2d", "[energies]") +{ + const Vector2d p0(0.5, 1); + std::vector> cells; + cells.push_back({{p0[0], p0[1], 0, 0}}); + cells.push_back({{p0[0], p0[1], 1, 0}}); + + optimization::DirichletEnergy2D energy(cells); + CHECK(energy.is_step_valid(p0, p0)); + + { + VectorXd g; + energy.gradient(p0, g); + CHECK(g[1] > 0); + + MatrixXd h; + CHECK_NOTHROW(energy.hessian(p0, h)); + } + auto x = energy.initial_position(); + const double e_before = energy.value(x); + { + auto m_solver = optimization::create_basic_solver(); + optimization::deactivate_opt_logger(); + + CHECK_NOTHROW(m_solver->minimize(energy, x)); + } + const double e_after = energy.value(x); + CHECK(e_after < e_before); +} + +TEST_CASE("amips_plus_dirichlet_energy_2d", "[energies]") +{ + const auto VF = utils::examples::tri::edge_region(); + MatrixXd V = VF.V; + MatrixXi F = VF.F; + + auto write = [&V, &F]() { + return; // comment out for debug print + static int wcount = 0; + MatrixXd V3; + V3.resize(V.rows(), 3); + V3.setZero(); + V3.block(0, 0, V.rows(), V.cols()) = V; + igl::writeOFF(fmt::format("debug_{}.off", wcount++), V3, F); + }; + + V.row(4) = Vector2d(1.8, 0); + + // make vertex 4 the first one in all triangles + for (int i = 0; i < F.rows(); ++i) { + int j = 0; + for (; j < 3; ++j) { + if (F(i, j) == 4) { + break; + } + } + if (j < 3) { + Vector3i f; + f[0] = F(i, j); + f[1] = F(i, (j + 1) % 3); + f[2] = F(i, (j + 2) % 3); + F.row(i) = f; + } + } + // logger().info("F:\n{}", F); + + // collect all triangles that contain vertex 4 + std::vector> amips_cells; + for (int i = 0; i < F.rows(); ++i) { + if (F(i, 0) != 4) { + continue; + } + std::array c; + for (size_t j = 0; j < 3; ++j) { + c[2 * j + 0] = V(F(i, j), 0); + c[2 * j + 1] = V(F(i, j), 1); + } + amips_cells.push_back(c); + } + + std::vector> smoothing_cells; + smoothing_cells.push_back({{V(4, 0), V(4, 1), V(5, 0), V(5, 1)}}); + smoothing_cells.push_back({{V(4, 0), V(4, 1), V(7, 0), V(7, 1)}}); + + write(); + + auto amips_energy = std::make_shared(amips_cells); + auto smooth_energy = std::make_shared(smoothing_cells); + + auto total_energy = std::make_shared(); + total_energy->add_energy(amips_energy); + total_energy->add_energy(smooth_energy, 10); + + auto x = amips_energy->initial_position(); + const double e_before = amips_energy->value(x); + { + auto m_solver = optimization::create_basic_solver(); + optimization::deactivate_opt_logger(); + // CHECK_NOTHROW(m_solver->minimize(*amips_energy, x)); + // CHECK_NOTHROW(m_solver->minimize(*smooth_energy, x)); + CHECK_NOTHROW(m_solver->minimize(*total_energy, x)); + } + const double e_after = amips_energy->value(x); + // logger().info("before: {:.4}, after: {:.4}", e_before, e_after); + CHECK(e_after < e_before); + V.row(4) = x; + write(); +} diff --git a/tests/test_orient.cpp b/tests/test_orient.cpp new file mode 100644 index 0000000000..411b5a26bf --- /dev/null +++ b/tests/test_orient.cpp @@ -0,0 +1,92 @@ +#include + +#include +#include + +using namespace wmtk; + +TEST_CASE("orient2d_double", "[orient][utils]") +{ + const Vector2d p0(0, 0); + const Vector2d p1(1, 0); + SECTION("positive") + { + const Vector2d p2(0.5, 1e-12); + CHECK(utils::orient2d(p0, p1, p2)); + } + SECTION("degenerate") + { + const Vector2d p2(0.5, 0); + CHECK_FALSE(utils::orient2d(p0, p1, p2)); + } + SECTION("negative") + { + const Vector2d p2(0.5, -1e-12); + CHECK_FALSE(utils::orient2d(p0, p1, p2)); + } +} + +TEST_CASE("orient2d_rational", "[orient][utils]") +{ + const Vector2r p0(0, 0); + const Vector2r p1(1, 0); + SECTION("positive") + { + const Vector2r p2(0.5, 1e-12); + CHECK(utils::orient2d(p0, p1, p2)); + } + SECTION("degenerate") + { + const Vector2r p2(0.5, 0); + CHECK_FALSE(utils::orient2d(p0, p1, p2)); + } + SECTION("negative") + { + const Vector2r p2(0.5, -1e-12); + CHECK_FALSE(utils::orient2d(p0, p1, p2)); + } +} + +TEST_CASE("orient3d_double", "[orient][utils]") +{ + const Vector3d p0(0, 0, 0); + const Vector3d p1(1, 0, 0); + const Vector3d p2(0.5, 1, 0); + SECTION("positive") + { + const Vector3d p3(0.5, 0.5, 1e-12); + CHECK(utils::orient3d(p0, p1, p2, p3)); + } + SECTION("degenerate") + { + const Vector3d p3(0.5, 0.5, 0); + CHECK_FALSE(utils::orient3d(p0, p1, p2, p3)); + } + SECTION("negative") + { + const Vector3d p3(0.5, 0.5, -1e-12); + CHECK_FALSE(utils::orient3d(p0, p1, p2, p3)); + } +} + +TEST_CASE("orient3d_rational", "[orient][utils]") +{ + const Vector3r p0(0, 0, 0); + const Vector3r p1(1, 0, 0); + const Vector3r p2(0.5, 1, 0); + SECTION("positive") + { + const Vector3r p3(0.5, 0.5, 1e-12); + CHECK(utils::orient3d(p0, p1, p2, p3)); + } + SECTION("degenerate") + { + const Vector3r p3(0.5, 0.5, 0); + CHECK_FALSE(utils::orient3d(p0, p1, p2, p3)); + } + SECTION("negative") + { + const Vector3r p3(0.5, 0.5, -1e-12); + CHECK_FALSE(utils::orient3d(p0, p1, p2, p3)); + } +} \ No newline at end of file From 9bb936af58436ea6c5c3eabd75fab57fb7324315 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 11 Mar 2026 09:52:01 +0100 Subject: [PATCH 02/19] Add envelope energy. --- .../ImageSimulationMeshTri.cpp | 95 +++++++++++-------- src/wmtk/envelope/Envelope.cpp | 8 ++ src/wmtk/envelope/Envelope.hpp | 1 + src/wmtk/optimization/EnvelopeEnergy.cpp | 55 +++++++++++ src/wmtk/optimization/EnvelopeEnergy.hpp | 41 ++++++++ 5 files changed, 158 insertions(+), 42 deletions(-) create mode 100644 src/wmtk/optimization/EnvelopeEnergy.cpp create mode 100644 src/wmtk/optimization/EnvelopeEnergy.hpp diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index a90b4b9f2a..d96394999d 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -13,6 +13,8 @@ #include #include #include +#include +#include #include #include #include @@ -1170,16 +1172,18 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) const Vector2d old_pos = VA[vid].m_pos; // call to polysolve - { - optimization::AMIPSEnergy2D energy(assembles); - auto x = energy.initial_position(); - try { - m_solver->minimize(energy, x); - } catch (const std::exception& e) { - // polysolve might throw errors that we want to ignore (e.g., line search failed) - } - m_vertex_attribute[vid].m_pos = x; - } + std::shared_ptr total_energy; + auto amips_energy = std::make_shared(assembles); + total_energy = amips_energy; + //{ + // auto x = amips_energy->initial_position(); + // try { + // m_solver->minimize(*amips_energy, x); + // } catch (const std::exception& e) { + // // polysolve might throw errors that we want to ignore (e.g., line search failed) + // } + // m_vertex_attribute[vid].m_pos = x; + //} // m_vertex_attribute[vid].m_pos = // newton_method_from_stack(assembles, AMIPS2D_energy, AMIPS2D_jacobian, AMIPS2D_hessian); @@ -1227,41 +1231,46 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) // m_vertex_attribute[vid].m_pos = project; //} - // call to polysolve - { - optimization::DirichletEnergy2D energy(surface_assemble); - auto x = energy.initial_position(); - try { - m_solver->minimize(energy, x); - } catch (const std::exception& e) { - // polysolve might throw errors that we want to ignore (e.g., line search failed) - } - m_vertex_attribute[vid].m_pos = x; + auto smooth_energy = std::make_shared(surface_assemble); + auto envelope_energy = + std::make_shared(m_envelope, surface_assemble); + auto energy_sum = std::make_shared(); + energy_sum->add_energy(amips_energy); + energy_sum->add_energy(smooth_energy, 10); + energy_sum->add_energy(envelope_energy, 100); + total_energy = energy_sum; + } + + { + auto x = amips_energy->initial_position(); + try { + m_solver->minimize(*total_energy, x); + } catch (const std::exception&) { + // polysolve might throw errors that we want to ignore (e.g., line search failed) + } + m_vertex_attribute[vid].m_pos = x; + } - return true; // DEBUG CODE!!!! + for (auto& n : surface_assemble) { + for (int kk = 0; kk < 2; kk++) { + n[kk] = VA[vid].m_pos[kk]; } + } - for (auto& n : surface_assemble) { + // check surface containment + for (const auto& n : surface_assemble) { + std::array edge; + for (int k = 0; k < 2; k++) { for (int kk = 0; kk < 2; kk++) { - n[kk] = VA[vid].m_pos[kk]; + edge[k][kk] = n[k * 2 + kk]; } } + if (m_envelope->is_outside(edge)) { + return false; + } } - //// check surface containment - // for (const auto& n : surface_assemble) { - // std::array edge; - // for (int k = 0; k < 2; k++) { - // for (int kk = 0; kk < 2; kk++) { - // edge[k][kk] = n[k * 2 + kk]; - // } - // } - // if (m_envelope->is_outside(edge)) { - // return false; - // } - // } - - // quality + // quality (only check if not on surface) auto max_after_quality = 0.; for (const size_t fid : locs) { if (is_inverted(fid)) { @@ -1271,8 +1280,10 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) m_face_attribute[fid].m_quality = q; max_after_quality = std::max(max_after_quality, q); } - if (max_after_quality > max_quality) { - return false; + if (!VA[vid].m_is_on_surface) { + if (max_after_quality > max_quality) { + return false; + } } return true; @@ -1515,9 +1526,9 @@ std::tuple ImageSimulationMeshTri::local_operations( for (int n = 0; n < ops[i]; n++) { wmtk::logger().info("==smoothing {}==", n); smooth_all_vertices(); - } - if (m_params.debug_output) { - write_vtu(fmt::format("debug_{}", debug_print_counter++)); + if (m_params.debug_output) { + write_vtu(fmt::format("debug_{}", debug_print_counter++)); + } } auto [max_energy, avg_energy] = get_max_avg_energy(); wmtk::logger().info("smooth max energy = {:.6} avg = {:.6}", max_energy, avg_energy); diff --git a/src/wmtk/envelope/Envelope.cpp b/src/wmtk/envelope/Envelope.cpp index 4514421ada..0e85c47473 100644 --- a/src/wmtk/envelope/Envelope.cpp +++ b/src/wmtk/envelope/Envelope.cpp @@ -312,4 +312,12 @@ double SampleEnvelope::squared_distance(const Eigen::Vector3d& p) const return d2; } +double SampleEnvelope::squared_distance(const Eigen::Vector2d& p) const +{ + double d2; + SimpleBVH::Vector2d out; + m_bvh->nearest_facet(p, out, d2); + return d2; +} + } // namespace wmtk \ No newline at end of file diff --git a/src/wmtk/envelope/Envelope.hpp b/src/wmtk/envelope/Envelope.hpp index 63495f3dfa..a5fc416870 100644 --- a/src/wmtk/envelope/Envelope.hpp +++ b/src/wmtk/envelope/Envelope.hpp @@ -82,6 +82,7 @@ class SampleEnvelope : public wmtk::Envelope bool initialized() { return m_bvh != nullptr; }; double squared_distance(const Eigen::Vector3d& p) const; + double squared_distance(const Eigen::Vector2d& p) const; private: std::vector geo_vertex_ind; diff --git a/src/wmtk/optimization/EnvelopeEnergy.cpp b/src/wmtk/optimization/EnvelopeEnergy.cpp new file mode 100644 index 0000000000..1549b2275d --- /dev/null +++ b/src/wmtk/optimization/EnvelopeEnergy.cpp @@ -0,0 +1,55 @@ +#include "EnvelopeEnergy.hpp" + +namespace wmtk::optimization { + +EnvelopeEnergy2D::EnvelopeEnergy2D( + const std::shared_ptr& envelope, + const std::vector>& cells) + : m_envelope(envelope) + , m_cells(cells) +{ + assert(m_envelope); +} + +double EnvelopeEnergy2D::value(const TVector& x) +{ + assert(x.size() == 2); + Vector2d r(x); + return m_envelope->squared_distance(r); +} + +void EnvelopeEnergy2D::gradient(const TVector& x, TVector& gradv) +{ + assert(x.size() == 2); + Vector2d r(x); + Vector2d n; + m_envelope->nearest_point(r, n); + gradv = r - n; +} + +void EnvelopeEnergy2D::hessian(const TVector& x, MatrixXd& hessian) +{ + hessian = 2 * Matrix2d::Identity(); +} + +void EnvelopeEnergy2D::solution_changed(const TVector& new_x) {} + +bool EnvelopeEnergy2D::is_step_valid(const TVector& x0, const TVector& x1) +{ + Vector2d r(x1); + if (m_envelope->is_outside(r)) { + return false; + } + // for (const auto& cell : m_cells) { + // std::array edge; + // edge[0] = r; + // edge[1] = Vector2d(cell[2], cell[3]); + // if (m_envelope->is_outside(edge)) { + // return false; + // } + // } + + return true; +} + +} // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/EnvelopeEnergy.hpp b/src/wmtk/optimization/EnvelopeEnergy.hpp new file mode 100644 index 0000000000..df8c5c936d --- /dev/null +++ b/src/wmtk/optimization/EnvelopeEnergy.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include +#include +#include + +namespace wmtk::optimization { + +class EnvelopeEnergy2D : public polysolve::nonlinear::Problem +{ +public: + using typename polysolve::nonlinear::Problem::Scalar; + using typename polysolve::nonlinear::Problem::THessian; + using typename polysolve::nonlinear::Problem::TVector; + + /** + * @brief The energy is the squared distance to an envelope. + * + */ + EnvelopeEnergy2D( + const std::shared_ptr& envelope, + const std::vector>& cells); + + 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; + +private: + std::shared_ptr m_envelope; + std::vector> m_cells; +}; + +} // namespace wmtk::optimization \ No newline at end of file From af4f68db2e875b54c8d5f28b08021e1d6913183d Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 11 Mar 2026 10:04:31 +0100 Subject: [PATCH 03/19] Do not collapse away from surface. --- .../components/image_simulation/ImageSimulationMeshTri.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index d96394999d..26feffc768 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -771,8 +771,11 @@ bool ImageSimulationMeshTri::collapse_edge_before(const Tuple& loc) // surface if (cache.edge_length > 0 && VA[v1_id].m_is_on_surface) { - if (!VA[v2_id].m_is_on_surface && m_envelope->is_outside(VA[v2_id].m_pos)) { - return false; + // if (!VA[v2_id].m_is_on_surface && m_envelope->is_outside(VA[v2_id].m_pos)) { + // return false; + // } + if (!VA[v2_id].m_is_on_surface) { + return false; // do not collapse away from surface } } From 7e99283ce154a51385069433f85690164b9828f2 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Fri, 13 Mar 2026 17:36:09 +0100 Subject: [PATCH 04/19] Implementing the smoothing energy. The construction is not final and needs to be adapted. --- .../ImageSimulationMeshTri.cpp | 21 +++++ .../ImageSimulationMeshTri.hpp | 1 + src/wmtk/optimization/DirichletEnergy.cpp | 89 ++++++++++++++++++- src/wmtk/optimization/DirichletEnergy.hpp | 44 ++++++++- tests/test_optimization.cpp | 36 ++++++++ 5 files changed, 189 insertions(+), 2 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 26feffc768..2380f907a3 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -487,6 +487,27 @@ std::vector> ImageSimulationMeshTri::get_edges_by_conditio void ImageSimulationMeshTri::split_all_edges() { + // build mass-matrix + { + const auto vs = get_vertices(); + m_surface_mass.resize(vert_capacity()); + for (const Tuple& t : vs) { + const size_t vid = t.vid(*this); + if (!m_vertex_attribute.at(vid).m_is_on_surface) { + continue; + } + const auto es = get_order1_edges_for_vertex(vid); + double mass = 0; + for (const simplex::Edge& e : es.edges()) { + const Vector2d& p0 = m_vertex_attribute.at(e.vertices()[0]).m_pos; + const Vector2d& p1 = m_vertex_attribute.at(e.vertices()[1]).m_pos; + mass += 0.5 * (p1 - p0).norm(); + } + m_surface_mass[vid] = mass; + } + } + + igl::Timer timer; double time; timer.start(); diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 47b186c088..c6cdb18bd4 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -104,6 +104,7 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool m_collapse_check_manifold = false; // manifoldness check after collapse std::unique_ptr m_solver; + std::vector m_surface_mass; // the mass matrix for surface vertices ImageSimulationMeshTri(Parameters& _m_params, double envelope_eps, int _num_threads = 0) : m_params(_m_params) diff --git a/src/wmtk/optimization/DirichletEnergy.cpp b/src/wmtk/optimization/DirichletEnergy.cpp index 98dbb35324..35a3499649 100644 --- a/src/wmtk/optimization/DirichletEnergy.cpp +++ b/src/wmtk/optimization/DirichletEnergy.cpp @@ -43,6 +43,93 @@ void DirichletEnergy2D::hessian(const TVector& x, MatrixXd& hessian) hessian = 2 * Matrix2d::Identity(); } -void DirichletEnergy2D::solution_changed(const TVector& new_x) {} + +SmoothingEnergy2D::SmoothingEnergy2D(std::vector> m_cells) +{ + assert(m_cells.size() == 2); + + const Vector2d p0(m_cells[0][2], m_cells[0][3]); + const Vector2d p1(m_cells[0][0], m_cells[0][1]); // this vertex is optimized + const Vector2d p2(m_cells[1][2], m_cells[1][3]); + + m_x.resize(3, 2); + m_x.row(0) = p0; + m_x.row(1) = p1; + m_x.row(2) = p2; + + // uniform laplacian + m_M = 2; + m_M_inv = 1 / m_M; + m_L_w[0] = 1; + m_L_w[1] = -2; + m_L_w[2] = 1; + + m_LTML_row1 = m_M_inv * m_L_w[1] * m_L_w; +} + +SmoothingEnergy2D::TVector SmoothingEnergy2D::initial_position() const +{ + return m_x.row(1); +} + +double SmoothingEnergy2D::value(const TVector& x) +{ + assert(x.size() == 2); + double energy = 0; + for (size_t i = 0; i < 2; ++i) { + Vector3d v(m_x(0, i), x[i], m_x(2, i)); + // energy += v.transpose() * m_LTML * v; + double Lwv = m_L_w.dot(v); + energy += m_M_inv * Lwv * Lwv; + } + + return energy; +} + +void SmoothingEnergy2D::gradient(const TVector& x, TVector& gradv) +{ + assert(x.size() == 2); + gradv.resize(2); + + for (size_t i = 0; i < 2; ++i) { + Vector3d v(m_x(0, i), x[i], m_x(2, i)); + gradv[i] = 2 * m_LTML_row1.dot(v); + } +} + +void SmoothingEnergy2D::hessian(const TVector& x, MatrixXd& hessian) +{ + assert(x.size() == 2); + hessian = Matrix2d::Identity() * 2 * m_LTML_row1[1]; +} + +void SmoothingEnergy2D::add_mass_and_stiffness_matrix(const double& M, const Vector3d& L_w) +{ + m_M = M; + m_M_inv = 1 / M; + m_L_w = L_w; + m_LTML_row1 = m_M_inv * m_L_w[1] * m_L_w; +} + +void SmoothingEnergy2D::compute_local_mass_and_stiffness( + std::vector> m_cells, + double& M, + Vector3d& L_w) +{ + assert(m_cells.size() == 2); + + const Vector2d p0(m_cells[0][2], m_cells[0][3]); + const Vector2d p1(m_cells[0][0], m_cells[0][1]); // this vertex is optimized + const Vector2d p2(m_cells[1][2], m_cells[1][3]); + + const double e0 = (p1 - p0).norm(); + const double e2 = (p2 - p1).norm(); + + M = 0.5 * (e0 + e2); + + L_w[0] = 1 / e0; + L_w[1] = -(1 / e0 + 1 / e2); + L_w[2] = 1 / e2; +} } // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/DirichletEnergy.hpp b/src/wmtk/optimization/DirichletEnergy.hpp index 0a09d97f97..2ac3f2141b 100644 --- a/src/wmtk/optimization/DirichletEnergy.hpp +++ b/src/wmtk/optimization/DirichletEnergy.hpp @@ -31,10 +31,52 @@ class DirichletEnergy2D : public polysolve::nonlinear::Problem } void hessian(const TVector& x, MatrixXd& hessian) override; - void solution_changed(const TVector& new_x) override; + void solution_changed(const TVector& new_x) override {} private: std::vector> m_cells; }; +class SmoothingEnergy2D : public polysolve::nonlinear::Problem +{ +public: + using typename polysolve::nonlinear::Problem::Scalar; + using typename polysolve::nonlinear::Problem::THessian; + using typename polysolve::nonlinear::Problem::TVector; + + /** + * @brief The smoothing energy of a vertex position on a polyline. + * + * Three positions must be provided, the optimized position (p1) and its two neighbors (p0, p2). + */ + SmoothingEnergy2D(std::vector> m_cells); + + TVector initial_position() const; + + 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 {} + + void add_mass_and_stiffness_matrix(const double& M, const Vector3d& L_w); + + static void compute_local_mass_and_stiffness( + std::vector> m_cells, + double& M, + Vector3d& L_w); + +private: + MatrixXd m_x; + Vector3d m_L_w; + double m_M; + double m_M_inv; + + Vector3d m_LTML_row1; // derived from the other matrices +}; + } // namespace wmtk::optimization \ No newline at end of file diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp index ca70df5a9c..549f589f24 100644 --- a/tests/test_optimization.cpp +++ b/tests/test_optimization.cpp @@ -151,3 +151,39 @@ TEST_CASE("amips_plus_dirichlet_energy_2d", "[energies]") V.row(4) = x; write(); } + +TEST_CASE("smoothing_energy_2d", "[energies]") +{ + const Vector2d p0(0.8, 1); + std::vector> cells; + cells.push_back({{p0[0], p0[1], 0, 0}}); + cells.push_back({{p0[0], p0[1], 1, 0}}); + + optimization::SmoothingEnergy2D energy(cells); + double M; + Vector3d L_w; + optimization::SmoothingEnergy2D::compute_local_mass_and_stiffness(cells, M, L_w); + energy.add_mass_and_stiffness_matrix(M, L_w); + + CHECK(energy.is_step_valid(p0, p0)); + + { + VectorXd g; + energy.gradient(p0, g); + CHECK(g[1] > 0); + + MatrixXd h; + CHECK_NOTHROW(energy.hessian(p0, h)); + } + auto x = energy.initial_position(); + const double e_before = energy.value(x); + { + auto m_solver = optimization::create_basic_solver(); + optimization::deactivate_opt_logger(); + + CHECK_NOTHROW(m_solver->minimize(energy, x)); + } + const double e_after = energy.value(x); + CHECK(e_after < e_before); + CHECK(e_after < 1e-5); +} \ No newline at end of file From cd5ca248146de7c8552e5bcb205e64bd607149db Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Mon, 16 Mar 2026 10:54:36 +0100 Subject: [PATCH 05/19] Change the smoothing energy such that it takes the entries of a global mass and stiffness matrix. --- .../ImageSimulationMeshTri.cpp | 76 +++++++++++------ .../ImageSimulationMeshTri.hpp | 1 + src/wmtk/optimization/DirichletEnergy.cpp | 81 ++++++++----------- src/wmtk/optimization/DirichletEnergy.hpp | 44 +++++++--- tests/test_optimization.cpp | 12 ++- 5 files changed, 123 insertions(+), 91 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 2380f907a3..7b0b54c21a 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -487,27 +487,6 @@ std::vector> ImageSimulationMeshTri::get_edges_by_conditio void ImageSimulationMeshTri::split_all_edges() { - // build mass-matrix - { - const auto vs = get_vertices(); - m_surface_mass.resize(vert_capacity()); - for (const Tuple& t : vs) { - const size_t vid = t.vid(*this); - if (!m_vertex_attribute.at(vid).m_is_on_surface) { - continue; - } - const auto es = get_order1_edges_for_vertex(vid); - double mass = 0; - for (const simplex::Edge& e : es.edges()) { - const Vector2d& p0 = m_vertex_attribute.at(e.vertices()[0]).m_pos; - const Vector2d& p1 = m_vertex_attribute.at(e.vertices()[1]).m_pos; - mass += 0.5 * (p1 - p0).norm(); - } - m_surface_mass[vid] = mass; - } - } - - igl::Timer timer; double time; timer.start(); @@ -1111,6 +1090,37 @@ void ImageSimulationMeshTri::smooth_all_vertices() { assert(m_solver); + // build mass-matrix + { + const auto vs = get_vertices(); + m_surface_mass.resize(vert_capacity()); + m_surface_stiffness.resize(vert_capacity()); + for (const Tuple& t : vs) { + const size_t vid = t.vid(*this); + if (!m_vertex_attribute.at(vid).m_is_on_surface) { + continue; + } + const auto es = get_order1_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::SmoothingEnergy2D::local_mass_and_stiffness(pts, M, L_w); + // optimization::SmoothingEnergy2D::uniform_mass_and_stiffness(pts, M, L_w); + } + } + igl::Timer timer; timer.start(); std::vector> collect_all_ops; @@ -1218,7 +1228,21 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) VA[vid].m_pos.transpose()); std::vector> surface_assemble; + std::array surface_pts; if (VA[vid].m_is_on_surface) { + const auto es = get_order1_edges_for_vertex(vid); + if (es.size() != 2) { + return false; // can only smooth vertices with 2 neighbors + } + + surface_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]; + surface_pts[i + 1] = m_vertex_attribute.at(neighbor_id).m_pos; + } + std::set unique_eid; for (const size_t fid : locs) { for (size_t j = 0; j < 3; j++) { @@ -1255,13 +1279,17 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) // m_vertex_attribute[vid].m_pos = project; //} - auto smooth_energy = std::make_shared(surface_assemble); + const auto& M = m_surface_mass[vid]; + const auto& L_w = m_surface_stiffness[vid]; + + // auto smooth_energy = std::make_shared(surface_assemble); + auto smooth_energy = std::make_shared(surface_pts, M, L_w); auto envelope_energy = std::make_shared(m_envelope, surface_assemble); auto energy_sum = std::make_shared(); energy_sum->add_energy(amips_energy); - energy_sum->add_energy(smooth_energy, 10); - energy_sum->add_energy(envelope_energy, 100); + energy_sum->add_energy(smooth_energy, 1e2); + energy_sum->add_energy(envelope_energy, 1e2 * M); total_energy = energy_sum; } diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index c6cdb18bd4..31923b05a7 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -105,6 +105,7 @@ class ImageSimulationMeshTri : public wmtk::TriMesh std::unique_ptr m_solver; std::vector m_surface_mass; // the mass matrix for surface vertices + std::vector m_surface_stiffness; // stiffness matrix for surface vertices ImageSimulationMeshTri(Parameters& _m_params, double envelope_eps, int _num_threads = 0) : m_params(_m_params) diff --git a/src/wmtk/optimization/DirichletEnergy.cpp b/src/wmtk/optimization/DirichletEnergy.cpp index 35a3499649..8ea3e29499 100644 --- a/src/wmtk/optimization/DirichletEnergy.cpp +++ b/src/wmtk/optimization/DirichletEnergy.cpp @@ -44,32 +44,21 @@ void DirichletEnergy2D::hessian(const TVector& x, MatrixXd& hessian) } -SmoothingEnergy2D::SmoothingEnergy2D(std::vector> m_cells) +SmoothingEnergy2D::SmoothingEnergy2D( + const std::array& pts, + const double& M, + const Vector3d& L_w) + : m_pts(pts) + , m_M(M) + , m_M_inv(1. / M) + , m_L_w_row(L_w) { - assert(m_cells.size() == 2); - - const Vector2d p0(m_cells[0][2], m_cells[0][3]); - const Vector2d p1(m_cells[0][0], m_cells[0][1]); // this vertex is optimized - const Vector2d p2(m_cells[1][2], m_cells[1][3]); - - m_x.resize(3, 2); - m_x.row(0) = p0; - m_x.row(1) = p1; - m_x.row(2) = p2; - - // uniform laplacian - m_M = 2; - m_M_inv = 1 / m_M; - m_L_w[0] = 1; - m_L_w[1] = -2; - m_L_w[2] = 1; - - m_LTML_row1 = m_M_inv * m_L_w[1] * m_L_w; + m_LTML_row = m_M_inv * m_L_w_row[0] * m_L_w_row; } SmoothingEnergy2D::TVector SmoothingEnergy2D::initial_position() const { - return m_x.row(1); + return m_pts[0]; } double SmoothingEnergy2D::value(const TVector& x) @@ -77,9 +66,8 @@ double SmoothingEnergy2D::value(const TVector& x) assert(x.size() == 2); double energy = 0; for (size_t i = 0; i < 2; ++i) { - Vector3d v(m_x(0, i), x[i], m_x(2, i)); - // energy += v.transpose() * m_LTML * v; - double Lwv = m_L_w.dot(v); + Vector3d v(x[i], m_pts[1][i], m_pts[2][i]); + double Lwv = m_L_w_row.dot(v); energy += m_M_inv * Lwv * Lwv; } @@ -92,44 +80,41 @@ void SmoothingEnergy2D::gradient(const TVector& x, TVector& gradv) gradv.resize(2); for (size_t i = 0; i < 2; ++i) { - Vector3d v(m_x(0, i), x[i], m_x(2, i)); - gradv[i] = 2 * m_LTML_row1.dot(v); + Vector3d v(x[i], m_pts[1][i], m_pts[2][i]); + gradv[i] = 2 * m_LTML_row.dot(v); } } void SmoothingEnergy2D::hessian(const TVector& x, MatrixXd& hessian) { assert(x.size() == 2); - hessian = Matrix2d::Identity() * 2 * m_LTML_row1[1]; -} - -void SmoothingEnergy2D::add_mass_and_stiffness_matrix(const double& M, const Vector3d& L_w) -{ - m_M = M; - m_M_inv = 1 / M; - m_L_w = L_w; - m_LTML_row1 = m_M_inv * m_L_w[1] * m_L_w; + hessian = Matrix2d::Identity() * 2 * m_LTML_row[0]; } -void SmoothingEnergy2D::compute_local_mass_and_stiffness( - std::vector> m_cells, +void SmoothingEnergy2D::local_mass_and_stiffness( + const std::array& pts, double& M, Vector3d& L_w) { - assert(m_cells.size() == 2); + const double e1 = (pts[1] - pts[0]).norm(); + const double e2 = (pts[2] - pts[0]).norm(); - const Vector2d p0(m_cells[0][2], m_cells[0][3]); - const Vector2d p1(m_cells[0][0], m_cells[0][1]); // this vertex is optimized - const Vector2d p2(m_cells[1][2], m_cells[1][3]); + M = 0.5 * (e1 + e2); - const double e0 = (p1 - p0).norm(); - const double e2 = (p2 - p1).norm(); - - M = 0.5 * (e0 + e2); - - L_w[0] = 1 / e0; - L_w[1] = -(1 / e0 + 1 / e2); + L_w[0] = -(1 / e1 + 1 / e2); + L_w[1] = 1 / e1; L_w[2] = 1 / e2; } +void SmoothingEnergy2D::uniform_mass_and_stiffness( + const std::array& pts, + double& M, + Vector3d& L_w) +{ + M = 2; + L_w[0] = -2; + L_w[1] = 1; + L_w[2] = 1; +} + } // namespace wmtk::optimization \ No newline at end of file diff --git a/src/wmtk/optimization/DirichletEnergy.hpp b/src/wmtk/optimization/DirichletEnergy.hpp index 2ac3f2141b..a93592c52e 100644 --- a/src/wmtk/optimization/DirichletEnergy.hpp +++ b/src/wmtk/optimization/DirichletEnergy.hpp @@ -47,9 +47,9 @@ class SmoothingEnergy2D : public polysolve::nonlinear::Problem /** * @brief The smoothing energy of a vertex position on a polyline. * - * Three positions must be provided, the optimized position (p1) and its two neighbors (p0, p2). + * Three positions must be provided, the optimized position (p0) and its two neighbors (p1, p2). */ - SmoothingEnergy2D(std::vector> m_cells); + SmoothingEnergy2D(const std::array& pts, const double& M, const Vector3d& L_w); TVector initial_position() const; @@ -63,20 +63,40 @@ class SmoothingEnergy2D : public polysolve::nonlinear::Problem void solution_changed(const TVector& new_x) override {} - void add_mass_and_stiffness_matrix(const double& M, const Vector3d& L_w); + static void + local_mass_and_stiffness(const std::array& pts, double& M, Vector3d& L_w); - static void compute_local_mass_and_stiffness( - std::vector> m_cells, - double& M, - Vector3d& L_w); + static void + uniform_mass_and_stiffness(const std::array& pts, double& M, Vector3d& L_w); private: - MatrixXd m_x; - Vector3d m_L_w; - double m_M; - double m_M_inv; + /** + * The optimized point and its (two) neighbors. For the optimization, we use the convention that + * m_pts[0] is the optimized vertex and the neighbors m_pts[1] and m_pts[2] are fixed. + */ + std::array m_pts; + /** + * The first row of stiffness matrix L_w which belongs to the optimized vertex. By convention, + * this needs to be sorted such that the first entry belongs to the optimized vertex. + * + * The full stiffness matrix L_w has shape NxN (more specifically in our case 3x3 as we have + * three points), but the rows of fixed vertices are empty. Therefore, we only need to store the + * one non-empty row of the optimized vertex. By convention of m_pts, this must be the first + * one, i.e., L_w.row(0). + */ + Vector3d m_L_w_row; + double m_M; // the mass of the optimized vertex + double m_M_inv; // the inverse mass, 1 / M - Vector3d m_LTML_row1; // derived from the other matrices + /** + * We only need the first row of LTML, because all other vertices are fixed. + * + * L = M_inv * L_w + * LTML = L^T * M * L = L_w^T * M_inv * L_w + * We only want to optimize one vertex, so we only need one row of LTML: + * LTML.row(0) = M_inv * L_w(0,0) * L_w.row(0) <-- this is what is stored in m_LTML_row + */ + Vector3d m_LTML_row; }; } // namespace wmtk::optimization \ No newline at end of file diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp index 549f589f24..2814ece5a2 100644 --- a/tests/test_optimization.cpp +++ b/tests/test_optimization.cpp @@ -154,16 +154,14 @@ TEST_CASE("amips_plus_dirichlet_energy_2d", "[energies]") TEST_CASE("smoothing_energy_2d", "[energies]") { - const Vector2d p0(0.8, 1); - std::vector> cells; - cells.push_back({{p0[0], p0[1], 0, 0}}); - cells.push_back({{p0[0], p0[1], 1, 0}}); + const Vector2d p0(0.8, 1); // this vertex is optimized + const Vector2d p1(0, 0); + const Vector2d p2(1, 0); - optimization::SmoothingEnergy2D energy(cells); double M; Vector3d L_w; - optimization::SmoothingEnergy2D::compute_local_mass_and_stiffness(cells, M, L_w); - energy.add_mass_and_stiffness_matrix(M, L_w); + optimization::SmoothingEnergy2D::local_mass_and_stiffness({{p0, p1, p2}}, M, L_w); + optimization::SmoothingEnergy2D energy({{p0, p1, p2}}, M, L_w); CHECK(energy.is_step_valid(p0, p0)); From 44212b2afd673ee6349b3ae53a2472172e6e5fc2 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Mon, 16 Mar 2026 11:03:15 +0100 Subject: [PATCH 06/19] clean up --- .../ImageSimulationMeshTri.cpp | 53 ++++--------------- src/wmtk/optimization/EnvelopeEnergy.cpp | 8 +-- src/wmtk/optimization/EnvelopeEnergy.hpp | 4 +- 3 files changed, 15 insertions(+), 50 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 7b0b54c21a..d94cb03ec9 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -1227,7 +1227,6 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) old_pos.transpose(), VA[vid].m_pos.transpose()); - std::vector> surface_assemble; std::array surface_pts; if (VA[vid].m_is_on_surface) { const auto es = get_order1_edges_for_vertex(vid); @@ -1243,33 +1242,6 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) surface_pts[i + 1] = m_vertex_attribute.at(neighbor_id).m_pos; } - std::set unique_eid; - for (const size_t fid : locs) { - for (size_t j = 0; j < 3; j++) { - Tuple f_t = tuple_from_edge(fid, j); - size_t eid = f_t.eid(*this); - auto [it, suc] = unique_eid.emplace(eid); - if (!suc) { - continue; - } - if (m_edge_attribute[eid].m_is_surface_fs) { - auto vs_id = get_edge_vids(f_t); - if (vs_id[1] == vid) { - std::swap(vs_id[1], vs_id[0]); - } - if (vs_id[0] != vid) { - continue; // does not contain point of interest - } - std::array coords; - for (int k = 0; k < 2; k++) { - for (int kk = 0; kk < 2; kk++) { - coords[k * 2 + kk] = VA[vs_id[k]].m_pos[kk]; - } - } - surface_assemble.emplace_back(coords); - } - } - } // project to surface //{ // assert(m_envelope->initialized()); @@ -1282,10 +1254,9 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) const auto& M = m_surface_mass[vid]; const auto& L_w = m_surface_stiffness[vid]; - // auto smooth_energy = std::make_shared(surface_assemble); auto smooth_energy = std::make_shared(surface_pts, M, L_w); auto envelope_energy = - std::make_shared(m_envelope, surface_assemble); + std::make_shared(m_envelope, surface_pts); auto energy_sum = std::make_shared(); energy_sum->add_energy(amips_energy); energy_sum->add_energy(smooth_energy, 1e2); @@ -1293,6 +1264,7 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) total_energy = energy_sum; } + // solve { auto x = amips_energy->initial_position(); try { @@ -1303,23 +1275,16 @@ bool ImageSimulationMeshTri::smooth_after(const Tuple& t) m_vertex_attribute[vid].m_pos = x; } - for (auto& n : surface_assemble) { - for (int kk = 0; kk < 2; kk++) { - n[kk] = VA[vid].m_pos[kk]; - } - } - // check surface containment - for (const auto& n : surface_assemble) { - std::array edge; - for (int k = 0; k < 2; k++) { - for (int kk = 0; kk < 2; kk++) { - edge[k][kk] = n[k * 2 + kk]; + if (VA[vid].m_is_on_surface) { + for (size_t i = 0; i < 2; ++i) { + std::array edge; + edge[0] = VA[vid].m_pos; + edge[1] = surface_pts[i + 1]; + if (m_envelope->is_outside(edge)) { + return false; } } - if (m_envelope->is_outside(edge)) { - return false; - } } // quality (only check if not on surface) diff --git a/src/wmtk/optimization/EnvelopeEnergy.cpp b/src/wmtk/optimization/EnvelopeEnergy.cpp index 1549b2275d..150bcef7fa 100644 --- a/src/wmtk/optimization/EnvelopeEnergy.cpp +++ b/src/wmtk/optimization/EnvelopeEnergy.cpp @@ -4,9 +4,9 @@ namespace wmtk::optimization { EnvelopeEnergy2D::EnvelopeEnergy2D( const std::shared_ptr& envelope, - const std::vector>& cells) + const const std::array& pts) : m_envelope(envelope) - , m_cells(cells) + , m_pts(pts) { assert(m_envelope); } @@ -40,10 +40,10 @@ bool EnvelopeEnergy2D::is_step_valid(const TVector& x0, const TVector& x1) if (m_envelope->is_outside(r)) { return false; } - // for (const auto& cell : m_cells) { + // for (size_t i = 0; i < 2; ++i) { // std::array edge; // edge[0] = r; - // edge[1] = Vector2d(cell[2], cell[3]); + // edge[1] = m_pts[i + 1]; // if (m_envelope->is_outside(edge)) { // return false; // } diff --git a/src/wmtk/optimization/EnvelopeEnergy.hpp b/src/wmtk/optimization/EnvelopeEnergy.hpp index df8c5c936d..85f75a6c6e 100644 --- a/src/wmtk/optimization/EnvelopeEnergy.hpp +++ b/src/wmtk/optimization/EnvelopeEnergy.hpp @@ -19,7 +19,7 @@ class EnvelopeEnergy2D : public polysolve::nonlinear::Problem */ EnvelopeEnergy2D( const std::shared_ptr& envelope, - const std::vector>& cells); + const std::array& pts); double value(const TVector& x) override; void gradient(const TVector& x, TVector& gradv) override; @@ -35,7 +35,7 @@ class EnvelopeEnergy2D : public polysolve::nonlinear::Problem private: std::shared_ptr m_envelope; - std::vector> m_cells; + std::array m_pts; }; } // namespace wmtk::optimization \ No newline at end of file From 6514918fa8006a35928286637169a50ce591ec1b Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Mon, 16 Mar 2026 14:23:40 +0100 Subject: [PATCH 07/19] Fix amips energy test. --- src/wmtk/optimization/EnvelopeEnergy.cpp | 2 +- tests/test_optimization.cpp | 10 +++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/wmtk/optimization/EnvelopeEnergy.cpp b/src/wmtk/optimization/EnvelopeEnergy.cpp index 150bcef7fa..084e42e989 100644 --- a/src/wmtk/optimization/EnvelopeEnergy.cpp +++ b/src/wmtk/optimization/EnvelopeEnergy.cpp @@ -4,7 +4,7 @@ namespace wmtk::optimization { EnvelopeEnergy2D::EnvelopeEnergy2D( const std::shared_ptr& envelope, - const const std::array& pts) + const std::array& pts) : m_envelope(envelope) , m_pts(pts) { diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp index 2814ece5a2..979b82793b 100644 --- a/tests/test_optimization.cpp +++ b/tests/test_optimization.cpp @@ -33,7 +33,15 @@ TEST_CASE("amips_energy_2d", "[energies]") auto x = energy.initial_position(); const double e_before = energy.value(x); { - auto m_solver = optimization::create_basic_solver(); + 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 m_solver = polysolve::nonlinear::Solver::create( + nonlinear_solver_params, + linear_solver_params, + 1, + opt_logger()); optimization::deactivate_opt_logger(); CHECK_NOTHROW(m_solver->minimize(energy, x)); From dba86ed984dce643486c77755a194fe9909e8e0b Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Mon, 16 Mar 2026 16:02:29 +0100 Subject: [PATCH 08/19] Update to newest gmp-mpfr on Windows. --- cmake/recipes/gmp_mpfr.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/recipes/gmp_mpfr.cmake b/cmake/recipes/gmp_mpfr.cmake index f2df0ba8ff..b21174b3da 100644 --- a/cmake/recipes/gmp_mpfr.cmake +++ b/cmake/recipes/gmp_mpfr.cmake @@ -26,7 +26,7 @@ if(WIN32) include(CPM) CPMAddPackage( NAME gmp_mpfr - URL "https://github.com/CGAL/cgal/releases/download/v5.2.1/CGAL-5.2.1-win64-auxiliary-libraries-gmp-mpfr.zip" + URL "https://github.com/CGAL/cgal/releases/download/v6.1.1/CGAL-6.1.1-win64-auxiliary-libraries-gmp-mpfr.zip" DOWNLOAD_ONLY YES ) From be715c88e710a5b9d43d8b0c75a364d4ed445000 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Tue, 17 Mar 2026 12:35:17 +0100 Subject: [PATCH 09/19] Use vcpkg in Windows CI to get GMP. --- .github/workflows/continuous.yml | 5 +++++ cmake/recipes/gmp.cmake | 4 ++-- components/cmake/add_component.cmake | 2 +- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml index 980eb2c6a4..381394272e 100644 --- a/.github/workflows/continuous.yml +++ b/.github/workflows/continuous.yml @@ -147,6 +147,10 @@ jobs: - name: Install Ninja uses: seanmiddleditch/gha-setup-ninja@master + - name: VCPKG + run: | + vcpkg install gmp:x64-windows + - name: Dependencies run: | choco install -y sccache @@ -176,6 +180,7 @@ jobs: -DCMAKE_POLICY_DEFAULT_CMP0141=NEW ^ -DCMAKE_MSVC_DEBUG_INFORMATION_FORMAT=Embedded ^ -DCMAKE_BUILD_TYPE=${{ matrix.config }} ^ + -DCMAKE_PREFIX_PATH=C:/vcpkg/installed/x64-windows ^ %WMTK_INTEGRATION_TESTS_FLAG% ^ -B build ^ -S . diff --git a/cmake/recipes/gmp.cmake b/cmake/recipes/gmp.cmake index 1bb8cac3c7..c318961abb 100644 --- a/cmake/recipes/gmp.cmake +++ b/cmake/recipes/gmp.cmake @@ -5,11 +5,11 @@ if(TARGET gmp::gmp) return() endif() +message(STATUS "Third-party: creating target 'gmp::gmp'") + # Download precompiled .dll on Windows include(gmp_mpfr) -message(STATUS "Third-party: creating target 'gmp::gmp'") - # Find_package will look for our downloaded lib on Windows, and system-wide on Linux/macOS find_package(GMP REQUIRED) diff --git a/components/cmake/add_component.cmake b/components/cmake/add_component.cmake index cb2bfa5355..9e07273b4d 100644 --- a/components/cmake/add_component.cmake +++ b/components/cmake/add_component.cmake @@ -25,7 +25,7 @@ function(add_component COMPONENT_NAME) "components_map[\"${COMPONENT_NAME}\"] = wmtk::components::${COMPONENT_NAME}::${COMPONENT_NAME};\n") # add component to wmtk app - target_link_libraries(wmtk_app PRIVATE wmtk::${COMPONENT_NAME}) + target_link_libraries(wmtk_app PUBLIC wmtk::${COMPONENT_NAME}) # add component to python bindings if(WMTK_PYBIND) target_link_libraries(pywmtk PUBLIC wmtk::${COMPONENT_NAME}) From dd417f0bc85dc3a55fcdd5291fa649ff2ff2a3e0 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Tue, 17 Mar 2026 14:36:09 +0100 Subject: [PATCH 10/19] missed a few important pieces --- cmake/FindGMP.cmake | 1 + cmake/recipes/gmp_mpfr.cmake | 14 +++++++++++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/cmake/FindGMP.cmake b/cmake/FindGMP.cmake index 6440df7702..79e9394b97 100644 --- a/cmake/FindGMP.cmake +++ b/cmake/FindGMP.cmake @@ -28,6 +28,7 @@ if(WIN32) find_file(GMP_RUNTIME_LIB NAMES gmp.dll + gmp-10.dll libgmp-10.dll PATHS ENV GMP_DIR diff --git a/cmake/recipes/gmp_mpfr.cmake b/cmake/recipes/gmp_mpfr.cmake index b21174b3da..cfa9129458 100644 --- a/cmake/recipes/gmp_mpfr.cmake +++ b/cmake/recipes/gmp_mpfr.cmake @@ -18,7 +18,19 @@ if(WIN32) lib ) - if(GMP_TRY_FIND_INCLUDES AND GMP_TRY_FIND_LIBRARIES) + find_file(GMP_TRY_RUNTIME_LIB + NAMES + gmp.dll + gmp-10.dll + libgmp-10.dll + PATHS + ENV GMP_DIR + ${LIB_INSTALL_DIR} + PATH_SUFFIXES + lib + ) + + if(GMP_TRY_FIND_INCLUDES AND GMP_TRY_FIND_LIBRARIES AND GMP_TRY_RUNTIME_LIB) message(STATUS "System GMP found: ${GMP_TRY_FIND_INCLUDES}") else() message(WARNING "Third-party: downloading gmp + mpfr. The code may be crazy slow! Please install GMP using your preferred package manager. Do not forget to delete your /CMakeCache.txt for the changes to take effect.") From 9734ecfeb13a243d3513cc1d905820592bc2be2c Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Tue, 17 Mar 2026 15:38:31 +0100 Subject: [PATCH 11/19] Update recipes. --- cmake/FindGMP.cmake | 3 +++ cmake/recipes/fenvelope.cmake | 18 +++++++++--------- cmake/recipes/gmp_mpfr.cmake | 7 ++++--- cmake/recipes/nanoflann.cmake | 5 +---- 4 files changed, 17 insertions(+), 16 deletions(-) diff --git a/cmake/FindGMP.cmake b/cmake/FindGMP.cmake index 79e9394b97..a852e60ed7 100644 --- a/cmake/FindGMP.cmake +++ b/cmake/FindGMP.cmake @@ -34,11 +34,14 @@ if(WIN32) ENV GMP_DIR ${LIB_INSTALL_DIR} PATH_SUFFIXES + bin lib ) list(APPEND GMP_EXTRA_VARS GMP_RUNTIME_LIB) endif() +# message(STATUS "GMP Paths: \n ${GMP_INCLUDES}\n ${GMP_LIBRARIES}\n ${GMP_RUNTIME_LIB}") + include(FindPackageHandleStandardArgs) find_package_handle_standard_args(GMP REQUIRED_VARS diff --git a/cmake/recipes/fenvelope.cmake b/cmake/recipes/fenvelope.cmake index e2b5c7b084..45e3a1cfdc 100644 --- a/cmake/recipes/fenvelope.cmake +++ b/cmake/recipes/fenvelope.cmake @@ -10,21 +10,21 @@ set(FAST_ENVELOPE_ENABLE_TBB OFF) include(cli11) -# HACK because there is a linker error on Windows otherwise -if(WIN32) - set(FAST_ENVELOPE_WITH_GEOGRAM_PREDICATES ON) - set(FAST_ENVELOPE_WITH_GEOGRAM_PSM_PREDICATES OFF) -else() - set(FAST_ENVELOPE_WITH_GEOGRAM_PREDICATES OFF) - set(FAST_ENVELOPE_WITH_GEOGRAM_PSM_PREDICATES ON) -endif() +# # HACK because there is a linker error on Windows otherwise +# if(WIN32) +# set(FAST_ENVELOPE_WITH_GEOGRAM_PREDICATES ON) +# set(FAST_ENVELOPE_WITH_GEOGRAM_PSM_PREDICATES OFF) +# else() +# set(FAST_ENVELOPE_WITH_GEOGRAM_PREDICATES OFF) +# set(FAST_ENVELOPE_WITH_GEOGRAM_PSM_PREDICATES ON) +# endif() message(STATUS "Third-party: creating target 'FastEnvelope::FastEnvelope'") include(CPM) -CPMAddPackage("gh:daniel-zint/fast-envelope#a631e3555ff7bd495a66d61a311f5aafe00a2f08") +CPMAddPackage("gh:daniel-zint/fast-envelope#5a826b77cdcf099dc693eac04612af04d32c8c03") set_target_properties(FastEnvelope PROPERTIES FOLDER third_party) add_library(FastEnvelope::FastEnvelope ALIAS FastEnvelope) \ No newline at end of file diff --git a/cmake/recipes/gmp_mpfr.cmake b/cmake/recipes/gmp_mpfr.cmake index cfa9129458..0d194f133a 100644 --- a/cmake/recipes/gmp_mpfr.cmake +++ b/cmake/recipes/gmp_mpfr.cmake @@ -18,7 +18,7 @@ if(WIN32) lib ) - find_file(GMP_TRY_RUNTIME_LIB + find_file(GMP_TRY_FIND_RUNTIME_LIB NAMES gmp.dll gmp-10.dll @@ -27,11 +27,12 @@ if(WIN32) ENV GMP_DIR ${LIB_INSTALL_DIR} PATH_SUFFIXES + bin lib ) - if(GMP_TRY_FIND_INCLUDES AND GMP_TRY_FIND_LIBRARIES AND GMP_TRY_RUNTIME_LIB) - message(STATUS "System GMP found: ${GMP_TRY_FIND_INCLUDES}") + if(GMP_TRY_FIND_INCLUDES AND GMP_TRY_FIND_LIBRARIES AND GMP_TRY_FIND_RUNTIME_LIB) + message(STATUS "System GMP found: \n ${GMP_TRY_FIND_INCLUDES}\n ${GMP_TRY_FIND_LIBRARIES}\n ${GMP_TRY_FIND_RUNTIME_LIB}") else() message(WARNING "Third-party: downloading gmp + mpfr. The code may be crazy slow! Please install GMP using your preferred package manager. Do not forget to delete your /CMakeCache.txt for the changes to take effect.") diff --git a/cmake/recipes/nanoflann.cmake b/cmake/recipes/nanoflann.cmake index d2709428f3..34a35ab7fd 100644 --- a/cmake/recipes/nanoflann.cmake +++ b/cmake/recipes/nanoflann.cmake @@ -15,7 +15,4 @@ CPMAddPackage( "NANOFLANN_USE_SYSTEM_GTEST OFF" "MASTER_PROJECT_HAS_TARGET_UNINSTALL ON" ) -set_target_properties(nanoflann PROPERTIES FOLDER third_party) -if(WIN32) - set_target_properties(nanoflann_uninstall PROPERTIES FOLDER third_party) -endif() \ No newline at end of file +set_target_properties(nanoflann PROPERTIES FOLDER third_party) \ No newline at end of file From 186a297315db8f45d63fc408619d6db9b639bc5d Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Tue, 17 Mar 2026 16:33:53 +0100 Subject: [PATCH 12/19] Trying again without vcpkg + using newer versions of GitHub actions. --- .github/workflows/continuous.yml | 23 +++++++++++------------ cmake/recipes/fenvelope.cmake | 2 +- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml index 381394272e..7ef28d2a76 100644 --- a/.github/workflows/continuous.yml +++ b/.github/workflows/continuous.yml @@ -27,7 +27,7 @@ jobs: name: Linux steps: - name: Checkout repository - uses: actions/checkout@v4.1.6 + uses: actions/checkout@v6 with: fetch-depth: 10 @@ -39,7 +39,7 @@ jobs: - name: Cache Build id: cache-build - uses: actions/cache@v4 + uses: actions/cache@v5 with: path: ${{ env.CACHE_PATH }} key: ${{ runner.os }}-${{ matrix.config }}-cache-${{ github.sha }} @@ -84,7 +84,7 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v4.1.6 + uses: actions/checkout@v6 with: fetch-depth: 10 @@ -95,7 +95,7 @@ jobs: - name: Cache Build id: cache-build - uses: actions/cache@v4 + uses: actions/cache@v5 with: path: ${{ env.CACHE_PATH }} key: ${{ matrix.os }}-${{ matrix.config }}-cache-${{ github.sha }} @@ -140,16 +140,16 @@ jobs: steps: - name: Checkout repository - uses: actions/checkout@v4.1.6 + uses: actions/checkout@v6 with: fetch-depth: 10 - - name: Install Ninja - uses: seanmiddleditch/gha-setup-ninja@master + # - name: Install Ninja + # uses: seanmiddleditch/gha-setup-ninja@master - - name: VCPKG - run: | - vcpkg install gmp:x64-windows + # - name: VCPKG + # run: | + # vcpkg install gmp:x64-windows - name: Dependencies run: | @@ -158,7 +158,7 @@ jobs: - name: Cache build id: cache-build - uses: actions/cache@v4 + uses: actions/cache@v5 with: path: ${{ env.SCCACHE_DIR }} key: ${{ runner.os }}-${{ matrix.config }}-cache-${{ github.sha }} @@ -180,7 +180,6 @@ jobs: -DCMAKE_POLICY_DEFAULT_CMP0141=NEW ^ -DCMAKE_MSVC_DEBUG_INFORMATION_FORMAT=Embedded ^ -DCMAKE_BUILD_TYPE=${{ matrix.config }} ^ - -DCMAKE_PREFIX_PATH=C:/vcpkg/installed/x64-windows ^ %WMTK_INTEGRATION_TESTS_FLAG% ^ -B build ^ -S . diff --git a/cmake/recipes/fenvelope.cmake b/cmake/recipes/fenvelope.cmake index 45e3a1cfdc..94d1229e52 100644 --- a/cmake/recipes/fenvelope.cmake +++ b/cmake/recipes/fenvelope.cmake @@ -24,7 +24,7 @@ message(STATUS "Third-party: creating target 'FastEnvelope::FastEnvelope'") include(CPM) -CPMAddPackage("gh:daniel-zint/fast-envelope#5a826b77cdcf099dc693eac04612af04d32c8c03") +CPMAddPackage("gh:daniel-zint/fast-envelope#0046a216171237aab654dd1580aa8a5eaef675ee") set_target_properties(FastEnvelope PROPERTIES FOLDER third_party) add_library(FastEnvelope::FastEnvelope ALIAS FastEnvelope) \ No newline at end of file From 72fd36c0a38cbd9a9612331f45783112e4cd0856 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 09:33:43 +0100 Subject: [PATCH 13/19] Add more debug messages for copying DLLs. --- app/CMakeLists.txt | 5 ++++- cmake/FindGMP.cmake | 3 ++- cmake/wmtk_copy_dll.cmake | 11 ++++++----- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 239c58957c..c1479518eb 100755 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -20,11 +20,14 @@ endif() # ############################################################################### add_executable(wmtk_app main.cpp) -target_link_libraries(wmtk_app PRIVATE +target_link_libraries(wmtk_app PUBLIC + wmtk::toolkit nlohmann_json::nlohmann_json CLI11::CLI11 ) +wmtk_copy_dll(wmtk_app) + if(WMTK_PYBIND) add_subdirectory(pywmtk) endif() \ No newline at end of file diff --git a/cmake/FindGMP.cmake b/cmake/FindGMP.cmake index a852e60ed7..37d526539f 100644 --- a/cmake/FindGMP.cmake +++ b/cmake/FindGMP.cmake @@ -38,9 +38,10 @@ if(WIN32) lib ) list(APPEND GMP_EXTRA_VARS GMP_RUNTIME_LIB) + + message(STATUS "Windows GMP Paths: \n ${GMP_INCLUDES}\n ${GMP_LIBRARIES}\n ${GMP_RUNTIME_LIB}") endif() -# message(STATUS "GMP Paths: \n ${GMP_INCLUDES}\n ${GMP_LIBRARIES}\n ${GMP_RUNTIME_LIB}") include(FindPackageHandleStandardArgs) find_package_handle_standard_args(GMP diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 9a74845709..8934c84ad9 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -80,11 +80,12 @@ function(wmtk_copy_dll target) # Instruction to copy target file if it exists string(APPEND COPY_SCRIPT_CONTENT - "if(EXISTS \"$\")\n " - "execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " - "\"$\" " - "\"$/$\")\n" - "endif()\n" + "message(STATUS \"For target ${target}, copy $ to $/$\")\n" + # "if(EXISTS \"$\")\n " + "execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " + "\"$\" " + "\"$/$\")\n" + # "endif()\n" ) endforeach() From f4e8334b4272a94ca35c7d95888aace7b60cfda7 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 10:12:35 +0100 Subject: [PATCH 14/19] Trying some stuff with Catch2 test discovery. --- cmake/wmtk_copy_dll.cmake | 25 ++++++++++++++++--------- tests/CMakeLists.txt | 4 ++-- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 8934c84ad9..64201b51b9 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -57,11 +57,12 @@ function(wmtk_copy_dll target) message(FATAL_ERROR "wmtk_copy_dll() was called on a non-executable target: ${target}") endif() - # Create a custom command to do the actual copy. This needs to be executed before Catch2's POST_BUILD command, - # so we define this as a PRE_LINK command for the executable target. + # Create a custom command to do the actual copy. We use POST_BUILD so the executable output directory + # always exists before attempting to copy dependencies. add_custom_command( TARGET ${target} - PRE_LINK + POST_BUILD + COMMAND ${CMAKE_COMMAND} -E make_directory "$" COMMAND ${CMAKE_COMMAND} -E touch "${CMAKE_BINARY_DIR}/runtime_deps/copy_dll_${target}_$.cmake" COMMAND ${CMAKE_COMMAND} -P "${CMAKE_BINARY_DIR}/runtime_deps/copy_dll_${target}_$.cmake" COMMENT "Copying dlls for target ${target}" @@ -80,12 +81,18 @@ function(wmtk_copy_dll target) # Instruction to copy target file if it exists string(APPEND COPY_SCRIPT_CONTENT - "message(STATUS \"For target ${target}, copy $ to $/$\")\n" - # "if(EXISTS \"$\")\n " - "execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " - "\"$\" " - "\"$/$\")\n" - # "endif()\n" + "if(EXISTS \"$\")\n" + " message(STATUS \"For target ${target}, copy $ to $/$\")\n" + " execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " + "\"$\" " # copy source + "\"$/$\" " # copy destination + "RESULT_VARIABLE COPY_RESULT)\n" + " if(NOT COPY_RESULT EQUAL 0)\n" + " message(FATAL_ERROR \"Failed to copy dependency $ for target ${target}.\")\n" + " endif()\n" + "else()\n " + "message(STATUS \"For target ${target}, cannot copy. Dependency $ does not exist.\")\n" + "endif()\n" ) endforeach() diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f744f73408..93df523b5f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,7 +35,7 @@ target_link_libraries(wmtk_tests PUBLIC wmtk_copy_dll(wmtk_tests) # Register unit tests -catch_discover_tests(wmtk_tests) +catch_discover_tests(wmtk_tests DISCOVERY_MODE PRE_TEST) if(WMTK_BUILD_INTEGRATION_TESTS) message(STATUS "Include integration tests.") @@ -54,6 +54,6 @@ if(WMTK_BUILD_INTEGRATION_TESTS) wmtk_copy_dll(wmtk_integration_tests) # Register unit tests - catch_discover_tests(wmtk_integration_tests) + catch_discover_tests(wmtk_integration_tests DISCOVERY_MODE PRE_TEST) endif() \ No newline at end of file From 0d358bc8519c20f5ee2b2caaed69fc7dc6cc94f2 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 10:43:55 +0100 Subject: [PATCH 15/19] Change debug message. --- cmake/wmtk_copy_dll.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 64201b51b9..676a315dac 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -87,7 +87,7 @@ function(wmtk_copy_dll target) "\"$\" " # copy source "\"$/$\" " # copy destination "RESULT_VARIABLE COPY_RESULT)\n" - " if(NOT COPY_RESULT EQUAL 0)\n" + " if(NOT EXISTS \"$/$\")\n" " message(FATAL_ERROR \"Failed to copy dependency $ for target ${target}.\")\n" " endif()\n" "else()\n " From c1184bbaae6ca605103ec47b5fc68b88febb092d Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 11:25:10 +0100 Subject: [PATCH 16/19] minor clean-up --- cmake/wmtk_copy_dll.cmake | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 676a315dac..0ff399d356 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -85,8 +85,7 @@ function(wmtk_copy_dll target) " message(STATUS \"For target ${target}, copy $ to $/$\")\n" " execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " "\"$\" " # copy source - "\"$/$\" " # copy destination - "RESULT_VARIABLE COPY_RESULT)\n" + "\"$/$\")\n" # copy destination " if(NOT EXISTS \"$/$\")\n" " message(FATAL_ERROR \"Failed to copy dependency $ for target ${target}.\")\n" " endif()\n" From 28cae85e52bc443f8a689bcc7d68e1ce5696de64 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 11:34:59 +0100 Subject: [PATCH 17/19] extending the copy log messages --- cmake/wmtk_copy_dll.cmake | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 0ff399d356..59ff9f0f12 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -82,15 +82,25 @@ function(wmtk_copy_dll target) # Instruction to copy target file if it exists string(APPEND COPY_SCRIPT_CONTENT "if(EXISTS \"$\")\n" - " message(STATUS \"For target ${target}, copy $ to $/$\")\n" + " set(DEST_DIR \"$\")\n" + " set(DEST_FILE \"$/$\")\n" + " message(STATUS \"Creating destination directory: \${DEST_DIR}\")\n" + " execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory \"\${DEST_DIR}\" RESULT_VARIABLE DIR_RESULT)\n" + " if(NOT DIR_RESULT EQUAL 0)\n" + " message(FATAL_ERROR \"Failed to create destination directory \${DEST_DIR} for target ${target}. Result: \${DIR_RESULT}\")\n" + " endif()\n" + " message(STATUS \"Copying: $ -> \${DEST_FILE}\")\n" " execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " - "\"$\" " # copy source - "\"$/$\")\n" # copy destination - " if(NOT EXISTS \"$/$\")\n" - " message(FATAL_ERROR \"Failed to copy dependency $ for target ${target}.\")\n" + "\"$\" \"\${DEST_FILE}\" " + "RESULT_VARIABLE COPY_RESULT ERROR_VARIABLE COPY_ERROR)\n" + " if(NOT COPY_RESULT EQUAL 0)\n" + " message(FATAL_ERROR \"Failed to copy $ to \${DEST_FILE} for target ${target}. Result: \${COPY_RESULT}, Error: \${COPY_ERROR}\")\n" + " endif()\n" + " if(NOT EXISTS \"\${DEST_FILE}\")\n" + " message(FATAL_ERROR \"Copy completed (exit 0) but destination file does not exist: \${DEST_FILE}\")\n" " endif()\n" "else()\n " - "message(STATUS \"For target ${target}, cannot copy. Dependency $ does not exist.\")\n" + "message(STATUS \"Dependency does not exist (will not copy): $\")\n" "endif()\n" ) endforeach() From f4af2b94651f051cefbde94d1925e00923c1c740 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 12:07:34 +0100 Subject: [PATCH 18/19] Fix command quoting in wmtk_copy_dll for directory creation and file copying --- cmake/wmtk_copy_dll.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 59ff9f0f12..5c56796d87 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -85,12 +85,12 @@ function(wmtk_copy_dll target) " set(DEST_DIR \"$\")\n" " set(DEST_FILE \"$/$\")\n" " message(STATUS \"Creating destination directory: \${DEST_DIR}\")\n" - " execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory \"\${DEST_DIR}\" RESULT_VARIABLE DIR_RESULT)\n" + " execute_process(COMMAND \"${CMAKE_COMMAND}\" -E make_directory \"\${DEST_DIR}\" RESULT_VARIABLE DIR_RESULT)\n" " if(NOT DIR_RESULT EQUAL 0)\n" " message(FATAL_ERROR \"Failed to create destination directory \${DEST_DIR} for target ${target}. Result: \${DIR_RESULT}\")\n" " endif()\n" " message(STATUS \"Copying: $ -> \${DEST_FILE}\")\n" - " execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different " + " execute_process(COMMAND \"${CMAKE_COMMAND}\" -E copy_if_different " "\"$\" \"\${DEST_FILE}\" " "RESULT_VARIABLE COPY_RESULT ERROR_VARIABLE COPY_ERROR)\n" " if(NOT COPY_RESULT EQUAL 0)\n" From 09ee3e254c067bff27b9608f22e01ad09bb792f8 Mon Sep 17 00:00:00 2001 From: Daniel Zint Date: Wed, 18 Mar 2026 12:54:16 +0100 Subject: [PATCH 19/19] Clean up copy dll and add Windows integration tests for TetWild. --- cmake/wmtk_copy_dll.cmake | 20 +++++++------------- tests/CMakeLists.txt | 4 ++-- tests/integration_tests.cpp | 4 +--- 3 files changed, 10 insertions(+), 18 deletions(-) diff --git a/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 5c56796d87..6a30c3ac93 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -61,8 +61,7 @@ function(wmtk_copy_dll target) # always exists before attempting to copy dependencies. add_custom_command( TARGET ${target} - POST_BUILD - COMMAND ${CMAKE_COMMAND} -E make_directory "$" + PRE_LINK COMMAND ${CMAKE_COMMAND} -E touch "${CMAKE_BINARY_DIR}/runtime_deps/copy_dll_${target}_$.cmake" COMMAND ${CMAKE_COMMAND} -P "${CMAKE_BINARY_DIR}/runtime_deps/copy_dll_${target}_$.cmake" COMMENT "Copying dlls for target ${target}" @@ -81,26 +80,21 @@ function(wmtk_copy_dll target) # Instruction to copy target file if it exists string(APPEND COPY_SCRIPT_CONTENT - "if(EXISTS \"$\")\n" - " set(DEST_DIR \"$\")\n" + "set(SRC_FILE \"$\")\n" + "if(EXISTS \"\${SRC_FILE}\")\n" " set(DEST_FILE \"$/$\")\n" - " message(STATUS \"Creating destination directory: \${DEST_DIR}\")\n" - " execute_process(COMMAND \"${CMAKE_COMMAND}\" -E make_directory \"\${DEST_DIR}\" RESULT_VARIABLE DIR_RESULT)\n" - " if(NOT DIR_RESULT EQUAL 0)\n" - " message(FATAL_ERROR \"Failed to create destination directory \${DEST_DIR} for target ${target}. Result: \${DIR_RESULT}\")\n" - " endif()\n" - " message(STATUS \"Copying: $ -> \${DEST_FILE}\")\n" + " message(STATUS \"Copying: \${SRC_FILE} -> \${DEST_FILE}\")\n" " execute_process(COMMAND \"${CMAKE_COMMAND}\" -E copy_if_different " - "\"$\" \"\${DEST_FILE}\" " + "\"\${SRC_FILE}\" \"\${DEST_FILE}\" " "RESULT_VARIABLE COPY_RESULT ERROR_VARIABLE COPY_ERROR)\n" " if(NOT COPY_RESULT EQUAL 0)\n" - " message(FATAL_ERROR \"Failed to copy $ to \${DEST_FILE} for target ${target}. Result: \${COPY_RESULT}, Error: \${COPY_ERROR}\")\n" + " message(FATAL_ERROR \"Failed to copy \${SRC_FILE} to \${DEST_FILE} for target ${target}. Result: \${COPY_RESULT}, Error: \${COPY_ERROR}\")\n" " endif()\n" " if(NOT EXISTS \"\${DEST_FILE}\")\n" " message(FATAL_ERROR \"Copy completed (exit 0) but destination file does not exist: \${DEST_FILE}\")\n" " endif()\n" "else()\n " - "message(STATUS \"Dependency does not exist (will not copy): $\")\n" + "message(STATUS \"Dependency does not exist (will not copy): \${SRC_FILE}\")\n" "endif()\n" ) endforeach() diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 93df523b5f..f744f73408 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,7 +35,7 @@ target_link_libraries(wmtk_tests PUBLIC wmtk_copy_dll(wmtk_tests) # Register unit tests -catch_discover_tests(wmtk_tests DISCOVERY_MODE PRE_TEST) +catch_discover_tests(wmtk_tests) if(WMTK_BUILD_INTEGRATION_TESTS) message(STATUS "Include integration tests.") @@ -54,6 +54,6 @@ if(WMTK_BUILD_INTEGRATION_TESTS) wmtk_copy_dll(wmtk_integration_tests) # Register unit tests - catch_discover_tests(wmtk_integration_tests DISCOVERY_MODE PRE_TEST) + catch_discover_tests(wmtk_integration_tests) endif() \ No newline at end of file diff --git a/tests/integration_tests.cpp b/tests/integration_tests.cpp index 9201c9ad83..a6fe21369c 100644 --- a/tests/integration_tests.cpp +++ b/tests/integration_tests.cpp @@ -27,10 +27,8 @@ std::vector input_files{ "qslim_octocat.json", "shortest_edge_collapse_octocat.json", "shortest_edge_collapse_sphere_with_env.json", -#ifndef _WIN32 - "tetwild_octocat.json", // segfaults on the windows CI (reason unknown) + "tetwild_octocat.json", "tetwild_sphere.json", -#endif "topological_offsets_bunny2d.json"};