diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml index 980eb2c6a4..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,12 +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: Dependencies run: | @@ -154,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 }} 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/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 6440df7702..37d526539f 100644 --- a/cmake/FindGMP.cmake +++ b/cmake/FindGMP.cmake @@ -28,16 +28,21 @@ if(WIN32) find_file(GMP_RUNTIME_LIB NAMES gmp.dll + gmp-10.dll libgmp-10.dll PATHS ENV GMP_DIR ${LIB_INSTALL_DIR} PATH_SUFFIXES + bin 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() + include(FindPackageHandleStandardArgs) find_package_handle_standard_args(GMP REQUIRED_VARS diff --git a/cmake/recipes/fenvelope.cmake b/cmake/recipes/fenvelope.cmake index e2b5c7b084..94d1229e52 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#0046a216171237aab654dd1580aa8a5eaef675ee") 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.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/cmake/recipes/gmp_mpfr.cmake b/cmake/recipes/gmp_mpfr.cmake index f2df0ba8ff..0d194f133a 100644 --- a/cmake/recipes/gmp_mpfr.cmake +++ b/cmake/recipes/gmp_mpfr.cmake @@ -18,15 +18,28 @@ if(WIN32) lib ) - if(GMP_TRY_FIND_INCLUDES AND GMP_TRY_FIND_LIBRARIES) - message(STATUS "System GMP found: ${GMP_TRY_FIND_INCLUDES}") + find_file(GMP_TRY_FIND_RUNTIME_LIB + NAMES + gmp.dll + gmp-10.dll + libgmp-10.dll + PATHS + ENV GMP_DIR + ${LIB_INSTALL_DIR} + PATH_SUFFIXES + bin + lib + ) + + 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.") 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 ) 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 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/cmake/wmtk_copy_dll.cmake b/cmake/wmtk_copy_dll.cmake index 9a74845709..6a30c3ac93 100644 --- a/cmake/wmtk_copy_dll.cmake +++ b/cmake/wmtk_copy_dll.cmake @@ -57,8 +57,8 @@ 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 @@ -80,10 +80,21 @@ 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" + "set(SRC_FILE \"$\")\n" + "if(EXISTS \"\${SRC_FILE}\")\n" + " set(DEST_FILE \"$/$\")\n" + " message(STATUS \"Copying: \${SRC_FILE} -> \${DEST_FILE}\")\n" + " execute_process(COMMAND \"${CMAKE_COMMAND}\" -E copy_if_different " + "\"\${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 \${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): \${SRC_FILE}\")\n" "endif()\n" ) endforeach() 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}) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 498167c05a..d94cb03ec9 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -11,6 +11,11 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include @@ -766,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 } } @@ -1080,6 +1088,39 @@ bool ImageSimulationMeshTri::swap_edge_after(const Tuple& t) 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; @@ -1164,72 +1205,89 @@ 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 + 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); wmtk::logger().trace( "old pos {} -> new pos {}", old_pos.transpose(), VA[vid].m_pos.transpose()); - std::vector> surface_assemble; + std::array surface_pts; if (VA[vid].m_is_on_surface) { - 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); - } - } + const auto es = get_order1_edges_for_vertex(vid); + if (es.size() != 2) { + return false; // can only smooth vertices with 2 neighbors } - { - assert(m_envelope->initialized()); - Vector2d project; - m_envelope->nearest_point(VA[vid].m_pos, project); - m_vertex_attribute[vid].m_pos = project; + 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; } - for (auto& n : surface_assemble) { - for (int kk = 0; kk < 2; kk++) { - n[kk] = VA[vid].m_pos[kk]; - } + // 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; + //} + + const auto& M = m_surface_mass[vid]; + const auto& L_w = m_surface_stiffness[vid]; + + auto smooth_energy = std::make_shared(surface_pts, M, L_w); + auto envelope_energy = + 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); + energy_sum->add_energy(envelope_energy, 1e2 * M); + total_energy = energy_sum; + } + + // solve + { + 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; } // 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 + // quality (only check if not on surface) auto max_after_quality = 0.; for (const size_t fid : locs) { if (is_inverted(fid)) { @@ -1239,8 +1297,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; @@ -1483,9 +1543,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/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 0a1caf4bc5..31923b05a7 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,10 @@ 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; + 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) , m_envelope_eps(envelope_eps) @@ -110,6 +115,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/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/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..8ea3e29499 --- /dev/null +++ b/src/wmtk/optimization/DirichletEnergy.cpp @@ -0,0 +1,120 @@ +#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(); +} + + +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) +{ + m_LTML_row = m_M_inv * m_L_w_row[0] * m_L_w_row; +} + +SmoothingEnergy2D::TVector SmoothingEnergy2D::initial_position() const +{ + return m_pts[0]; +} + +double SmoothingEnergy2D::value(const TVector& x) +{ + assert(x.size() == 2); + double energy = 0; + for (size_t i = 0; i < 2; ++i) { + 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; + } + + 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(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_row[0]; +} + +void SmoothingEnergy2D::local_mass_and_stiffness( + const std::array& pts, + double& M, + Vector3d& L_w) +{ + const double e1 = (pts[1] - pts[0]).norm(); + const double e2 = (pts[2] - pts[0]).norm(); + + M = 0.5 * (e1 + 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 new file mode 100644 index 0000000000..a93592c52e --- /dev/null +++ b/src/wmtk/optimization/DirichletEnergy.hpp @@ -0,0 +1,102 @@ +#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; +}; + +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 (p0) and its two neighbors (p1, p2). + */ + SmoothingEnergy2D(const std::array& pts, const double& M, const Vector3d& L_w); + + 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 {} + + static void + local_mass_and_stiffness(const std::array& pts, double& M, Vector3d& L_w); + + static void + uniform_mass_and_stiffness(const std::array& pts, double& M, Vector3d& L_w); + +private: + /** + * 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 + + /** + * 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/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/EnvelopeEnergy.cpp b/src/wmtk/optimization/EnvelopeEnergy.cpp new file mode 100644 index 0000000000..084e42e989 --- /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::array& pts) + : m_envelope(envelope) + , m_pts(pts) +{ + 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 (size_t i = 0; i < 2; ++i) { + // std::array edge; + // edge[0] = r; + // edge[1] = m_pts[i + 1]; + // 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..85f75a6c6e --- /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::array& pts); + + 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::array m_pts; +}; + +} // 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/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"}; diff --git a/tests/test_optimization.cpp b/tests/test_optimization.cpp new file mode 100644 index 0000000000..979b82793b --- /dev/null +++ b/tests/test_optimization.cpp @@ -0,0 +1,195 @@ +#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 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)); + } + 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(); +} + +TEST_CASE("smoothing_energy_2d", "[energies]") +{ + const Vector2d p0(0.8, 1); // this vertex is optimized + const Vector2d p1(0, 0); + const Vector2d p2(1, 0); + + double M; + Vector3d 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)); + + { + 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 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