From cc625946786b8e2fa39d0585738b82dc884207b5 Mon Sep 17 00:00:00 2001 From: mattoverby Date: Wed, 8 Oct 2025 14:08:24 -0500 Subject: [PATCH 1/4] update lbfgs and cmake --- CMakeLists.txt | 58 +-- LICENSE.txt | 2 +- cmake/libigl.cmake | 11 + deps/mcloptlib/CMakeLists.txt | 42 -- deps/mcloptlib/LICENSE | 21 - deps/mcloptlib/README.md | 29 -- deps/mcloptlib/cmake/FindEigen3.cmake | 81 ---- deps/mcloptlib/include/MCL/Backtracking.hpp | 150 ------ deps/mcloptlib/include/MCL/LBFGS.hpp | 158 ------- deps/mcloptlib/include/MCL/Minimizer.hpp | 113 ----- deps/mcloptlib/include/MCL/MoreThuente.hpp | 330 ------------- deps/mcloptlib/include/MCL/Newton.hpp | 79 ---- deps/mcloptlib/include/MCL/NonLinearCG.hpp | 86 ---- deps/mcloptlib/include/MCL/Problem.hpp | 135 ------ deps/mcloptlib/include/MCL/TrustRegion.hpp | 208 --------- deps/mcloptlib/include/MCL/WolfeBisection.hpp | 110 ----- deps/mcloptlib/test/TestProblem.hpp | 130 ------ deps/mcloptlib/test/testSolvers.cpp | 210 --------- src/Interp.hpp | 2 +- src/LBFGS.hpp | 440 ++++++++++++++++++ src/MPM.cpp | 1 + src/Particle.hpp | 1 + src/Solver.cpp | 47 +- src/Solver.hpp | 9 +- test/solver.cpp | 9 + 25 files changed, 507 insertions(+), 1955 deletions(-) create mode 100644 cmake/libigl.cmake delete mode 100644 deps/mcloptlib/CMakeLists.txt delete mode 100644 deps/mcloptlib/LICENSE delete mode 100644 deps/mcloptlib/README.md delete mode 100644 deps/mcloptlib/cmake/FindEigen3.cmake delete mode 100644 deps/mcloptlib/include/MCL/Backtracking.hpp delete mode 100644 deps/mcloptlib/include/MCL/LBFGS.hpp delete mode 100644 deps/mcloptlib/include/MCL/Minimizer.hpp delete mode 100644 deps/mcloptlib/include/MCL/MoreThuente.hpp delete mode 100644 deps/mcloptlib/include/MCL/Newton.hpp delete mode 100644 deps/mcloptlib/include/MCL/NonLinearCG.hpp delete mode 100644 deps/mcloptlib/include/MCL/Problem.hpp delete mode 100644 deps/mcloptlib/include/MCL/TrustRegion.hpp delete mode 100644 deps/mcloptlib/include/MCL/WolfeBisection.hpp delete mode 100644 deps/mcloptlib/test/TestProblem.hpp delete mode 100644 deps/mcloptlib/test/testSolvers.cpp create mode 100644 src/LBFGS.hpp create mode 100644 test/solver.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f925de..8fc6554 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,6 +25,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") set(CMAKE_BUILD_TYPE Release) add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") +option(MPM_TESTING_ONLY "Only compile lib and tests" OFF) ############################################################ # @@ -32,25 +33,8 @@ add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") # ############################################################ -# libigl for rendering -include(DownloadProject) -option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" OFF) # off = header only -option(LIBIGL_WITH_COMISO "Use CoMiso" OFF) -option(LIBIGL_WITH_EMBREE "Use Embree" OFF) -option(LIBIGL_WITH_OPENGL "Use OpenGL" ON) -option(LIBIGL_WITH_OPENGL_GLFW "Use GLFW" ON) -option(LIBIGL_WITH_OPENGL_GLFW_IMGUI "Use ImGui" OFF) -option(LIBIGL_WITH_PNG "Use PNG" OFF) -option(LIBIGL_WITH_TETGEN "Use Tetgen" OFF) -option(LIBIGL_WITH_TRIANGLE "Use Triangle" OFF) -option(LIBIGL_WITH_PREDICATES "Use exact predicates" OFF) -option(LIBIGL_WITH_XML "Use XML" OFF) -download_project(PROJ libigl - GIT_REPOSITORY https://github.com/libigl/libigl.git - GIT_TAG main - UPDATE_DISCONNECTED 1 - QUIET) -add_subdirectory(${libigl_SOURCE_DIR} ${libigl_BINARY_DIR}) +# Libigl also includes Eigen, etc +include(libigl) # OpenMP find_package(OpenMP) @@ -60,20 +44,18 @@ if (OPENMP_FOUND) add_definitions(-DOMP_NESTED) endif() -# Include headers -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src) -include_directories(SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/deps/mcloptlib/include) -include_directories(SYSTEM ${glad_SOURCE_DIR}/include) -include_directories(SYSTEM ${libigl_SOURCE_DIR}/include) -include_directories(SYSTEM ${libigl_SOURCE_DIR}/external/eigen) +set(MPM_SRC + ${CMAKE_CURRENT_SOURCE_DIR}/src/LBFGS.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Interp.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Particle.hpp) -set(MPM_SRCS - src/Solver.hpp - src/Solver.cpp - src/MPM.hpp - src/MPM.cpp - src/Interp.hpp - src/Particle.hpp) + +add_library(mpm_optimization ${MPM_SRC}) +target_link_libraries(mpm_optimization PUBLIC igl::core) ############################################################ # @@ -81,5 +63,13 @@ set(MPM_SRCS # ############################################################ -add_executable(sphere src/sphere.cpp ${MPM_SRCS}) -target_link_libraries(sphere glad glfw) +if (NOT MPM_TESTING_ONLY) + igl_include(glfw) + add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/src/sphere.cpp) + target_link_libraries(sphere PUBLIC mpm_optimization igl::glfw) +endif() + +enable_testing() +add_executable(test_solver ${CMAKE_CURRENT_SOURCE_DIR}/test/solver.cpp) +target_link_libraries(test_solver PUBLIC mpm_optimization) +add_test(NAME TestSolver COMMAND test_solver) \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt index 56ef4ca..13de209 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -10,7 +10,7 @@ of conditions and the following disclaimer in the documentation and/or other mat provided with the distribution. THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER diff --git a/cmake/libigl.cmake b/cmake/libigl.cmake new file mode 100644 index 0000000..1aae978 --- /dev/null +++ b/cmake/libigl.cmake @@ -0,0 +1,11 @@ +if(TARGET igl::core) + return() +endif() + +include(FetchContent) +FetchContent_Declare( + libigl + GIT_REPOSITORY https://github.com/libigl/libigl.git + GIT_TAG v2.5.0 +) +FetchContent_MakeAvailable(libigl) \ No newline at end of file diff --git a/deps/mcloptlib/CMakeLists.txt b/deps/mcloptlib/CMakeLists.txt deleted file mode 100644 index d33d831..0000000 --- a/deps/mcloptlib/CMakeLists.txt +++ /dev/null @@ -1,42 +0,0 @@ -# The MIT License (MIT) -# Copyright (c) 2017 Matt Overby -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. - -cmake_minimum_required(VERSION 3.1) -project(mcloptlib C CXX) -set(CMAKE_CXX_STANDARD 11) -set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) -set(CMAKE_BUILD_TYPE Debug) -add_definitions( -DMCL_DEBUG=1 ) - -if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-long-long") -endif() - -find_package(Eigen3 REQUIRED) -include_directories(SYSTEM ${EIGEN3_INCLUDE_DIR}) -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include) - -enable_testing() -add_executable(testSolvers test/testSolvers.cpp) -add_test(testLBFGS testSolvers lbfgs) -add_test(testCG testSolvers cg) -add_test(testNewton testSolvers newton) diff --git a/deps/mcloptlib/LICENSE b/deps/mcloptlib/LICENSE deleted file mode 100644 index c95244f..0000000 --- a/deps/mcloptlib/LICENSE +++ /dev/null @@ -1,21 +0,0 @@ -The MIT License (MIT) - -Copyright (c) 2017 Matt Overby - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. diff --git a/deps/mcloptlib/README.md b/deps/mcloptlib/README.md deleted file mode 100644 index bc89cbb..0000000 --- a/deps/mcloptlib/README.md +++ /dev/null @@ -1,29 +0,0 @@ -# mcloptlib - -By Matt Overby -[http://www.mattoverby.net](http://www.mattoverby.net) - -mcloptlib is a header-only optimization library for C++ using Eigen and is geared towards lower-dimension graphics problems. -Originally a fork of [Patrick Wieschollek's CppOptimizationLibrary](https://github.com/PatWie/CppNumericalSolvers), but has diverged considerably. - -## Contents: - -Optimization algorithms: -- Newton's -- Non-linear conjugate gradient -- L-BFGS -- Trust Region with - - Cauchy Point - - Dog Leg - -Linesearch methods: -- Backtracking (Armijo) -- Backtracking with cubic interpolation -- Bisection -- MoreThuente - -## To-do: - -- Option of std::function for value/gradient instead of Problem class -- Sparse Hessians -- Auto-diff diff --git a/deps/mcloptlib/cmake/FindEigen3.cmake b/deps/mcloptlib/cmake/FindEigen3.cmake deleted file mode 100644 index 9c546a0..0000000 --- a/deps/mcloptlib/cmake/FindEigen3.cmake +++ /dev/null @@ -1,81 +0,0 @@ -# - Try to find Eigen3 lib -# -# This module supports requiring a minimum version, e.g. you can do -# find_package(Eigen3 3.1.2) -# to require version 3.1.2 or newer of Eigen3. -# -# Once done this will define -# -# EIGEN3_FOUND - system has eigen lib with correct version -# EIGEN3_INCLUDE_DIR - the eigen include directory -# EIGEN3_VERSION - eigen version - -# Copyright (c) 2006, 2007 Montel Laurent, -# Copyright (c) 2008, 2009 Gael Guennebaud, -# Copyright (c) 2009 Benoit Jacob -# Redistribution and use is allowed according to the terms of the 2-clause BSD license. - -if(NOT Eigen3_FIND_VERSION) - if(NOT Eigen3_FIND_VERSION_MAJOR) - set(Eigen3_FIND_VERSION_MAJOR 2) - endif(NOT Eigen3_FIND_VERSION_MAJOR) - if(NOT Eigen3_FIND_VERSION_MINOR) - set(Eigen3_FIND_VERSION_MINOR 91) - endif(NOT Eigen3_FIND_VERSION_MINOR) - if(NOT Eigen3_FIND_VERSION_PATCH) - set(Eigen3_FIND_VERSION_PATCH 0) - endif(NOT Eigen3_FIND_VERSION_PATCH) - - set(Eigen3_FIND_VERSION "${Eigen3_FIND_VERSION_MAJOR}.${Eigen3_FIND_VERSION_MINOR}.${Eigen3_FIND_VERSION_PATCH}") -endif(NOT Eigen3_FIND_VERSION) - -macro(_eigen3_check_version) - file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header) - - string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}") - set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen3_major_version_match "${_eigen3_version_header}") - set(EIGEN3_MAJOR_VERSION "${CMAKE_MATCH_1}") - string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen3_minor_version_match "${_eigen3_version_header}") - set(EIGEN3_MINOR_VERSION "${CMAKE_MATCH_1}") - - set(EIGEN3_VERSION ${EIGEN3_WORLD_VERSION}.${EIGEN3_MAJOR_VERSION}.${EIGEN3_MINOR_VERSION}) - if(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - set(EIGEN3_VERSION_OK FALSE) - else(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - set(EIGEN3_VERSION_OK TRUE) - endif(${EIGEN3_VERSION} VERSION_LESS ${Eigen3_FIND_VERSION}) - - if(NOT EIGEN3_VERSION_OK) - - message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in ${EIGEN3_INCLUDE_DIR}, " - "but at least version ${Eigen3_FIND_VERSION} is required") - endif(NOT EIGEN3_VERSION_OK) -endmacro(_eigen3_check_version) - -if (EIGEN3_INCLUDE_DIR) - - # in cache already - _eigen3_check_version() - set(EIGEN3_FOUND ${EIGEN3_VERSION_OK}) - -else (EIGEN3_INCLUDE_DIR) - - find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library - PATHS - ${CMAKE_INSTALL_PREFIX}/include - ${KDE4_INCLUDE_DIR} - PATH_SUFFIXES eigen3 eigen - ) - - if(EIGEN3_INCLUDE_DIR) - _eigen3_check_version() - endif(EIGEN3_INCLUDE_DIR) - - include(FindPackageHandleStandardArgs) - find_package_handle_standard_args(Eigen3 DEFAULT_MSG EIGEN3_INCLUDE_DIR EIGEN3_VERSION_OK) - - mark_as_advanced(EIGEN3_INCLUDE_DIR) - -endif(EIGEN3_INCLUDE_DIR) - diff --git a/deps/mcloptlib/include/MCL/Backtracking.hpp b/deps/mcloptlib/include/MCL/Backtracking.hpp deleted file mode 100644 index 90d0764..0000000 --- a/deps/mcloptlib/include/MCL/Backtracking.hpp +++ /dev/null @@ -1,150 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_BACKTRACKING_H -#define MCL_BACKTRACKING_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - - -// -// Old reliable backtracking Armijo -// -template -class Backtracking { -public: - typedef Eigen::Matrix VecX; - - static inline Scalar search(int verbose, int max_iters, Scalar decrease, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - // First things first, check descent norm - const Scalar t_eps = std::numeric_limits::epsilon(); - if( p.norm() <= t_eps ){ return decrease; } - - const Scalar tau = 0.7; - Scalar alpha = alpha0; - VecX grad; - if( DIM == Eigen::Dynamic ){ grad = VecX::Zero(x.rows()); } - Scalar fx0 = problem.gradient(x, grad); - Scalar gtp = grad.dot(p); - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - Scalar fxa = problem.value(x + alpha*p); - Scalar fx0_fxa = fx0 + alpha*decrease*gtp; // Armijo condition I - if( fxa <= fx0_fxa ){ break; } // sufficient decrease - alpha *= tau; - } - - if( iter >= max_iters ){ - if( verbose > 0 ){ printf("Backtracking::search Error: Reached max_iters\n"); } - return -1; - } - - return alpha; - } - -}; // end class Backtracking - - -// -// Backtracking Armijo with cubic interpolation -// -template -class BacktrackingCubic { -public: - typedef Eigen::Matrix VecX; - - static inline Scalar search(int verbose, int max_iters, Scalar decrease, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - // First things first, check descent norm - const Scalar t_eps = std::numeric_limits::epsilon(); - if( p.norm() <= t_eps ){ return decrease; } - - Scalar alpha = alpha0; - VecX grad; - if( DIM == Eigen::Dynamic ){ grad = VecX::Zero(x.rows()); } - Scalar fx0 = problem.gradient(x, grad); - Scalar gtp = grad.dot(p); - Scalar fxp = fx0; - Scalar alphap = alpha; - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - Scalar fxa = problem.value(x + alpha*p); - Scalar fx0_fxa = fx0 + alpha*decrease*gtp; // Armijo condition I - if( fxa <= fx0_fxa ){ break; } // sufficient decrease - - Scalar alpha_tmp = iter == 0 ? - ( gtp / (2.0 * (fx0 + gtp - fxa)) ) : - cubic( fx0, gtp, fxa, alpha, fxp, alphap ); - fxp = fxa; - alphap = alpha; - alpha = range( alpha_tmp, 0.1*alpha, 0.5*alpha ); - } - - if( iter >= max_iters ){ - if( verbose > 0 ){ printf("BacktrackingCubic::search Error: Reached max_iters\n"); } - return -1; - } - - return alpha; - } - -private: - static inline Scalar range( Scalar alpha, Scalar low, Scalar high ){ - if( alpha < low ){ return low; } - else if( alpha > high ){ return high; } - return alpha; - } - - // Cubic interpolation - // fx0 = f(x0) - // gtp = f'(x0)^T p - // fxa = f(x0 + alpha*p) - // alpha = step length - // fxp = previous fxa - // alphap = previous alpha - static inline Scalar cubic( Scalar fx0, Scalar gtp, Scalar fxa, Scalar alpha, Scalar fxp, Scalar alphap ){ - typedef Eigen::Matrix Vec2; - typedef Eigen::Matrix Mat2; - - Scalar mult = 1.0 / ( alpha*alpha * alphap*alphap * (alpha-alphap) ); - Mat2 A; - A(0,0) = alphap*alphap; A(0,1) = -alpha*alpha; - A(1,0) = -alphap*alphap*alphap; A(1,1) = alpha*alpha*alpha; - Vec2 B; - B[0] = fxa - fx0 - alpha*gtp; B[1] = fxp - fx0 - alphap*gtp; - Vec2 r = mult * A * B; - if( std::abs(r[0]) <= 0.0 ){ return -gtp / (2.0*r[1]); } // if quadratic - Scalar d = std::sqrt( r[1]*r[1] - 3.0*r[0]*gtp ); // discrim - return (-r[1] + d) / (3.0*r[0]); - } - -}; // end class BacktrackingCubic - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/LBFGS.hpp b/deps/mcloptlib/include/MCL/LBFGS.hpp deleted file mode 100644 index bd06032..0000000 --- a/deps/mcloptlib/include/MCL/LBFGS.hpp +++ /dev/null @@ -1,158 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 University of Minnesota -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_LBFGS_H -#define MCL_LBFGS_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -// L-BFGS implementation based on Nocedal & Wright Numerical Optimization book (Section 7.2) -// DIM = dimension of the problem -// M = history window -// -// Original Author: Ioannis Karamouzas -// -template -class LBFGS : public Minimizer { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatM; - typedef Eigen::Matrix VecM; - -public: - bool show_denom_warning; // Print out warning for zero denominators - - LBFGS() : show_denom_warning(false) { - this->m_settings.max_iters = 50; - show_denom_warning = this->m_settings.verbose > 0 ? true : false; - } - - // Returns number of iterations used - int minimize(Problem &problem, VecX &x){ - - MatM s, y; - VecM alpha, rho; - VecX grad, q, grad_old, x_old, x_last; - - if( DIM==Eigen::Dynamic ){ - int dim = x.rows(); - s = MatM::Zero(dim,M); - y = MatM::Zero(dim,M); - alpha = VecM::Zero(M); - rho = VecM::Zero(M); - grad = VecX::Zero(dim); - q = VecX::Zero(dim); - grad_old = VecX::Zero(dim); - x_old = VecX::Zero(dim); - x_last = VecX::Zero(dim); - } - - problem.gradient(x, grad); - - Scalar gamma_k = 1.0; - Scalar alpha_init = 1.0; - - int global_iter = 0; - int max_iters = this->m_settings.max_iters; - int verbose = this->m_settings.verbose; - - for( int k=0; k= 0; --i){ - rho(i) = 1.0 / ((s.col(i)).dot(y.col(i))); - alpha(i) = rho(i)*(s.col(i)).dot(q); - q = q - alpha(i)*y.col(i); - } - - // L-BFGS second - loop recursion - q = gamma_k*q; - for(int i = 0; i < iter; ++i){ - Scalar beta = rho(i)*q.dot(y.col(i)); - q = q + (alpha(i) - beta)*s.col(i); - } - - // is there a descent - Scalar dir = q.dot(grad); - if(dir <= 0 ){ - q = grad; - max_iters -= k; - k = 0; - alpha_init = std::min(1.0, 1.0 / grad.template lpNorm() ); - } - - Scalar rate = this->linesearch(x, -q, problem, alpha_init); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("LBFGS::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x -= rate * q; - if( problem.converged(x_last,x,grad) ){ break; } - - problem.gradient(x,grad); - VecX s_temp = x - x_old; - VecX y_temp = grad - grad_old; - - // update the history - if(k < M){ - s.col(k) = s_temp; - y.col(k) = y_temp; - } - else { - s.leftCols(M - 1) = s.rightCols(M - 1).eval(); - s.rightCols(1) = s_temp; - y.leftCols(M - 1) = y.rightCols(M - 1).eval(); - y.rightCols(1) = y_temp; - } - - Scalar denom = y_temp.dot(y_temp); - if( std::abs(denom) <= 0 ){ - if( show_denom_warning ){ - printf("LBFGS::minimize Warning: Encountered a zero denominator\n"); - } - break; - } - gamma_k = s_temp.dot(y_temp) / denom; - alpha_init = 1.0; - - } - - return global_iter; - - } // end minimize -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/Minimizer.hpp b/deps/mcloptlib/include/MCL/Minimizer.hpp deleted file mode 100644 index 70aa7f8..0000000 --- a/deps/mcloptlib/include/MCL/Minimizer.hpp +++ /dev/null @@ -1,113 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_MINIMIZER_H -#define MCL_MINIMIZER_H - -#include "Problem.hpp" -#include "Backtracking.hpp" -#include "MoreThuente.hpp" -#include "WolfeBisection.hpp" -#include - -namespace mcl { -namespace optlib { - -// The different line search methods currently implemented -enum class LSMethod { - None = 0, // use step length = 1, not recommended - MoreThuente, // TODO test this one for correctness - Backtracking, // basic backtracking with sufficient decrease - BacktrackingCubic, // backtracking with cubic interpolation - WeakWolfeBisection // slow -}; - -// Trust region subproblem method (see TrustRegion.hpp) -enum class TRMethod { - CauchyPoint, - DogLeg -}; - -// -// Base class for optimization algs -// -template -class Minimizer { -public: - typedef Eigen::Matrix VecX; - static const int FAILURE = -1; // returned by minimize if an error is encountered - - struct Settings { - int verbose; // higher = more printouts - int max_iters; // usually changed by derived constructors - int ls_max_iters; // max line search iters - Scalar ls_decrease; // sufficient decrease param - LSMethod ls_method; // see LSMethod (above) - TRMethod tr_method; // see TRMethod (above) - - Settings() : verbose(0), max_iters(100), - ls_max_iters(100000), ls_decrease(1e-4), - ls_method(LSMethod::BacktrackingCubic), - tr_method(TRMethod::DogLeg) - {} - } m_settings; - - // - // Performs optimization - // - virtual int minimize(Problem &problem, VecX &x) = 0; - - -protected: - - // Line search method/options can be changed through m_settings. - Scalar linesearch(const VecX &x, const VecX &p, Problem &prob, double alpha0) const { - double alpha = alpha0; - int mi = m_settings.ls_max_iters; - int v = m_settings.verbose; - Scalar sd = m_settings.ls_decrease; - switch( m_settings.ls_method ){ - default:{ - alpha = Backtracking::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::None: { alpha = 1.0; } break; - case LSMethod::MoreThuente: { - alpha = MoreThuente::search(x, p, prob, alpha0); - } break; - case LSMethod::Backtracking: { - alpha = Backtracking::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::BacktrackingCubic: { - alpha = BacktrackingCubic::search(v, mi, sd, x, p, prob, alpha0); - } break; - case LSMethod::WeakWolfeBisection: { - alpha = WolfeBisection::search(v, mi, x, p, prob, alpha0); - } break; - } - return alpha; - } // end do linesearch - -}; // class minimizer - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/MoreThuente.hpp b/deps/mcloptlib/include/MCL/MoreThuente.hpp deleted file mode 100644 index dd3ff15..0000000 --- a/deps/mcloptlib/include/MCL/MoreThuente.hpp +++ /dev/null @@ -1,330 +0,0 @@ -// The MIT License (MIT) -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. -// -// From https://github.com/PatWie/CppNumericalSolvers -// - -#ifndef MCL_MORETHUENTE_H -#define MCL_MORETHUENTE_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - -template -class MoreThuente { -private: - typedef Eigen::Matrix VectorX; - -public: - - static inline Scalar search(const VectorX &x, const VectorX &p, Problem &problem, Scalar alpha0){ - Scalar alpha = alpha0; - cvsrch(problem, x, alpha, p); - return alpha; - } - - static void cvsrch(Problem &problem, const VectorX &x0, Scalar &stp, const VectorX &s) { - int info = 0; - int infoc = 1; - const Scalar xtol = 1e-15; - const Scalar ftol = 1e-4; - const Scalar gtol = 1e-2; - const Scalar stpmin = 1e-15; - const Scalar stpmax = 1e15; - const Scalar xtrapf = 4; - const int maxfev = 20; - int nfev = 0; - int dim = x0.rows(); - - VectorX g; - if( DIM == Eigen::Dynamic ){ g = VectorX::Zero(dim); } - else{ g.setZero(); } - - Scalar f = problem.gradient(x0, g); - Scalar dginit = g.dot(s); - if (dginit >= 0.0) { - // no descent direction - return; - } - - bool brackt = false; - bool stage1 = true; - - Scalar finit = f; - Scalar dgtest = ftol * dginit; - Scalar width = stpmax - stpmin; - Scalar width1 = 2 * width; - VectorX x = x0.eval(); - - Scalar stx = 0.0; - Scalar fx = finit; - Scalar dgx = dginit; - Scalar sty = 0.0; - Scalar fy = finit; - Scalar dgy = dginit; - - Scalar stmin = 0.0; - Scalar stmax = 0.0; - - const int max_iters = 100000; - int iter = 0; - for( ; iter= stmax))) - || (nfev >= maxfev - 1 ) || (infoc == 0) - || (brackt & (stmax - stmin <= xtol * stmax))) { - stp = stx; - } - - // test new point - x = x0 + stp * s; - f = problem.gradient(x, g); - nfev++; - Scalar dg = g.dot(s); - Scalar ftest1 = finit + stp * dgtest; - - // all possible convergence tests - if ((brackt & ((stp <= stmin) | (stp >= stmax))) | (infoc == 0)) - info = 6; - - if ((stp == stpmax) & (f <= ftest1) & (dg <= dgtest)) - info = 5; - - if ((stp == stpmin) & ((f > ftest1) | (dg >= dgtest))) - info = 4; - - if (nfev >= maxfev) - info = 3; - - if (brackt & (stmax - stmin <= xtol * stmax)) - info = 2; - - if ((f <= ftest1) & (fabs(dg) <= gtol * (-dginit))) - info = 1; - - // terminate when convergence reached - if (info != 0) - return; - - if (stage1 & (f <= ftest1) & (dg >= std::min(ftol, gtol)*dginit)) - stage1 = false; - - if (stage1 & (f <= fx) & (f > ftest1)) { - Scalar fm = f - stp * dgtest; - Scalar fxm = fx - stx * dgtest; - Scalar fym = fy - sty * dgtest; - Scalar dgm = dg - dgtest; - Scalar dgxm = dgx - dgtest; - Scalar dgym = dgy - dgtest; - - cstep( stx, fxm, dgxm, sty, fym, dgym, stp, fm, dgm, brackt, stmin, stmax, infoc); - - fx = fxm + stx * dgtest; - fy = fym + sty * dgtest; - dgx = dgxm + dgtest; - dgy = dgym + dgtest; - } else { - // this is ugly and some variables should be moved to the class scope - cstep( stx, fx, dgx, sty, fy, dgy, stp, f, dg, brackt, stmin, stmax, infoc); - } - - if (brackt) { - if (fabs(sty - stx) >= 0.66 * width1) - stp = stx + 0.5 * (sty - stx); - - width1 = width; - width = fabs(sty - stx); - } - - } // end while true - - if( iter == max_iters ){ - throw std::runtime_error("MoreThuente::linesearch Error: Reached max_iter"); - } - - return; - } - - static void cstep(Scalar& stx, Scalar& fx, Scalar& dx, Scalar& sty, Scalar& fy, Scalar& dy, Scalar& stp, - Scalar& fp, Scalar& dp, bool& brackt, Scalar& stpmin, Scalar& stpmax, int& info) { - info = 0; - bool bound = false; - - // Check the input parameters for errors. - if ((brackt & ((stp <= std::min(stx, sty) ) | (stp >= std::max(stx, sty)))) | (dx * (stp - stx) >= 0.0) - | (stpmax < stpmin)) { - return; - } - - Scalar sgnd = dp * (dx / fabs(dx)); - - Scalar stpf = 0; - Scalar stpc = 0; - Scalar stpq = 0; - - if (fp > fx) { - info = 1; - bound = true; - Scalar theta = 3. * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max(dx, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dx / s) * (dp / s)); - if (stp < stx) - gamma = -gamma; - - Scalar p = (gamma - dx) + theta; - Scalar q = ((gamma - dx) + gamma) + dp; - Scalar r = p / q; - stpc = stx + r * (stp - stx); - stpq = stx + ((dx / ((fx - fp) / (stp - stx) + dx)) / 2.) * (stp - stx); - if (fabs(stpc - stx) < fabs(stpq - stx)) - stpf = stpc; - else - stpf = stpc + (stpq - stpc) / 2; - - brackt = true; - } else if (sgnd < 0.0) { - info = 2; - bound = false; - Scalar theta = 3 * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max(dx, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dx / s) * (dp / s)); - if (stp > stx) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = ((gamma - dp) + gamma) + dx; - Scalar r = p / q; - stpc = stp + r * (stx - stp); - stpq = stp + (dp / (dp - dx)) * (stx - stp); - if (fabs(stpc - stp) > fabs(stpq - stp)) - stpf = stpc; - else - stpf = stpq; - - brackt = true; - } else if (fabs(dp) < fabs(dx)) { - info = 3; - bound = 1; - Scalar theta = 3 * (fx - fp) / (stp - stx) + dx + dp; - Scalar s = std::max(theta, std::max( dx, dp)); - Scalar gamma = s * sqrt(std::max(static_cast(0.), (theta / s) * (theta / s) - (dx / s) * (dp / s))); - if (stp > stx) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = (gamma + (dx - dp)) + gamma; - Scalar r = p / q; - if ((r < 0.0) & (gamma != 0.0)) { - stpc = stp + r * (stx - stp); - } else if (stp > stx) { - stpc = stpmax; - } else { - stpc = stpmin; - } - stpq = stp + (dp / (dp - dx)) * (stx - stp); - if (brackt) { - if (fabs(stp - stpc) < fabs(stp - stpq)) { - stpf = stpc; - } else { - stpf = stpq; - } - } else { - if (fabs(stp - stpc) > fabs(stp - stpq)) { - stpf = stpc; - } else { - stpf = stpq; - } - } - } else { - info = 4; - bound = false; - if (brackt) { - Scalar theta = 3 * (fp - fy) / (sty - stp) + dy + dp; - Scalar s = std::max(theta, std::max(dy, dp)); - Scalar gamma = s * sqrt((theta / s) * (theta / s) - (dy / s) * (dp / s)); - if (stp > sty) - gamma = -gamma; - - Scalar p = (gamma - dp) + theta; - Scalar q = ((gamma - dp) + gamma) + dy; - Scalar r = p / q; - stpc = stp + r * (sty - stp); - stpf = stpc; - } else if (stp > stx) - stpf = stpmax; - else { - stpf = stpmin; - } - } - - if (fp > fx) { - sty = stp; - fy = fp; - dy = dp; - } else { - if (sgnd < 0.0) { - sty = stx; - fy = fx; - dy = dx; - } - stx = stp; - fx = fp; - dx = dp; - } - - stpf = std::min(stpmax, stpf); - stpf = std::max(stpmin, stpf); - stp = stpf; - - if (brackt & bound) { - if (sty > stx) { - stp = std::min(stx + static_cast(0.66) * (sty - stx), stp); - } else { - stp = std::max(stx + static_cast(0.66) * (sty - stx), stp); - } - } - - return; - - } // end cstep - -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/Newton.hpp b/deps/mcloptlib/include/MCL/Newton.hpp deleted file mode 100644 index f1632b8..0000000 --- a/deps/mcloptlib/include/MCL/Newton.hpp +++ /dev/null @@ -1,79 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_NEWTON_H -#define MCL_NEWTON_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class Newton : public Minimizer { -private: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - -public: - Newton() { - this->m_settings.max_iters = 20; - } - - int minimize(Problem &problem, VectorX &x){ - - VectorX grad, delta_x, x_last; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.resize(dim); - grad.resize(dim); - delta_x.resize(dim); - } - - int verbose = this->m_settings.verbose; - int max_iters = this->m_settings.max_iters; - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - problem.gradient(x,grad); - problem.solve_hessian(x,grad,delta_x); - - Scalar rate = this->linesearch(x, delta_x, problem, 1.0); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("Newton::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x += rate * delta_x; - if( problem.converged(x_last,x,grad) ){ break; } - } - - return iter; - } - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/NonLinearCG.hpp b/deps/mcloptlib/include/MCL/NonLinearCG.hpp deleted file mode 100644 index 064fd99..0000000 --- a/deps/mcloptlib/include/MCL/NonLinearCG.hpp +++ /dev/null @@ -1,86 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_NONLINEARCG_H -#define MCL_NONLINEARCG_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class NonLinearCG : public Minimizer { -private: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - -public: - NonLinearCG() { - this->m_settings.max_iters = 100; - } - - int minimize(Problem &problem, VectorX &x){ - - VectorX grad, grad_old, p, x_last; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.setZero(dim); - grad.resize(dim); - grad_old.resize(dim); - p.resize(dim); - } - - int verbose = this->m_settings.verbose; - int max_iters = this->m_settings.max_iters; - int iter=0; - for( ; iterlinesearch(x, p, problem, 1.0); - - if( rate <= 0 ){ - if( verbose > 0 ){ printf("NonLinearCG::minimize: Failure in linesearch\n"); } - return Minimizer::FAILURE; - } - - x_last = x; - x += rate*p; - grad_old = grad; - - if( problem.converged(x_last,x,grad) ){ break; } - } - return iter; - } // end minimize - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/Problem.hpp b/deps/mcloptlib/include/MCL/Problem.hpp deleted file mode 100644 index 54ebc3c..0000000 --- a/deps/mcloptlib/include/MCL/Problem.hpp +++ /dev/null @@ -1,135 +0,0 @@ -// The MIT License (MIT) -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_PROBLEM_H -#define MCL_PROBLEM_H - -#if MCL_DEBUG == 1 -#include -#endif - -#include - -namespace mcl { -namespace optlib { - -template -class Problem { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - -public: - // Returns true if the solver has converged - // x0 is the result of the previous iteration - // x1 is the result at the current iteration - // grad is the gradient at the last iteration - virtual bool converged(const VecX &x0, const VecX &x1, const VecX &grad) = 0; - - // Compute just the value - virtual Scalar value(const VecX &x) = 0; - - // Compute the objective value and the gradient - virtual Scalar gradient(const VecX &x, VecX &grad){ - finiteGradient(x, grad); - return value(x); - } - - // Compute hessian - virtual void hessian(const VecX &x, MatX &hessian){ - finiteHessian(x, hessian); - } - - // Solve dx = H^-1 -g (used by Newton's) - virtual void solve_hessian(const VecX &x, const VecX &grad, VecX &dx){ - MatX hess; - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - hess = MatX::Zero(dim,dim); - } - hessian(x,hess); // hessian at x_n - - // Going with with high-accurate, low requirements as default factorization for lin-solve - // Copied from https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html - // Method Requirements Spd (sm) Spd (lg) Accuracy - // partialPivLu() Invertible ++ ++ + - // fullPivLu() None - - - +++ - // householderQr() None ++ ++ + - // colPivHouseholderQr() None ++ - +++ - // fullPivHouseholderQr() None - - - +++ - // llt() PD +++ +++ + - // ldlt() P/N SD +++ + ++ - if( DIM == Eigen::Dynamic || DIM > 4 ){ - dx = hess.fullPivLu().solve(-grad); - } - else { - dx = -hess.inverse()*grad; - } - } - - // Gradient with finite differences - inline void finiteGradient(const VecX &x, VecX &grad){ - const int accuracy = 0; // accuracy can be 0, 1, 2, 3 - const Scalar eps = 2.2204e-6; - const std::vector< std::vector > coeff = - { {1, -1}, {1, -8, 8, -1}, {-1, 9, -45, 45, -9, 1}, {3, -32, 168, -672, 672, -168, 32, -3} }; - const std::vector< std::vector > coeff2 = - { {1, -1}, {-2, -1, 1, 2}, {-3, -2, -1, 1, 2, 3}, {-4, -3, -2, -1, 1, 2, 3, 4} }; - const std::vector dd = {2, 12, 60, 840}; - int dim = x.rows(); - if( grad.rows() != dim ){ grad = VecX::Zero(dim); } - else{ grad.setZero(); } - for(int d = 0; d < dim; ++d){ - for (int s = 0; s < 2*(accuracy+1); ++s){ - VecX xx = x.eval(); - xx[d] += coeff2[accuracy][s]*eps; - grad[d] += coeff[accuracy][s]*value(xx); - } - grad[d] /= (dd[accuracy]* eps); - } - } // end finite grad - - // Hessian with finite differences - inline void finiteHessian(const VecX &x, MatX &hess){ - const Scalar eps = std::numeric_limits::epsilon()*10e7; - int dim = x.rows(); - if( hess.rows() != dim || hess.cols() != dim ){ hess = MatX::Zero(dim,dim); } - for(int i = 0; i < dim; ++i){ - for(int j = 0; j < dim; ++j){ - VecX xx = x; - Scalar f4 = value(xx); - xx[i] += eps; - xx[j] += eps; - Scalar f1 = value(xx); - xx[j] -= eps; - Scalar f2 = value(xx); - xx[j] += eps; - xx[i] -= eps; - Scalar f3 = value(xx); - hess(i, j) = (f1 - f2 - f3 + f4) / (eps * eps); - } - } - } // end finite hess -}; - -} -} - -#endif diff --git a/deps/mcloptlib/include/MCL/TrustRegion.hpp b/deps/mcloptlib/include/MCL/TrustRegion.hpp deleted file mode 100644 index 52bc961..0000000 --- a/deps/mcloptlib/include/MCL/TrustRegion.hpp +++ /dev/null @@ -1,208 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_TRUSTREGION_H -#define MCL_TRUSTREGION_H - -#include "Minimizer.hpp" - -namespace mcl { -namespace optlib { - -template -class TrustRegion : public Minimizer { -private: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - -public: - // - // The trust region method operates on an approximate hessian (B). - // To keep the interface simple, Problem::hessian is used to obtain - // this approximation. It shouldn't be much of a problem, since you - // can still use approximations for Newton's (the only other method - // implemented that requires a hessian evaluation). - // - // TODO Replace B with possibly sparse representation - // - TrustRegion() { - this->m_settings.max_iters = 100; - } - - int minimize(Problem &problem, VecX &x){ - - Scalar delta_k = 2.0; // trust region radius - const Scalar delta_max = 8.0; // max trust region radius - const Scalar eta = 0.125; // min reduction ratio allowed (0m_settings.tr_method; - int max_iters = this->m_settings.max_iters; - int verbose = this->m_settings.verbose; - - VecX grad, dx, x_last; - MatX B; // Approximate hessian - if( DIM == Eigen::Dynamic ){ - int dim = x.rows(); - x_last.resize(dim); // last variable - grad.resize(dim); // grad at k - dx.resize(dim); // descent direction - B.resize(dim,dim); - } - - // Init gradient and hessian - Scalar fxk = problem.gradient(x,grad); // gradient and objective - problem.hessian(x,B); // get hessian (or approximation) - problem.solve_hessian(x,grad,dx); // attempt with newtons - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - // If it's outside the trust region, pick new descent - if( dx.norm() > delta_k ){ - eval_subproblem(m, delta_k, grad, B, dx); - } - - Scalar fxdx = problem.value(x+dx); - - // Compute reduction ratio - Scalar rho_k = eval_reduction(fxk, fxdx, dx, grad, B); - Scalar dx_norm = dx.norm(); - - // Update trust region radius - if( rho_k < 0.25 ){ delta_k = 0.25*delta_k; } - - // Full step, good approximation - else if( rho_k > 0.75 && std::abs(dx_norm-delta_k) <= 0.0 ){ - delta_k = std::min( 2.0*delta_k, delta_max ); - } - - // Take a step, otherwise need to re-eval sub problem - if( rho_k > eta ){ - - x_last = x; - x = x + dx; - - // Only need to compute gradient and hessian - // if x has actually changed. - fxk = problem.gradient(x,grad); // gradient and objective - if( problem.converged(x_last,x,grad) ){ break; } - - // I think I should improve this as to not call both hessian - // and solve_hessian, which likely causes redundant computation. - problem.hessian(x,B); // get hessian (or approximation) - problem.solve_hessian(x,grad,dx); // attempt with newtons - } - - if( std::isnan(rho_k) ){ - if( verbose ){ printf("\n**TrustRegion Error: NaN reduction"); } - return Minimizer::FAILURE; - } - - } - - return iter; - - } // end minimize - -protected: - - // Assumes coeffs size 3 - static inline Scalar max_roots(Scalar a, Scalar b, Scalar c){ - - Scalar d = (b*b) - (4.0*a*c); - if( d > 0 ){ - Scalar sqrt_d = std::sqrt(d); - Scalar r1 = (-b + sqrt_d) / (2.0*a); - Scalar r2 = (-b - sqrt_d) / (2.0*a); - return std::max(r1,r2); - } - // Should I do something with the real/imaginary parts? - // Scalar real_part = -b/(2.0*a); - // Scalar imaginary = std::sqrt(-d)/(2.0*a); - throw std::runtime_error("TrustRegion Error: Problem in quadratic roots"); - return 0; - } - - // fxk = objective at x_k - // fxdx = objective at x_k + dx - // dx = descent direction - // grad_k = gradient at x_k - // B_k = hessian guess. - static inline Scalar eval_reduction( Scalar fxk, Scalar fxdx, const VecX &dx, - const VecX &grad_k, const MatX &B_k ){ - // rho = ( f(x) - f(x-dx) ) / ( model(0) - model(dx) ) - // with model = f(x) + dx^T grad + 0.5 dx^T B dx - Scalar num = fxk - fxdx; - Scalar denom = fxk - ( fxk + dx.dot(grad_k) + 0.5 * dx.dot( B_k * dx ) ); - return num/denom; - } - - static inline void eval_subproblem( - const TRMethod &m, Scalar delta_k, - const VecX &grad, const MatX &B, - VecX &dx ){ - - Scalar gTBg = grad.dot(B*grad); - - switch( m ){ - - // Generally requires a large number of outer solver iterations - case TRMethod::CauchyPoint: { - - Scalar grad_norm = grad.norm(); - Scalar tau = 1.0; - if( gTBg > 0.0 ){ tau = std::min( 1.0, std::pow(grad_norm,3) / (delta_k*gTBg) ); } - dx = ( -tau * delta_k / grad_norm ) * grad; - - } break; - - // Uses steepest descent if possible - case TRMethod::DogLeg: { - - Scalar gTg = grad.dot(grad); - VecX dx_U = ( -gTg / gTBg ) * grad; - Scalar dx_U_norm = dx_U.norm(); - - // Use steepest descent - if( dx_U_norm >= delta_k ){ dx = delta_k/dx_U_norm * dx_U; } - else{ - // Compute tau and update descent - VecX dx_C = dx - dx_U; - Scalar dx_C_norm = dx_C.norm(); - Scalar tau = max_roots( // Ax^2 + Bx + c - dx_C_norm*dx_C_norm, - 2.0*dx_C.dot(dx_U), - dx_U_norm*dx_U_norm - delta_k*delta_k - ); - dx = dx_U + tau*dx_C; - } - - } break; - - } // end swithc method - - } // end eval sub problem - -}; - -} // ns optlib -} // ns mcl - -#endif diff --git a/deps/mcloptlib/include/MCL/WolfeBisection.hpp b/deps/mcloptlib/include/MCL/WolfeBisection.hpp deleted file mode 100644 index b9e6a06..0000000 --- a/deps/mcloptlib/include/MCL/WolfeBisection.hpp +++ /dev/null @@ -1,110 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2018 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#ifndef MCL_WOLFEBISECTION_H -#define MCL_WOLFEBISECTION_H - -#include "Problem.hpp" - -namespace mcl { -namespace optlib { - -// Bisection method for Weak Wolfe conditions -template -class WolfeBisection { -private: - // Strong wolfe conditions: the armijo rule and a stronger curvature condition - // alpha = step length - // fx_ap = f(x + alpha p) - // fx = f(x) - // pT_gx = p^T ( grad f(x) ) - // pT_gx_ap = p^T ( grad f(x + alpha p) ) - static inline bool strong_wolfe( Scalar alpha, - Scalar fx, Scalar pT_gx, - Scalar fx_ap, Scalar pT_gx_ap, - Scalar wolfe_c1, Scalar wolfe_c2 ){ - if( !(fx_ap <= fx + wolfe_c1 * alpha * pT_gx) ){ return false; } // armijo rule - return std::abs( pT_gx_ap ) <= wolfe_c2 * std::abs( pT_gx ); - } - -public: - typedef Eigen::Matrix VecX; - typedef Eigen::Matrix MatX; - - static inline Scalar search(int verbose, int max_iters, const VecX &x, const VecX &p, Problem &problem, Scalar alpha0) { - - const Scalar t_eps = std::numeric_limits::epsilon(); - const Scalar wolfe_c1 = 0.0001; - const Scalar wolfe_c2 = 0.8; // should be 0.1 for CG! - const int dim = x.rows(); - double alpha = alpha0; - double alpha_min = 1e-8; - double alpha_max = 1; - - VecX grad0, grad_new; - if( DIM == Eigen::Dynamic ){ - grad0 = VecX::Zero(dim); - grad_new = VecX::Zero(dim); - } - Scalar fx0 = problem.gradient(x, grad0); - const Scalar gtp = grad0.dot(p); - bool min_set = false; - - int iter = 0; - for( ; iter < max_iters; ++iter ){ - - // Should we stop iterating? - if( std::abs(alpha_max-alpha_min) <= t_eps ){ break; } - - // Step halfway - alpha = ( alpha_max + alpha_min ) * 0.5; - grad_new.setZero(); - Scalar fx_ap = problem.gradient(x + alpha*p, grad_new); - Scalar gt_ap = grad_new.dot( p ); - - // Check the wolfe conditions - bool happy_wolfe = strong_wolfe( alpha, fx0, gtp, fx_ap, gt_ap, wolfe_c1, wolfe_c2 ); - if( happy_wolfe ){ - alpha_min = alpha; - min_set = true; - } - else { alpha_max = alpha; } - - } // end bs iters - - if( iter == max_iters ){ - if( verbose > 0 ){ printf("WolfeBisection::linesearch Error: Reached max_iters\n"); } - return -1; - } - if( !min_set ){ - if( verbose > 0 ){ printf("WolfeBisection::linesearch Error: LS blocked\n"); } - return -1; - } - - return alpha_min; - - } -}; - -} -} - -#endif diff --git a/deps/mcloptlib/test/TestProblem.hpp b/deps/mcloptlib/test/TestProblem.hpp deleted file mode 100644 index c4c8a4f..0000000 --- a/deps/mcloptlib/test/TestProblem.hpp +++ /dev/null @@ -1,130 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#include "MCL/Problem.hpp" - -// min |Ax-b| -class DynProblem : public mcl::optlib::Problem { -public: - typedef Eigen::Matrix VectorX; - typedef Eigen::Matrix MatrixX; - - MatrixX A; - VectorX b; - DynProblem( int dim_ ){ - - // Test on random SPD - A = MatrixX::Random(dim_,dim_); - A = A.transpose() * A; - A = A + MatrixX::Identity(dim_,dim_); - - b = VectorX::Random(dim_); - } - - int dim() const { return b.rows(); } - bool converged(const VectorX &x0, const VectorX &x1, const VectorX &grad){ - - // Check sizes of input - int m_dim = dim(); - if( x0.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: x0 wrong dimension"); - } - if( x1.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: x1 wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::converged: gradient wrong dimension"); - } - - return grad.norm() < 1e-10 || (x0-x1).norm() < 1e-10; - } - - double value(const VectorX &x){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::value: x wrong dimension"); - } - - return (A*x-b).norm(); - } - - double gradient(const VectorX &x, VectorX &grad){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::gradient: x wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::gradient: gradient wrong dimension"); - } - - grad = A*x-b; return value(x); - } - void hessian(const VectorX &x, MatrixX &hess){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::hessian: x wrong dimension"); - } - if( hess.rows() != m_dim || hess.cols() != m_dim ){ - throw std::runtime_error("Error in Problem::hessian: hessian wrong dimension"); - } - hess = A; - } - - void solve_hessian(const VectorX &x, const VectorX &grad, VectorX &dx){ - - // Check sizes of input - int m_dim = dim(); - if( x.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: x wrong dimension"); - } - if( dx.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: dx wrong dimension"); - } - if( grad.rows() != m_dim ){ - throw std::runtime_error("Error in Problem::solve_hessian: gradient wrong dimension"); - } - - // Check to make sure base class function works as expected - Problem::solve_hessian(x,grad,dx); - } -}; - -class Rosenbrock : public mcl::optlib::Problem { -public: - typedef Eigen::Matrix VectorX; - bool converged(const VectorX &x0, const VectorX &x1, const VectorX &grad){ - (void)(x1); (void)(x0); - return grad.norm() < 1e-10; - } - double value(const VectorX &x){ - double a = 1.0 - x[0]; - double b = x[1] - x[0]*x[0]; - return a*a + b*b*100.0; - } - // Test finite diff as well I guess -}; - diff --git a/deps/mcloptlib/test/testSolvers.cpp b/deps/mcloptlib/test/testSolvers.cpp deleted file mode 100644 index 4f047bf..0000000 --- a/deps/mcloptlib/test/testSolvers.cpp +++ /dev/null @@ -1,210 +0,0 @@ -// The MIT License (MIT) -// Copyright (c) 2017 Matt Overby -// -// Permission is hereby granted, free of charge, to any person obtaining a copy -// of this software and associated documentation files (the "Software"), to deal -// in the Software without restriction, including without limitation the rights -// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the Software is -// furnished to do so, subject to the following conditions: -// -// The above copyright notice and this permission notice shall be included in all -// copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -// SOFTWARE. - -#include -#include "TestProblem.hpp" -#include "MCL/LBFGS.hpp" -#include "MCL/NonLinearCG.hpp" -#include "MCL/Newton.hpp" -#include "MCL/TrustRegion.hpp" -#include - -using namespace mcl::optlib; -typedef std::shared_ptr< Minimizer > MinPtr2; // rb -typedef std::shared_ptr< Minimizer > MinPtrD; // linear - -bool test_linear( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest linear:" << std::endl; - typedef Eigen::Matrix VecX; - bool success = true; - - // Using multiple dimensions in the linear case, which changes - // how Newton's performs a solve. - std::vector test_dims = { 4, 16, 64 }; - - int n_test_dims = test_dims.size(); - for( int i=0; im_settings.max_iters = 1; } - else { solvers[i]->m_settings.max_iters = 100; } - solvers[i]->m_settings.verbose = 1; - VecX x = VecX::Zero(dim); - solvers[i]->minimize( cp, x ); - for( int i=0; i 1e-3 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: |Ax-b| = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Linear (" << dim << "): Success" << std::endl; } - else{ success = false; } - - } // end loop solvers - - } // end loop test dims - - return success; - -} - - -bool test_rb( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest Rosenbrock:" << std::endl; - Rosenbrock rb; // Dim = 2, also tests finite gradient/hessian - bool success = true; - - int n_solvers = solvers.size(); - for( int i=0; im_settings.max_iters = 1000; - solvers[i]->m_settings.verbose = 1; - Eigen::Vector2d x = Eigen::Vector2d::Zero(); - solvers[i]->minimize( rb, x ); - - for( int i=0; i<2; ++i ){ - if( std::isnan(x[i]) || std::isinf(x[i]) ){ - std::cerr << "(" << names[i] << ") Bad values in x: " << x[i] << std::endl; - curr_success = false; - } - } - - double rn = (Eigen::Vector2d(1,1) - x).norm(); - if( rn > 1e-8 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: Rosenbrock = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Rosenbrock: Success" << std::endl; } - else{ success = false; } - } - - return success; -} - -// Test what happens when the energy is already minimized -bool test_zero( std::vector &solvers, std::vector &names ){ - - std::cout << "\nTest zero energy:" << std::endl; - typedef Eigen::Matrix VecX; - bool success = true; - int dim = 2; - - // High/Low dimensions, linear, and Eigen::Dynamic. - // The solvers should work, since the Vecs/Mats are resized at run time. - DynProblem cp(dim); - - int n_solvers = solvers.size(); - for( int i=0; im_settings.max_iters = 1; - solvers[i]->m_settings.verbose = 1; - VecX x = cp.A.inverse() * cp.b; - solvers[i]->minimize( cp, x ); - - for( int i=0; i 1e-10 ){ - std::cerr << "(" << names[i] << ") Failed to minimize: |Ax-b| = " << rn << std::endl; - curr_success = false; - } - - if( curr_success ){ std::cout << "(" << names[i] << ") Linear (" << dim << "): Success" << std::endl; } - else{ success = false; } - - } // end loop solvers - - return success; -} - - -int main(int argc, char *argv[] ){ - srand(100); - std::vector< std::string > names; - std::vector< MinPtr2 > min2; - std::vector< MinPtrD > minD; - - std::string mode = "all"; - if( argc == 2 ){ mode = std::string(argv[1]); } - - if( mode=="lbfgs" || mode=="all" ){ - min2.emplace_back( std::make_shared< LBFGS >( LBFGS() ) ); - minD.emplace_back( std::make_shared< LBFGS >( LBFGS() ) ); - names.emplace_back( "lbfgs" ); - } - if( mode=="cg" || mode=="all" ){ - min2.emplace_back( std::make_shared< NonLinearCG >( NonLinearCG() ) ); - minD.emplace_back( std::make_shared< NonLinearCG >( NonLinearCG() ) ); - names.emplace_back( "cg" ); - } - if( mode=="newton" || mode=="all" ){ - min2.emplace_back( std::make_shared< Newton >( Newton() ) ); - minD.emplace_back( std::make_shared< Newton >( Newton() ) ); - names.emplace_back( "newton" ); - } - if( mode=="trustregion" || mode=="all" ){ - min2.emplace_back( std::make_shared< TrustRegion >( TrustRegion() ) ); - minD.emplace_back( std::make_shared< TrustRegion >( TrustRegion() ) ); - names.emplace_back( "trustregion" ); - } - - bool success = true; - success &= test_linear( minD, names ); - success &= test_rb( min2, names ); - success &= test_zero( minD, names ); - if( success ){ - std::cout << "\nSUCCESS!" << std::endl; - return EXIT_SUCCESS; - } - else{ std::cout << "\n**FAILURE!" << std::endl; } - return EXIT_FAILURE; -} - - - diff --git a/src/Interp.hpp b/src/Interp.hpp index 4269055..9520d60 100644 --- a/src/Interp.hpp +++ b/src/Interp.hpp @@ -20,7 +20,7 @@ #ifndef INTERP_HPP #define INTERP_HPP 1 -#include +#include #include namespace mpm diff --git a/src/LBFGS.hpp b/src/LBFGS.hpp new file mode 100644 index 0000000..839ebb3 --- /dev/null +++ b/src/LBFGS.hpp @@ -0,0 +1,440 @@ +// Copyright Matt Overby 2021. +// Distributed under the MIT License. +// From: https://github.com/mattoverby/mclgeom + +#ifndef MCL_LBFGS_HPP +#define MCL_LBFGS_HPP 1 + +#include +#include +#include + +namespace mcl +{ + +// L-BFGS implementation based on Nocedal & Wright Numerical Optimization book (Section 7.2) +// adapted from code by Ioannis Karamouzas. +// +// Function pointers are meant to be used with lambdas, e.g. +// LBFGS lbfgs; +// lbfgs.gradient = [&](const MatrixXd &x, MatrixXd &g)->Scalar { return ... }; +// double obj = lbfgs.minimize(x); +// +// If the g arg is not sized, don't compute gradient. This is to avoid +// extra work when only the objective is needed (line search) and to +// avoid redundant code. Example: +// lbfgs.gradient = [&](const MatrixXd &x, MatrixXd &g)->Scalar +// { +// objective = ... +// if (g.rows() == x.rows()) +// { +// g = ... (reuse computation for objective) +// } +// return objective; +// }; +// +template +class LBFGS +{ +public: + typedef typename MatrixType::Scalar Scalar; + + struct Options + { + int min_iters; + int max_iters; + Scalar abs_tol; // absolute tol if converged(...) not set + Scalar rel_tol; // relative tol if converged(...) not set + int M; // history window size + Scalar gamma; // init Hessian = gamma * I + Options() : + min_iters(0), + max_iters(100), + abs_tol(1e-5), + rel_tol(0), + M(6), + gamma(1) + {} + } options; + + LBFGS(); + + // Output from the last call to minimize(...) + int iters() const { return num_iters; } + Scalar gamma() const { return gamma_k; } + + // Resizes buffers and sets to zero. + // Called during minimize(...) ONLY if there is a change in dof or M. + // Otherwise it's assumed you're picking up where you left off + // on the previous call to minimize(...) + void reset(int rows, int cols); + + // Required: + // computes objective value and gradient + // obj = gradient(x, g) + // If the g arg is not sized, don't compute gradient + std::function gradient; + + // Optional: + // Returns true if the solver should exit, default uses ||g|| converged; + + // Optional: + // Linesearch function, default uses bracketing weak wolfe (slow!) + // obj_k1 = (x, grad, descent, alpha) + // Returns new objective value and updates both x AND gradient + std::function linesearch; + + // Optional: + // Filter descent direction, p = B(p) + // Otherwise p = gamma_k * p is used. + std::function filter; + + // Calls initialize(x) once and iterate(x) until converged + Scalar minimize(MatrixType& x); + + // Initialize the solver + void initialize(MatrixType& x); + + // Take an iteration + // Returns objective + Scalar iterate(MatrixType& x); + + // i.e. bisection with weak Wolfe conditions + Scalar bracketing_weakwolfe( + MatrixType& x, + MatrixType& grad, + const MatrixType& p, + Scalar &alpha) const; + + // Used if converged not set + // Returns true if: + // grad.norm <= abs_tol + // or + // grad.norm() <= rel_tol * x.norm() + bool default_converged( + Scalar curr_obj, + const MatrixType& xprev, + const MatrixType& x, + const MatrixType& grad) const; + + Scalar inner( + const MatrixType &a, + const MatrixType &b) const; + +protected: + bool initialized; + int num_iters, max_iters, k; + Scalar gamma_k, obj_0, obj_k; + std::vector s; + std::vector y; + Eigen::Matrix alpha; + Eigen::Matrix rho; + MatrixType grad; + MatrixType q; + MatrixType descent; + MatrixType grad_old; + MatrixType x_old; + MatrixType x_last; + MatrixType s_temp; + MatrixType y_temp; + +}; // end class LBFGS + +// +// Implementation +// + +template +LBFGS::LBFGS() : + initialized(false), + num_iters(0), + max_iters(0), + k(0), + gamma_k(1), + obj_0(0), + obj_k(0) + {} + +template +void LBFGS::reset(int rows, int cols) +{ + using namespace Eigen; + num_iters = 0; + max_iters = 0; + k = 0; + gamma_k = 1; + obj_0 = std::numeric_limits::max(); + obj_k = std::numeric_limits::max(); + int M = options.M; + s = std::vector(M, MatrixType::Zero(rows,cols)); + y = std::vector(M, MatrixType::Zero(rows,cols)); + alpha = VectorXd::Zero(M); + rho = VectorXd::Zero(M); + grad = MatrixType::Zero(rows, cols); + q = MatrixType::Zero(rows, cols); // inv descent + descent = q; + grad_old = MatrixType::Zero(rows, cols); + x_old = MatrixType::Zero(rows, cols); + x_last = MatrixType::Zero(rows, cols); + s_temp = MatrixType::Zero(rows, cols); + y_temp = MatrixType::Zero(rows, cols); +} + +// Returns number of iterations used +template +typename LBFGS::Scalar +LBFGS::minimize(MatrixType &x) +{ + initialize(x); + + // Did we start at the initializer? + if (num_iters >= options.min_iters && + converged(obj_k,x_last,x,grad)) + { + num_iters = 1; + return obj_k; + } + + for (; k +void LBFGS::initialize(MatrixType& x) +{ + initialized = true; + + if (gradient == nullptr) { + throw std::runtime_error("no gradient function"); + } + + if (converged == nullptr) + { + using namespace std::placeholders; + converged = std::bind(&LBFGS::default_converged, this, _1, _2, _3, _4); + } + + if (linesearch == nullptr) + { + using namespace std::placeholders; + linesearch = std::bind(&LBFGS::bracketing_weakwolfe, this, _1, _2, _3, _4); + } + + // Resize/resize variables? + if (alpha.rows() != options.M || + grad.rows() != x.rows() || + grad.cols() != x.cols()) { + reset(x.rows(), x.cols()); + } + + num_iters = 0; + max_iters = std::max(options.min_iters, options.max_iters); + gamma_k = options.gamma; + obj_0 = gradient(x, grad); + obj_k = obj_0; + k = 0; +} + +template +typename LBFGS::Scalar +LBFGS::iterate(MatrixType& x) +{ + if (!initialized) { + throw std::runtime_error("not initialized"); + } + + x_old = x; + grad_old = grad; + q = grad; + num_iters++; + + // + // Two-loop recursion + // + { + // L-BFGS first - loop recursion + int iter = std::min(options.M, k); + for(int i = iter - 1; i >= 0; --i) + { + Scalar denom = inner(s[i], y[i]); + if (std::abs(denom) <= 0.0) + { + rho(i) = 0; + alpha(i) = 0; + continue; + } + rho(i) = 1.0 / denom; + alpha(i) = rho(i)*inner(s[i], q); + q -= alpha(i) * y[i]; + } + + if (filter != nullptr) { filter(q); } + else { q = gamma_k*q; } + + // L-BFGS second - loop recursion + for(int i = 0; i < iter; ++i) + { + Scalar beta = rho(i)*inner(q, y[i]); + q += (alpha(i) - beta)*s[i]; + } + } + + // + // Perform step + // + { + // If our hess approx is bad and we start going in + // the wrong direction, restart memory + Scalar step_size = 1.0; + Scalar dir = inner(q, grad); + if (dir <= 0) + { + q = grad; + max_iters -= k; // Restart memory + k = 0; + step_size = std::min(1.0, 1.0 / grad.template lpNorm() ); + } + + // We've hit local minima, we have to exit + if (q.squaredNorm() <= 0.0) { + return obj_k; + } + + descent = -q; + x_last = x; + obj_k = linesearch(x, grad, descent, step_size); + if (num_iters >= options.min_iters && converged(obj_k,x_last,x,grad)) { + return obj_k; + } + } + + // + // Correction term + // + { + s_temp = x - x_old; + y_temp = grad - grad_old; + + // update the history + if (k < options.M) + { + s[k] = s_temp; + y[k] = y_temp; + } + else + { + for (int i=0; i 0.0) + { + gamma_k = inner(s_temp, y_temp) / denom; + } + } + + return obj_k; +} + +template +typename LBFGS::Scalar +LBFGS::bracketing_weakwolfe( + MatrixType& x, + MatrixType& g, + const MatrixType &drt, + Scalar &step) const +{ + using namespace Eigen; + + Scalar c1 = 10e-4; + Scalar c2 = 0.9; + if (step <= 0.0) { step = 1.0; } + int x_cols = x.cols(); + + g = MatrixType::Zero(x.rows(), x_cols); + MatrixType xp = x; + const Scalar fx_init = gradient(x,g); + Scalar fx = fx_init; + + Scalar dg_init = inner(g, drt); + if (dg_init > 0) { + throw std::runtime_error("direction increases objective"); + } + + const Scalar test_decr = c1 * dg_init; + Scalar lower = 0; + Scalar upper = std::numeric_limits::infinity(); + int maxiter = 200; + int iter = 0; + for (iter = 0; iter < maxiter; iter++) + { + x = xp + step * drt; + fx = gradient(x, g); + if (fx > fx_init + step * test_decr) { // Armijo rule + upper = step; + } + else + { + Scalar dg = inner(g, drt); + if(dg < c2 * dg_init){ // Weak wolfe + lower = step; + } else { + break; // both met + } + } + step = std::isinf(upper) ? 2*step : lower/2 + upper/2; + } + return fx; + +} // end linesearch + +template +bool LBFGS::default_converged( + Scalar curr_obj, + const MatrixType& xprev, + const MatrixType& x, + const MatrixType& g) const +{ + (void)(curr_obj); + (void)(xprev); + Scalar gnorm = g.norm(); + if (gnorm <= options.abs_tol) { + return true; + } + Scalar xnorm = x.norm(); + if (gnorm < options.rel_tol * xnorm) { + return true; + } + return false; +} + +template +typename LBFGS::Scalar +LBFGS::inner( + const MatrixType &a, + const MatrixType &b) const +{ + int cols = std::min(a.cols(), b.cols()); + Scalar dot = 0; + for (int i=0; i namespace mpm { diff --git a/src/Particle.hpp b/src/Particle.hpp index 1dd9693..6be8521 100644 --- a/src/Particle.hpp +++ b/src/Particle.hpp @@ -22,6 +22,7 @@ #include #include +#include #include "Interp.hpp" namespace mpm diff --git a/src/Solver.cpp b/src/Solver.cpp index 847fb31..13be7ad 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -22,30 +22,6 @@ namespace mpm { -// The system as represented by an optimization problem for generic solvers -class Objective : public mcl::optlib::Problem -{ -public: - Objective(Solver *solver_) : solver(solver_) {} - Solver *solver; - - double value(const Eigen::VectorXd &v) - { - Eigen::VectorXd grad; - return gradient(v,grad); - } - - bool converged(const Eigen::VectorXd &x0, const Eigen::VectorXd &x1, const Eigen::VectorXd &grad) - { - if (grad.norm() < 1e-2) { return true; } - if ((x0-x1).norm() < 1e-2) { return true; } - return false; - } - - double gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad); - -}; // end class Objective - bool Solver::initialize() { // Fill a sphere with particles @@ -78,8 +54,6 @@ bool Solver::initialize() std::cout << "Num Particles: " << m_particles.size() << std::endl; std::cout << "Num Nodes: " << m_grid.size() << std::endl; - // Set up the solver - optimizer.m_settings.ls_method = mcl::optlib::LSMethod::Backtracking; return true; } @@ -169,7 +143,7 @@ void Solver::explicit_solve() } // end compute forces -double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) +double Solver::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) { using namespace Eigen; double tot_energy = 0.0; @@ -178,8 +152,8 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) // Compute Particle Deformation Gradients for new grid velocities // #pragma omp parallel for reduction (+:tot_energy) - for(int i=0; im_particles.size(); ++i){ - Particle *p = solver->m_particles[i]; + for(int i=0; iget_deform_grad(); Matrix3d vel_grad; vel_grad.fill(0.f); @@ -190,7 +164,7 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) MPM::make_grid_loop(p, grid_nodes, wip, dwip); for(int j=0; jm_grid[ grid_nodes[j] ]->active_idx; + int idx = m_grid[ grid_nodes[j] ]->active_idx; Eigen::Vector3d curr_v(v[idx*3+0], v[idx*3+1], v[idx*3+2]); vel_grad += (MPM::timestep_s*curr_v) * dwip[j].transpose(); } @@ -209,9 +183,9 @@ double Objective::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) // Compute energy gradient // #pragma omp parallel for reduction (+:tot_energy) - for(int i=0; iactive_grid.size(); ++i) + for(int i=0; iactive_grid[i]; + GridNode *node = active_grid[i]; Eigen::Vector3d curr_v(v[i*3+0], v[i*3+1], v[i*3+2]); Eigen::Vector3d momentum_grad = node->m * (curr_v - node->v); @@ -249,8 +223,13 @@ void Solver::implicit_solve() } // Minimize - Objective obj(this); - optimizer.minimize(obj, v); + optimizer.gradient = [&](const VectorXd &x, VectorXd &g)->double + { + double objective = gradient(x, g); + return objective; + }; + + optimizer.minimize(v); // Copy solver results back to grid #pragma omp parallel for diff --git a/src/Solver.hpp b/src/Solver.hpp index 127cf4c..53ebdbc 100644 --- a/src/Solver.hpp +++ b/src/Solver.hpp @@ -21,7 +21,7 @@ #define MPM_SOLVER_HPP 1 #include "MPM.hpp" -#include "MCL/LBFGS.hpp" +#include "LBFGS.hpp" namespace mpm { @@ -41,7 +41,7 @@ class Solver std::vector active_grid; // resized each time step std::vector m_particles; void get_vertices(Eigen::MatrixXd &X) const; - mcl::optlib::LBFGS optimizer; + mcl::LBFGS optimizer; // // Settings @@ -61,6 +61,9 @@ class Solver // Returns true on success. bool step(float screen_dt); + // Gradient of the objective. + double gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad); + private: double compute_timestep(double screen_dt); @@ -72,7 +75,7 @@ class Solver void explicit_solve(); void implicit_solve(); -}; // end class system +}; // end class solver } // end namespace mpm diff --git a/test/solver.cpp b/test/solver.cpp new file mode 100644 index 0000000..ddffc6f --- /dev/null +++ b/test/solver.cpp @@ -0,0 +1,9 @@ +#include + +// placeholder +int main(int argc, char *argv[]) +{ + (void)(argc); + (void)(argv); + return EXIT_SUCCESS; +} \ No newline at end of file From fa862dbb319e6cabaf5c4878f1a3fb1c80de2189 Mon Sep 17 00:00:00 2001 From: mattoverby Date: Wed, 8 Oct 2025 14:19:27 -0500 Subject: [PATCH 2/4] mv sphere to test --- .github/workflows/build_and_test.yml | 30 ++++++++++++++++ CMakeLists.txt | 51 +++++++--------------------- test/solver.cpp | 3 ++ {src => test}/sphere.cpp | 0 4 files changed, 45 insertions(+), 39 deletions(-) create mode 100644 .github/workflows/build_and_test.yml rename {src => test}/sphere.cpp (100%) diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml new file mode 100644 index 0000000..2d619d7 --- /dev/null +++ b/.github/workflows/build_and_test.yml @@ -0,0 +1,30 @@ +name: Build and Test + +on: + push: + branches: [ master ] + pull_request: + branches: [ master ] + +jobs: + build_and_test: + runs-on: ubuntu-latest + + steps: + - name: clone + uses: actions/checkout@v4 # Action to check out your repository + + - name: dependencies + run: | + sudo apt-get update + sudo apt-get install -y build-essential cmake + + - name: build + working-directory: ${{ github.workspace }} + run: | + cmake -DMPM_TESTING_ONLY=ON . + make -j + + - name: test + working-directory: ${{ github.workspace }} + run: make test \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 8fc6554..179767f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,42 +1,18 @@ -# Copyright (c) 2016 Matt Overby -# -# MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) -# Redistribution and use in source and binary forms, with or without modification, are -# permitted provided that the following conditions are met: -# 1. Redistributions of source code must retain the above copyright notice, this list of -# conditions and the following disclaimer. -# 2. Redistributions in binary form must reproduce the above copyright notice, this list -# of conditions and the following disclaimer in the documentation and/or other materials -# provided with the distribution. -# THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE -# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -# OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER -# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +cmake_minimum_required(VERSION 3.27) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -cmake_minimum_required(VERSION 3.0) project(mpm C CXX) +set(CMAKE_CXX_STANDARD 17) +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") -set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") -set(CMAKE_BUILD_TYPE Release) - -add_definitions(-DMPM_SRC_DIR="${CMAKE_CURRENT_SOURCE_DIR}") +# Options option(MPM_TESTING_ONLY "Only compile lib and tests" OFF) -############################################################ -# -# Libraries -# -############################################################ - # Libigl also includes Eigen, etc include(libigl) # OpenMP +# (TO-DO: replace with TBB) find_package(OpenMP) if (OPENMP_FOUND) set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") @@ -44,6 +20,7 @@ if (OPENMP_FOUND) add_definitions(-DOMP_NESTED) endif() +# Simulation library set(MPM_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/LBFGS.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Solver.hpp @@ -52,24 +29,20 @@ set(MPM_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/MPM.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interp.hpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Particle.hpp) - - add_library(mpm_optimization ${MPM_SRC}) target_link_libraries(mpm_optimization PUBLIC igl::core) -############################################################ -# -# Binaries -# -############################################################ - +# Interactive viewer if (NOT MPM_TESTING_ONLY) igl_include(glfw) - add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/src/sphere.cpp) + add_executable(sphere ${CMAKE_CURRENT_SOURCE_DIR}/test/sphere.cpp) target_link_libraries(sphere PUBLIC mpm_optimization igl::glfw) + target_include_directories(sphere PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src) endif() +# Unit testing (TO-DO) enable_testing() add_executable(test_solver ${CMAKE_CURRENT_SOURCE_DIR}/test/solver.cpp) target_link_libraries(test_solver PUBLIC mpm_optimization) +target_include_directories(test_solver PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src) add_test(NAME TestSolver COMMAND test_solver) \ No newline at end of file diff --git a/test/solver.cpp b/test/solver.cpp index ddffc6f..aa06e89 100644 --- a/test/solver.cpp +++ b/test/solver.cpp @@ -1,8 +1,11 @@ #include +#include "Solver.hpp" // placeholder int main(int argc, char *argv[]) { + mpm::Solver solver; + (void)(argc); (void)(argv); return EXIT_SUCCESS; diff --git a/src/sphere.cpp b/test/sphere.cpp similarity index 100% rename from src/sphere.cpp rename to test/sphere.cpp From 9ec65f8c408b03355c15fbbecd574ffcef648661 Mon Sep 17 00:00:00 2001 From: mattoverby Date: Thu, 11 Dec 2025 12:25:38 -0600 Subject: [PATCH 3/4] update yml --- .github/workflows/build_and_test.yml | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 2d619d7..9079757 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -17,7 +17,7 @@ jobs: - name: dependencies run: | sudo apt-get update - sudo apt-get install -y build-essential cmake + sudo apt-get install -y build-essential cmake libtbb-dev - name: build working-directory: ${{ github.workspace }} diff --git a/README.md b/README.md index fd36516..1a6cfbf 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ MPM is great for continuum substances like fluids and soft bodies, as well as fo It has received a lot of attention in computer graphics, spurred mostly by the [snow simulation](https://doi.org/10.1145/2461912.2461948) paper. Some time around 2016 I coded up MPM to test out some research ideas, and wrote a little text to go along with it. -My ideas didn't work out, but I put this code up on github instead of letting it collect metaphorical dust on my hard drive. +My ideas didn't work out, but I put this code up on github instead of letting it collect dust on my hard drive. It is not efficient nor bug free. The demo is nothing more than a big neo-Hookean blob dropped on a ground plane. From 99cd8137f1e1589f4325a41d727dc11abe5e7eea Mon Sep 17 00:00:00 2001 From: mattoverby Date: Sat, 13 Dec 2025 10:53:49 -0600 Subject: [PATCH 4/4] minor updates --- .clang-format | 69 +++ .github/workflows/build_and_test.yml | 4 +- CMakeLists.txt | 11 +- README.md | 20 +- cmake/DownloadProject.CMakeLists.cmake.in | 17 - cmake/DownloadProject.cmake | 182 ------ src/Interp.hpp | 134 +++-- src/LBFGS.hpp | 701 +++++++++++----------- src/MPM.cpp | 538 +++++++++-------- src/MPM.hpp | 158 +++-- src/Particle.hpp | 229 +++---- src/Solver.cpp | 490 ++++++++------- src/Solver.hpp | 80 +-- test/solver.cpp | 10 +- test/sphere.cpp | 83 ++- 15 files changed, 1365 insertions(+), 1361 deletions(-) create mode 100644 .clang-format delete mode 100644 cmake/DownloadProject.CMakeLists.cmake.in delete mode 100644 cmake/DownloadProject.cmake diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..b5e8e0a --- /dev/null +++ b/.clang-format @@ -0,0 +1,69 @@ +BasedOnStyle: LLVM +AccessModifierOffset: '-2' +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: 'true' +AlignConsecutiveDeclarations: 'true' +AlignOperands: 'true' +AlignTrailingComments: 'true' +AllowAllParametersOfDeclarationOnNextLine: true +AllowAllArgumentsOnNextLine: true +AllowShortBlocksOnASingleLine: 'false' +AllowShortCaseLabelsOnASingleLine: 'false' +AllowShortFunctionsOnASingleLine: Inline +AllowShortIfStatementsOnASingleLine: 'false' +AllowShortLoopsOnASingleLine: 'false' +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: 'true' +AlwaysBreakTemplateDeclarations: 'true' +BinPackArguments: false +BinPackParameters: false +ExperimentalAutoDetectBinPacking: 'false' +BreakBeforeBinaryOperators: NonAssignment +BreakBeforeBraces: Custom +BreakBeforeTernaryOperators: 'false' +BreakConstructorInitializersBeforeComma: 'true' +ColumnLimit: '100' +ConstructorInitializerAllOnOneLineOrOnePerLine: 'false' +Cpp11BracedListStyle: 'true' +IndentCaseLabels: 'true' +IndentWidth: '4' +KeepEmptyLinesAtTheStartOfBlocks: 'true' +Language: Cpp +MaxEmptyLinesToKeep: '2' +NamespaceIndentation: Inner +ObjCSpaceBeforeProtocolList: 'true' +PointerAlignment: Left +SpaceAfterCStyleCast: 'false' +SpaceBeforeAssignmentOperators: 'true' +SpaceBeforeParens: Never +SpaceInEmptyParentheses: 'false' +SpacesBeforeTrailingComments: '2' +SpacesInAngles: 'false' +SpacesInCStyleCastParentheses: 'false' +SpacesInParentheses: 'false' +SpacesInSquareBrackets: 'false' +Standard: Cpp11 +TabWidth: '4' +UseTab: Never +SortIncludes: 'false' +ReflowComments: 'false' +BraceWrapping: { + AfterClass: 'true' + AfterControlStatement: 'true' + AfterEnum: 'true' + AfterFunction: 'true' + AfterNamespace: 'true' + AfterStruct: 'true' + AfterUnion: 'true' + BeforeCatch: 'true' + BeforeElse: 'true' + IndentBraces: 'false' + BeforeLambdaBody: 'true' +} +PenaltyExcessCharacter: 1 +PenaltyBreakBeforeFirstCallParameter: 40 +PenaltyBreakFirstLessLess: 1 +PenaltyBreakComment: 30 +PenaltyBreakString: 30 +PenaltyReturnTypeOnItsOwnLine: 9999 +BreakStringLiterals: false diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 9079757..7203552 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -12,12 +12,12 @@ jobs: steps: - name: clone - uses: actions/checkout@v4 # Action to check out your repository + uses: actions/checkout@v4 - name: dependencies run: | sudo apt-get update - sudo apt-get install -y build-essential cmake libtbb-dev + sudo apt-get install -y build-essential cmake - name: build working-directory: ${{ github.workspace }} diff --git a/CMakeLists.txt b/CMakeLists.txt index 179767f..a920589 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,18 +8,9 @@ list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") # Options option(MPM_TESTING_ONLY "Only compile lib and tests" OFF) -# Libigl also includes Eigen, etc +# Libigl also includes Eigen include(libigl) -# OpenMP -# (TO-DO: replace with TBB) -find_package(OpenMP) -if (OPENMP_FOUND) - set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - add_definitions(-DOMP_NESTED) -endif() - # Simulation library set(MPM_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src/LBFGS.hpp diff --git a/README.md b/README.md index 1a6cfbf..86bbbe2 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,15 @@ implemented from: [SIGGRAPH course notes](https://dl.acm.org/doi/abs/10.1145/289 [Optimization Integrator for Large Time Steps](https://dl.acm.org/doi/10.5555/2849517.2849523), and [The Affine Particle-In-Cell Method](https://dl.acm.org/doi/10.1145/2766996). -***2021 Update***: I've copied the original blog post to this readme, updated code formatting, added rendering with libigl, and switched to [mcloptlib](https://github.com/mattoverby/mcloptlib)'s LBFGS. -I wrote this in 2016 as an exercise, so please don't hold it against me :smile:. +***2021 Update***: I've copied the original blog post to this readme, updated code formatting, added rendering with libigl, and switched to a new LBFGS implementation. +I wrote this in 2016 as a learning exercise; no promises on efficiency or correctness :wink:. + +## Compiling + +``` +mkdir build && cd build && cmake -DCMAKE_BUILD_TYPE=Release .. && make -j +./sphere +``` ## Optimization-Based Material Point Method @@ -16,7 +23,6 @@ It has received a lot of attention in computer graphics, spurred mostly by the [ Some time around 2016 I coded up MPM to test out some research ideas, and wrote a little text to go along with it. My ideas didn't work out, but I put this code up on github instead of letting it collect dust on my hard drive. -It is not efficient nor bug free. The demo is nothing more than a big neo-Hookean blob dropped on a ground plane. ### Overview @@ -50,14 +56,14 @@ One thing to keep in mind: in MPM the velocity calculations happen on grid nodes Mapping to/from the grid is an important part of the material point method. MPM was originally an extension of Fluid Implicit Particle (FLIP), which in turn was an extension of the Particle in Cell (PIC) method. PIC and FLIP are the standard in fluid simulation, and were used in the above Gast paper. -Mike Seymour wrote [a good article on fxguide](https://www.fxguide.com/fxfeatured/the-science-of-fluid-sims/) that covers what PIC/FLIP is and why they are useful. +Mike Seymour wrote [a good article on fxguide](https://www.fxguide.com/fxfeatured/the-science-of-fluid-sims/) that covers what PIC/FLIP are and why they are useful. However, they suffer from a few well known limitations. Namely, PIC causes things to "clump" together, and FLIP can be unstable. A technique called the [affine particle in cell (APIC)](https://dl.acm.org/doi/10.1145/2766996) was recently introduced to combat these limitations. It's more stable and does a better job at conserving momentum. -In short, it's worth using. To do so, you only need to make a few subtle changes to the MPM algorithm. -The SIGGRAPH course notes explain what to change, but the Gast paper does not. +Implementation-wise, it's only a few subtle changes to the MPM algorithm, namely adding a per-particle affine matrix. +The SIGGRAPH course notes explain this and other changes. ## Implementation @@ -66,6 +72,6 @@ If you've ever done fluid simulation with PIC/FLIP, it should be straight forwar 1. It seemed to help with stability to have a cell size large enough to fit 8-16 particles. 2. Keep a fully-allocated grid as the domain boundary, but do minimization on a separate "active grid" list made up of pointers to grid nodes. -3. A "make grid loop" function for particle/grid mapping: given a particle, return the indices of grid nodes within the kernel radius, as well as kernel smoothing values/derivatives. I felt this was the trickiest part of the whole method. Much of this can be preallocated if you have a lot of memory, but I found (in 3D) allocating/deallocating on the spot was helpful and easier to debug. +3. A "make grid loop" function for particle/grid mapping: given a particle, return the indices of grid nodes within the kernel radius, as well as kernel smoothing values/derivatives. I felt this was the trickiest part of the whole method. Much of this can be preallocated but I found (in 3D) allocating/deallocating on the spot was helpful and easier to debug. The code uses libigl for rendering and Eigen for matrix/vector calculations. These are automatically downloaded with cmake. diff --git a/cmake/DownloadProject.CMakeLists.cmake.in b/cmake/DownloadProject.CMakeLists.cmake.in deleted file mode 100644 index 89be4fd..0000000 --- a/cmake/DownloadProject.CMakeLists.cmake.in +++ /dev/null @@ -1,17 +0,0 @@ -# Distributed under the OSI-approved MIT License. See accompanying -# file LICENSE or https://github.com/Crascit/DownloadProject for details. - -cmake_minimum_required(VERSION 2.8.2) - -project(${DL_ARGS_PROJ}-download NONE) - -include(ExternalProject) -ExternalProject_Add(${DL_ARGS_PROJ}-download - ${DL_ARGS_UNPARSED_ARGUMENTS} - SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" - BINARY_DIR "${DL_ARGS_BINARY_DIR}" - CONFIGURE_COMMAND "" - BUILD_COMMAND "" - INSTALL_COMMAND "" - TEST_COMMAND "" -) diff --git a/cmake/DownloadProject.cmake b/cmake/DownloadProject.cmake deleted file mode 100644 index e300f42..0000000 --- a/cmake/DownloadProject.cmake +++ /dev/null @@ -1,182 +0,0 @@ -# Distributed under the OSI-approved MIT License. See accompanying -# file LICENSE or https://github.com/Crascit/DownloadProject for details. -# -# MODULE: DownloadProject -# -# PROVIDES: -# download_project( PROJ projectName -# [PREFIX prefixDir] -# [DOWNLOAD_DIR downloadDir] -# [SOURCE_DIR srcDir] -# [BINARY_DIR binDir] -# [QUIET] -# ... -# ) -# -# Provides the ability to download and unpack a tarball, zip file, git repository, -# etc. at configure time (i.e. when the cmake command is run). How the downloaded -# and unpacked contents are used is up to the caller, but the motivating case is -# to download source code which can then be included directly in the build with -# add_subdirectory() after the call to download_project(). Source and build -# directories are set up with this in mind. -# -# The PROJ argument is required. The projectName value will be used to construct -# the following variables upon exit (obviously replace projectName with its actual -# value): -# -# projectName_SOURCE_DIR -# projectName_BINARY_DIR -# -# The SOURCE_DIR and BINARY_DIR arguments are optional and would not typically -# need to be provided. They can be specified if you want the downloaded source -# and build directories to be located in a specific place. The contents of -# projectName_SOURCE_DIR and projectName_BINARY_DIR will be populated with the -# locations used whether you provide SOURCE_DIR/BINARY_DIR or not. -# -# The DOWNLOAD_DIR argument does not normally need to be set. It controls the -# location of the temporary CMake build used to perform the download. -# -# The PREFIX argument can be provided to change the base location of the default -# values of DOWNLOAD_DIR, SOURCE_DIR and BINARY_DIR. If all of those three arguments -# are provided, then PREFIX will have no effect. The default value for PREFIX is -# CMAKE_BINARY_DIR. -# -# The QUIET option can be given if you do not want to show the output associated -# with downloading the specified project. -# -# In addition to the above, any other options are passed through unmodified to -# ExternalProject_Add() to perform the actual download, patch and update steps. -# The following ExternalProject_Add() options are explicitly prohibited (they -# are reserved for use by the download_project() command): -# -# CONFIGURE_COMMAND -# BUILD_COMMAND -# INSTALL_COMMAND -# TEST_COMMAND -# -# Only those ExternalProject_Add() arguments which relate to downloading, patching -# and updating of the project sources are intended to be used. Also note that at -# least one set of download-related arguments are required. -# -# If using CMake 3.2 or later, the UPDATE_DISCONNECTED option can be used to -# prevent a check at the remote end for changes every time CMake is run -# after the first successful download. See the documentation of the ExternalProject -# module for more information. It is likely you will want to use this option if it -# is available to you. Note, however, that the ExternalProject implementation contains -# bugs which result in incorrect handling of the UPDATE_DISCONNECTED option when -# using the URL download method or when specifying a SOURCE_DIR with no download -# method. Fixes for these have been created, the last of which is scheduled for -# inclusion in CMake 3.8.0. Details can be found here: -# -# https://gitlab.kitware.com/cmake/cmake/commit/bdca68388bd57f8302d3c1d83d691034b7ffa70c -# https://gitlab.kitware.com/cmake/cmake/issues/16428 -# -# If you experience build errors related to the update step, consider avoiding -# the use of UPDATE_DISCONNECTED. -# -# EXAMPLE USAGE: -# -# include(DownloadProject) -# download_project(PROJ googletest -# GIT_REPOSITORY https://github.com/google/googletest.git -# GIT_TAG master -# UPDATE_DISCONNECTED 1 -# QUIET -# ) -# -# add_subdirectory(${googletest_SOURCE_DIR} ${googletest_BINARY_DIR}) -# -#======================================================================================== - - -set(_DownloadProjectDir "${CMAKE_CURRENT_LIST_DIR}") - -include(CMakeParseArguments) - -function(download_project) - - set(options QUIET) - set(oneValueArgs - PROJ - PREFIX - DOWNLOAD_DIR - SOURCE_DIR - BINARY_DIR - # Prevent the following from being passed through - CONFIGURE_COMMAND - BUILD_COMMAND - INSTALL_COMMAND - TEST_COMMAND - ) - set(multiValueArgs "") - - cmake_parse_arguments(DL_ARGS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) - - # Hide output if requested - if (DL_ARGS_QUIET) - set(OUTPUT_QUIET "OUTPUT_QUIET") - else() - unset(OUTPUT_QUIET) - message(STATUS "Downloading/updating ${DL_ARGS_PROJ}") - endif() - - # Set up where we will put our temporary CMakeLists.txt file and also - # the base point below which the default source and binary dirs will be. - # The prefix must always be an absolute path. - if (NOT DL_ARGS_PREFIX) - set(DL_ARGS_PREFIX "${CMAKE_BINARY_DIR}") - else() - get_filename_component(DL_ARGS_PREFIX "${DL_ARGS_PREFIX}" ABSOLUTE - BASE_DIR "${CMAKE_CURRENT_BINARY_DIR}") - endif() - if (NOT DL_ARGS_DOWNLOAD_DIR) - set(DL_ARGS_DOWNLOAD_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-download") - endif() - - # Ensure the caller can know where to find the source and build directories - if (NOT DL_ARGS_SOURCE_DIR) - set(DL_ARGS_SOURCE_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-src") - endif() - if (NOT DL_ARGS_BINARY_DIR) - set(DL_ARGS_BINARY_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-build") - endif() - set(${DL_ARGS_PROJ}_SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" PARENT_SCOPE) - set(${DL_ARGS_PROJ}_BINARY_DIR "${DL_ARGS_BINARY_DIR}" PARENT_SCOPE) - - # The way that CLion manages multiple configurations, it causes a copy of - # the CMakeCache.txt to be copied across due to it not expecting there to - # be a project within a project. This causes the hard-coded paths in the - # cache to be copied and builds to fail. To mitigate this, we simply - # remove the cache if it exists before we configure the new project. It - # is safe to do so because it will be re-generated. Since this is only - # executed at the configure step, it should not cause additional builds or - # downloads. - file(REMOVE "${DL_ARGS_DOWNLOAD_DIR}/CMakeCache.txt") - - # Create and build a separate CMake project to carry out the download. - # If we've already previously done these steps, they will not cause - # anything to be updated, so extra rebuilds of the project won't occur. - # Make sure to pass through CMAKE_MAKE_PROGRAM in case the main project - # has this set to something not findable on the PATH. - configure_file("${_DownloadProjectDir}/DownloadProject.CMakeLists.cmake.in" - "${DL_ARGS_DOWNLOAD_DIR}/CMakeLists.txt") - execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" - -D "CMAKE_MAKE_PROGRAM:FILE=${CMAKE_MAKE_PROGRAM}" - . - RESULT_VARIABLE result - ${OUTPUT_QUIET} - WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" - ) - if(result) - message(FATAL_ERROR "CMake step for ${DL_ARGS_PROJ} failed: ${result}") - endif() - execute_process(COMMAND ${CMAKE_COMMAND} --build . - RESULT_VARIABLE result - ${OUTPUT_QUIET} - WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" - ) - if(result) - message(FATAL_ERROR "Build step for ${DL_ARGS_PROJ} failed: ${result}") - endif() - -endfunction() diff --git a/src/Interp.hpp b/src/Interp.hpp index 9520d60..1f374c6 100644 --- a/src/Interp.hpp +++ b/src/Interp.hpp @@ -1,78 +1,110 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -#ifndef INTERP_HPP -#define INTERP_HPP 1 +#ifndef MPM_INTERP_HPP +#define MPM_INTERP_HPP 1 #include -#include namespace mpm { - // Cubic B-spline - static inline double cspline(double x) - { - if (x==0.0) { x=1e-12; } - x = fabs(x); - if (x < 1.0) { return x*x*(x*0.5 - 1.0) + 2.0/3.0; } - else if (x < 2.0) { return x*(x*(-x/6.0 + 1.0) - 2.0) + 4.0/3.0; } - return 0.0; - } +// Cubic B-spline +static inline double cspline(double x) +{ + if(x == 0.0) + { + x = 1e-12; + } + x = std::abs(x); + if(x < 1.0) + { + return x * x * (x * 0.5 - 1.0) + 2.0 / 3.0; + } + else if(x < 2.0) + { + return x * (x * (-x / 6.0 + 1.0) - 2.0) + 4.0 / 3.0; + } + return 0.0; +} - // Slope of cubic spline - static inline double d_cspline(double x) - { - if (x==0.0) { x=1e-12; } - double abs_x = fabs(x); - if (abs_x < 1.0) { return 1.5*x*abs_x - 2.0*x; } - else if (x < 2.0) { return -x*abs_x*0.5 + 2.0*x - 2.0*x/abs_x; } - return 0.0; - } +// Slope of cubic spline +static inline double d_cspline(double x) +{ + if(x == 0.0) + { + x = 1e-12; + } + double abs_x = std::abs(x); + if(abs_x < 1.0) + { + return 1.5 * x * abs_x - 2.0 * x; + } + else if(x < 2.0) + { + return -x * abs_x * 0.5 + 2.0 * x - 2.0 * x / abs_x; + } + return 0.0; +} - // Quadratic spline - static inline double qspline( double x ) - { - double fx = fabs(x); - if (fx < 0.5) { return ( 0.75 - x*x ); } - else if (fx < 1.5) { return ( 0.5*x*x - 1.5*fx + 9.0/8.0 ); } - return 0.0; - } +// Quadratic spline +static inline double qspline(double x) +{ + double fx = std::abs(x); + if(fx < 0.5) + { + return (0.75 - x * x); + } + else if(fx < 1.5) + { + return (0.5 * x * x - 1.5 * fx + 9.0 / 8.0); + } + return 0.0; +} - // Slope of quadratic spline - static inline double d_qspline( double x ) - { - double fx = fabs(x); - if (fx < 0.5) { return ( -2.0 * x ); } - else if (fx < 1.5) - { - if (x < 0.0) { return fabs( fx - 1.5 ); } - else { return (fx - 1.5); } - } - return 0.0; - } +// Slope of quadratic spline +static inline double d_qspline(double x) +{ + double fx = std::abs(x); + if(fx < 0.5) + { + return (-2.0 * x); + } + else if(fx < 1.5) + { + if(x < 0.0) + { + return std::abs(fx - 1.5); + } + else + { + return (fx - 1.5); + } + } + return 0.0; +} - static inline bool isreal(double n) - { - return !std::isnan(n) && !std::isinf(n); - } +static inline bool isreal(double n) +{ + return !std::isnan(n) && !std::isinf(n); +} -} // end namespace mpm +} // end namespace mpm #endif diff --git a/src/LBFGS.hpp b/src/LBFGS.hpp index 839ebb3..2381215 100644 --- a/src/LBFGS.hpp +++ b/src/LBFGS.hpp @@ -31,410 +31,413 @@ namespace mcl // g = ... (reuse computation for objective) // } // return objective; -// }; +// }; // template class LBFGS { -public: + public: typedef typename MatrixType::Scalar Scalar; - struct Options - { - int min_iters; - int max_iters; - Scalar abs_tol; // absolute tol if converged(...) not set - Scalar rel_tol; // relative tol if converged(...) not set - int M; // history window size - Scalar gamma; // init Hessian = gamma * I - Options() : - min_iters(0), - max_iters(100), - abs_tol(1e-5), - rel_tol(0), - M(6), - gamma(1) - {} - } options; - - LBFGS(); - - // Output from the last call to minimize(...) - int iters() const { return num_iters; } - Scalar gamma() const { return gamma_k; } - - // Resizes buffers and sets to zero. - // Called during minimize(...) ONLY if there is a change in dof or M. - // Otherwise it's assumed you're picking up where you left off - // on the previous call to minimize(...) - void reset(int rows, int cols); - - // Required: - // computes objective value and gradient - // obj = gradient(x, g) - // If the g arg is not sized, don't compute gradient - std::function gradient; - - // Optional: - // Returns true if the solver should exit, default uses ||g|| converged; - - // Optional: - // Linesearch function, default uses bracketing weak wolfe (slow!) - // obj_k1 = (x, grad, descent, alpha) - // Returns new objective value and updates both x AND gradient - std::function linesearch; - - // Optional: - // Filter descent direction, p = B(p) - // Otherwise p = gamma_k * p is used. - std::function filter; - - // Calls initialize(x) once and iterate(x) until converged - Scalar minimize(MatrixType& x); - - // Initialize the solver - void initialize(MatrixType& x); - - // Take an iteration - // Returns objective - Scalar iterate(MatrixType& x); - - // i.e. bisection with weak Wolfe conditions - Scalar bracketing_weakwolfe( - MatrixType& x, - MatrixType& grad, - const MatrixType& p, - Scalar &alpha) const; - - // Used if converged not set - // Returns true if: - // grad.norm <= abs_tol - // or - // grad.norm() <= rel_tol * x.norm() - bool default_converged( - Scalar curr_obj, - const MatrixType& xprev, - const MatrixType& x, - const MatrixType& grad) const; - - Scalar inner( - const MatrixType &a, - const MatrixType &b) const; - -protected: - bool initialized; - int num_iters, max_iters, k; - Scalar gamma_k, obj_0, obj_k; - std::vector s; - std::vector y; - Eigen::Matrix alpha; - Eigen::Matrix rho; - MatrixType grad; - MatrixType q; - MatrixType descent; - MatrixType grad_old; - MatrixType x_old; - MatrixType x_last; - MatrixType s_temp; - MatrixType y_temp; - -}; // end class LBFGS + struct Options + { + int min_iters; + int max_iters; + Scalar abs_tol; // absolute tol if converged(...) not set + Scalar rel_tol; // relative tol if converged(...) not set + int M; // history window size + Scalar gamma; // init Hessian = gamma * I + Options() + : min_iters(0) + , max_iters(100) + , abs_tol(1e-5) + , rel_tol(0) + , M(6) + , gamma(1) + { + } + } options; + + LBFGS(); + + // Output from the last call to minimize(...) + int iters() const { return num_iters; } + Scalar gamma() const { return gamma_k; } + + // Resizes buffers and sets to zero. + // Called during minimize(...) ONLY if there is a change in dof or M. + // Otherwise it's assumed you're picking up where you left off + // on the previous call to minimize(...) + void reset(int rows, int cols); + + // Required: + // computes objective value and gradient + // obj = gradient(x, g) + // If the g arg is not sized, don't compute gradient + std::function gradient; + + // Optional: + // Returns true if the solver should exit, default uses ||g|| converged; + + // Optional: + // Linesearch function, default uses bracketing weak wolfe (slow!) + // obj_k1 = (x, grad, descent, alpha) + // Returns new objective value and updates both x AND gradient + std::function linesearch; + + // Optional: + // Filter descent direction, p = B(p) + // Otherwise p = gamma_k * p is used. + std::function filter; + + // Calls initialize(x) once and iterate(x) until converged + Scalar minimize(MatrixType& x); + + // Initialize the solver + void initialize(MatrixType& x); + + // Take an iteration + // Returns objective + Scalar iterate(MatrixType& x); + + // i.e. bisection with weak Wolfe conditions + Scalar bracketing_weakwolfe(MatrixType& x, MatrixType& grad, const MatrixType& p, Scalar& alpha) const; + + // Used if converged not set + // Returns true if: + // grad.norm <= abs_tol + // or + // grad.norm() <= rel_tol * x.norm() + bool default_converged(Scalar curr_obj, const MatrixType& xprev, const MatrixType& x, const MatrixType& grad) const; + + Scalar inner(const MatrixType& a, const MatrixType& b) const; + + protected: + bool initialized; + int num_iters, max_iters, k; + Scalar gamma_k, obj_0, obj_k; + std::vector s; + std::vector y; + Eigen::Matrix alpha; + Eigen::Matrix rho; + MatrixType grad; + MatrixType q; + MatrixType descent; + MatrixType grad_old; + MatrixType x_old; + MatrixType x_last; + MatrixType s_temp; + MatrixType y_temp; + +}; // end class LBFGS // // Implementation // template -LBFGS::LBFGS() : - initialized(false), - num_iters(0), - max_iters(0), - k(0), - gamma_k(1), - obj_0(0), - obj_k(0) - {} +LBFGS::LBFGS() + : initialized(false) + , num_iters(0) + , max_iters(0) + , k(0) + , gamma_k(1) + , obj_0(0) + , obj_k(0) +{ +} template void LBFGS::reset(int rows, int cols) { using namespace Eigen; - num_iters = 0; - max_iters = 0; - k = 0; - gamma_k = 1; - obj_0 = std::numeric_limits::max(); - obj_k = std::numeric_limits::max(); - int M = options.M; - s = std::vector(M, MatrixType::Zero(rows,cols)); - y = std::vector(M, MatrixType::Zero(rows,cols)); - alpha = VectorXd::Zero(M); - rho = VectorXd::Zero(M); - grad = MatrixType::Zero(rows, cols); - q = MatrixType::Zero(rows, cols); // inv descent - descent = q; - grad_old = MatrixType::Zero(rows, cols); - x_old = MatrixType::Zero(rows, cols); - x_last = MatrixType::Zero(rows, cols); - s_temp = MatrixType::Zero(rows, cols); - y_temp = MatrixType::Zero(rows, cols); + num_iters = 0; + max_iters = 0; + k = 0; + gamma_k = 1; + obj_0 = std::numeric_limits::max(); + obj_k = std::numeric_limits::max(); + int M = options.M; + s = std::vector(M, MatrixType::Zero(rows, cols)); + y = std::vector(M, MatrixType::Zero(rows, cols)); + alpha = VectorXd::Zero(M); + rho = VectorXd::Zero(M); + grad = MatrixType::Zero(rows, cols); + q = MatrixType::Zero(rows, cols); // inv descent + descent = q; + grad_old = MatrixType::Zero(rows, cols); + x_old = MatrixType::Zero(rows, cols); + x_last = MatrixType::Zero(rows, cols); + s_temp = MatrixType::Zero(rows, cols); + y_temp = MatrixType::Zero(rows, cols); } // Returns number of iterations used template -typename LBFGS::Scalar -LBFGS::minimize(MatrixType &x) +typename LBFGS::Scalar LBFGS::minimize(MatrixType& x) { - initialize(x); + initialize(x); - // Did we start at the initializer? - if (num_iters >= options.min_iters && - converged(obj_k,x_last,x,grad)) - { - num_iters = 1; - return obj_k; - } + // Did we start at the initializer? + if(num_iters >= options.min_iters && converged(obj_k, x_last, x, grad)) + { + num_iters = 1; + return obj_k; + } - for (; k void LBFGS::initialize(MatrixType& x) { - initialized = true; - - if (gradient == nullptr) { - throw std::runtime_error("no gradient function"); - } - - if (converged == nullptr) - { - using namespace std::placeholders; - converged = std::bind(&LBFGS::default_converged, this, _1, _2, _3, _4); - } - - if (linesearch == nullptr) - { - using namespace std::placeholders; - linesearch = std::bind(&LBFGS::bracketing_weakwolfe, this, _1, _2, _3, _4); - } - - // Resize/resize variables? - if (alpha.rows() != options.M || - grad.rows() != x.rows() || - grad.cols() != x.cols()) { - reset(x.rows(), x.cols()); - } - - num_iters = 0; - max_iters = std::max(options.min_iters, options.max_iters); - gamma_k = options.gamma; - obj_0 = gradient(x, grad); - obj_k = obj_0; - k = 0; + initialized = true; + + if(gradient == nullptr) + { + throw std::runtime_error("no gradient function"); + } + + if(converged == nullptr) + { + using namespace std::placeholders; + converged = std::bind(&LBFGS::default_converged, this, _1, _2, _3, _4); + } + + if(linesearch == nullptr) + { + using namespace std::placeholders; + linesearch = std::bind(&LBFGS::bracketing_weakwolfe, this, _1, _2, _3, _4); + } + + // Resize/resize variables? + if(alpha.rows() != options.M || grad.rows() != x.rows() || grad.cols() != x.cols()) + { + reset(x.rows(), x.cols()); + } + + num_iters = 0; + max_iters = std::max(options.min_iters, options.max_iters); + gamma_k = options.gamma; + obj_0 = gradient(x, grad); + obj_k = obj_0; + k = 0; } template -typename LBFGS::Scalar -LBFGS::iterate(MatrixType& x) +typename LBFGS::Scalar LBFGS::iterate(MatrixType& x) { - if (!initialized) { - throw std::runtime_error("not initialized"); - } - - x_old = x; - grad_old = grad; - q = grad; - num_iters++; - - // - // Two-loop recursion - // - { - // L-BFGS first - loop recursion - int iter = std::min(options.M, k); - for(int i = iter - 1; i >= 0; --i) - { - Scalar denom = inner(s[i], y[i]); - if (std::abs(denom) <= 0.0) - { - rho(i) = 0; - alpha(i) = 0; - continue; - } - rho(i) = 1.0 / denom; - alpha(i) = rho(i)*inner(s[i], q); - q -= alpha(i) * y[i]; - } - - if (filter != nullptr) { filter(q); } - else { q = gamma_k*q; } - - // L-BFGS second - loop recursion - for(int i = 0; i < iter; ++i) - { - Scalar beta = rho(i)*inner(q, y[i]); - q += (alpha(i) - beta)*s[i]; - } - } - - // - // Perform step - // - { - // If our hess approx is bad and we start going in - // the wrong direction, restart memory - Scalar step_size = 1.0; - Scalar dir = inner(q, grad); - if (dir <= 0) - { - q = grad; - max_iters -= k; // Restart memory - k = 0; - step_size = std::min(1.0, 1.0 / grad.template lpNorm() ); - } - - // We've hit local minima, we have to exit - if (q.squaredNorm() <= 0.0) { - return obj_k; - } - - descent = -q; - x_last = x; - obj_k = linesearch(x, grad, descent, step_size); - if (num_iters >= options.min_iters && converged(obj_k,x_last,x,grad)) { - return obj_k; - } - } - - // - // Correction term - // - { - s_temp = x - x_old; - y_temp = grad - grad_old; - - // update the history - if (k < options.M) - { - s[k] = s_temp; - y[k] = y_temp; - } - else - { - for (int i=0; i 0.0) - { - gamma_k = inner(s_temp, y_temp) / denom; - } - } - - return obj_k; + if(!initialized) + { + throw std::runtime_error("not initialized"); + } + + x_old = x; + grad_old = grad; + q = grad; + num_iters++; + + // + // Two-loop recursion + // + { + // L-BFGS first - loop recursion + int iter = std::min(options.M, k); + for(int i = iter - 1; i >= 0; --i) + { + Scalar denom = inner(s[i], y[i]); + if(std::abs(denom) <= 0.0) + { + rho(i) = 0; + alpha(i) = 0; + continue; + } + rho(i) = 1.0 / denom; + alpha(i) = rho(i) * inner(s[i], q); + q -= alpha(i) * y[i]; + } + + if(filter != nullptr) + { + filter(q); + } + else + { + q = gamma_k * q; + } + + // L-BFGS second - loop recursion + for(int i = 0; i < iter; ++i) + { + Scalar beta = rho(i) * inner(q, y[i]); + q += (alpha(i) - beta) * s[i]; + } + } + + // + // Perform step + // + { + // If our hess approx is bad and we start going in + // the wrong direction, restart memory + Scalar step_size = 1.0; + Scalar dir = inner(q, grad); + if(dir <= 0) + { + q = grad; + max_iters -= k; // Restart memory + k = 0; + step_size = std::min(1.0, 1.0 / grad.template lpNorm()); + } + + // We've hit local minima, we have to exit + if(q.squaredNorm() <= 0.0) + { + return obj_k; + } + + descent = -q; + x_last = x; + obj_k = linesearch(x, grad, descent, step_size); + if(num_iters >= options.min_iters && converged(obj_k, x_last, x, grad)) + { + return obj_k; + } + } + + // + // Correction term + // + { + s_temp = x - x_old; + y_temp = grad - grad_old; + + // update the history + if(k < options.M) + { + s[k] = s_temp; + y[k] = y_temp; + } + else + { + for(int i = 0; i < options.M - 1; ++i) + { + s[i] = s[i + 1]; + y[i] = y[i + 1]; + } + s.back() = s_temp; + y.back() = y_temp; + } + + Scalar denom = inner(y_temp, y_temp); + if(std::abs(denom) > 0.0) + { + gamma_k = inner(s_temp, y_temp) / denom; + } + } + + return obj_k; } template -typename LBFGS::Scalar -LBFGS::bracketing_weakwolfe( - MatrixType& x, - MatrixType& g, - const MatrixType &drt, - Scalar &step) const +typename LBFGS::Scalar LBFGS::bracketing_weakwolfe(MatrixType& x, + MatrixType& g, + const MatrixType& drt, + Scalar& step) const { using namespace Eigen; - Scalar c1 = 10e-4; - Scalar c2 = 0.9; - if (step <= 0.0) { step = 1.0; } + Scalar c1 = 10e-4; + Scalar c2 = 0.9; + if(step <= 0.0) + { + step = 1.0; + } int x_cols = x.cols(); - g = MatrixType::Zero(x.rows(), x_cols); - MatrixType xp = x; - const Scalar fx_init = gradient(x,g); - Scalar fx = fx_init; - - Scalar dg_init = inner(g, drt); - if (dg_init > 0) { + g = MatrixType::Zero(x.rows(), x_cols); + MatrixType xp = x; + const Scalar fx_init = gradient(x, g); + Scalar fx = fx_init; + + Scalar dg_init = inner(g, drt); + if(dg_init > 0) + { throw std::runtime_error("direction increases objective"); } - const Scalar test_decr = c1 * dg_init; - Scalar lower = 0; - Scalar upper = std::numeric_limits::infinity(); - int maxiter = 200; - int iter = 0; - for (iter = 0; iter < maxiter; iter++) - { - x = xp + step * drt; - fx = gradient(x, g); - if (fx > fx_init + step * test_decr) { // Armijo rule - upper = step; - } - else - { - Scalar dg = inner(g, drt); - if(dg < c2 * dg_init){ // Weak wolfe - lower = step; - } else { - break; // both met - } - } - step = std::isinf(upper) ? 2*step : lower/2 + upper/2; - } - return fx; - -} // end linesearch + const Scalar test_decr = c1 * dg_init; + Scalar lower = 0; + Scalar upper = std::numeric_limits::infinity(); + int maxiter = 200; + int iter = 0; + for(iter = 0; iter < maxiter; iter++) + { + x = xp + step * drt; + fx = gradient(x, g); + if(fx > fx_init + step * test_decr) + { // Armijo rule + upper = step; + } + else + { + Scalar dg = inner(g, drt); + if(dg < c2 * dg_init) + { // Weak wolfe + lower = step; + } + else + { + break; // both met + } + } + step = std::isinf(upper) ? 2 * step : lower / 2 + upper / 2; + } + return fx; + +} // end linesearch template -bool LBFGS::default_converged( - Scalar curr_obj, - const MatrixType& xprev, - const MatrixType& x, - const MatrixType& g) const +bool LBFGS::default_converged(Scalar curr_obj, + const MatrixType& xprev, + const MatrixType& x, + const MatrixType& g) const { - (void)(curr_obj); - (void)(xprev); - Scalar gnorm = g.norm(); - if (gnorm <= options.abs_tol) { - return true; - } - Scalar xnorm = x.norm(); - if (gnorm < options.rel_tol * xnorm) { - return true; - } - return false; + (void)(curr_obj); + (void)(xprev); + Scalar gnorm = g.norm(); + if(gnorm <= options.abs_tol) + { + return true; + } + Scalar xnorm = x.norm(); + if(gnorm < options.rel_tol * xnorm) + { + return true; + } + return false; } template -typename LBFGS::Scalar -LBFGS::inner( - const MatrixType &a, - const MatrixType &b) const +typename LBFGS::Scalar LBFGS::inner(const MatrixType& a, const MatrixType& b) const { - int cols = std::min(a.cols(), b.cols()); - Scalar dot = 0; - for (int i=0; i +#include namespace mpm { -// -// Declare static members -// These are default, usually overwritten by the solver. -// -Eigen::Vector3d MPM::cellsize(0.25,0.25,0.25); -Eigen::Vector3i MPM::gridsize(32,32,32); -double MPM::timestep_s = 0.04; - // Pass in (normalized) grid index and returns its index in the grid vector. -int MPM::grid_index(const Eigen::Vector3i &pos) +int MPM::grid_index(const Settings& settings, const Eigen::Vector3i& pos) { - return (pos[2]*gridsize[2]*gridsize[1]) + (pos[1]*gridsize[0]) + pos[0]; + return (pos[2] * settings.gridsize[2] * settings.gridsize[1]) + (pos[1] * settings.gridsize[0]) + pos[0]; } // Inverse of grid_index (position in grid, NOT world space!) -Eigen::Vector3i MPM::grid_pos(int grid_idx) +Eigen::Vector3i MPM::grid_pos(const Settings& settings, int grid_idx) { - int z = grid_idx / (gridsize[0]*gridsize[1]); - grid_idx -= (z * gridsize[0]*gridsize[1]); - int y = grid_idx / gridsize[0]; - int x = grid_idx % gridsize[0]; - return Eigen::Vector3i(x, y, z); + int z = grid_idx / (settings.gridsize[0] * settings.gridsize[1]); + grid_idx -= (z * settings.gridsize[0] * settings.gridsize[1]); + int y = grid_idx / settings.gridsize[0]; + int x = grid_idx % settings.gridsize[0]; + return Eigen::Vector3i(x, y, z); } - // // Clear the grid // -void MPM::clear_grid(std::vector &active_grid) +void MPM::clear_grid(std::vector& active_grid) { + igl::parallel_for( + active_grid.size(), + [&](int i) + { + GridNode* g = active_grid[i]; -#pragma omp parallel for - for(int i=0; iv = Eigen::Vector3d(0, 0, 0); + g->m = 0.0; + g->active_idx = 0; - g->v = Eigen::Vector3d(0,0,0); - g->m = 0.0; - g->active_idx = 0; + g->particles.clear(); + g->wip.clear(); + g->dwip.clear(); + }, + MPM_MIN_PARALLEL); - g->particles.clear(); - g->wip.clear(); - g->dwip.clear(); - } + active_grid.clear(); - active_grid.clear(); - -} // end clear grid +} // end clear grid // // Particle to grid: masses and interpolation weights // -void MPM::p2g_mass(std::vector &particles, std::vector &grid, std::vector &active_grid) +void MPM::p2g_mass(const Settings& settings, + std::vector& particles, + std::vector& grid, + std::vector& active_grid) { - using namespace Eigen; - -#pragma omp parallel for - for(int i=0; iD.setZero(); - - std::vector grid_nodes; - std::vector wip; - std::vector dwip; - MPM::make_grid_loop(p, grid_nodes, wip, dwip); - for(int j=0; jD += wip[j] * (xixp) * (xixp).transpose(); // eq 174 course notes - - // Atomics are bad and all, but this way works fine for now. -#pragma omp critical - { // add mass to grid - int grid_idx = grid_nodes[j]; - grid[ grid_idx ]->m += wip[j]*p->m; - grid[ grid_idx ]->particles.push_back(p); - grid[ grid_idx ]->wip.push_back(wip[j]); - grid[ grid_idx ]->dwip.push_back(dwip[j]); - } - - } // end loop grid nodes - - } // end loop particles - - // Create the active grid - for(int i=0; im > 0.0) - { - grid[i]->active_idx = active_grid.size(); - active_grid.push_back(grid[i]); - } - } - -} // end p2g mass and weights + using namespace Eigen; + + igl::parallel_for( + particles.size(), + [&](int i) + { + Particle* p = particles[i]; + p->D.setZero(); + + std::vector grid_nodes; + std::vector wip; + std::vector dwip; + MPM::make_grid_loop(settings, p, grid_nodes, wip, dwip); + for(size_t j = 0; j < grid_nodes.size(); ++j) + { + // Compute helper matrix + Eigen::Vector3d xixp = xi_xp(settings, grid_nodes[j], p); + p->D += wip[j] * (xixp) * (xixp).transpose(); // eq 174 course notes + + { // add mass to grid + int grid_idx = grid_nodes[j]; + std::lock_guard lock(grid[grid_idx]->mutex); + grid[grid_idx]->m += wip[j] * p->m; + grid[grid_idx]->particles.push_back(p); + grid[grid_idx]->wip.push_back(wip[j]); + grid[grid_idx]->dwip.push_back(dwip[j]); + } + + } // end loop grid nodes + }, + MPM_MIN_PARALLEL); + + // Create the active grid + for(size_t i = 0; i < grid.size(); ++i) + { + if(grid[i]->m > 0.0) + { + grid[i]->active_idx = active_grid.size(); + active_grid.push_back(grid[i]); + } + } + +} // end p2g mass and weights // // Calculate (initial) particle volumes // -void MPM::init_volumes(std::vector &particles, const std::vector &grid) +void MPM::init_volumes(const Settings& settings, std::vector& particles, const std::vector& grid) { -#pragma omp parallel for - for(int i=0; i grid_nodes; - std::vector wip; - std::vector dwip; - MPM::make_grid_loop(p, grid_nodes, wip, dwip); - for(int j=0; jm; - } + igl::parallel_for( + particles.size(), + [&](int i) + { + Particle* p = particles[i]; + double dens = 0.0; - dens /= (cellsize[0]*cellsize[1]*cellsize[2]); // density - assert(dens > 0.0); - p->vol = p->m / dens; // volume + std::vector grid_nodes; + std::vector wip; + std::vector dwip; + MPM::make_grid_loop(settings, p, grid_nodes, wip, dwip); + for(size_t j = 0; j < grid_nodes.size(); ++j) + { + dens += wip[j] * grid[grid_nodes[j]]->m; + } - } // end loop particles + dens /= (settings.cellsize[0] * settings.cellsize[1] * settings.cellsize[2]); // density + assert(dens > 0.0); + p->vol = p->m / dens; // volume + }, + MPM_MIN_PARALLEL); -} // end compute particle volumes +} // end compute particle volumes // // Particle to Grid Velocity // -void MPM::p2g_velocity(const std::vector &particles, std::vector &active_grid) +void MPM::p2g_velocity(const Settings& settings, + const std::vector& particles, + std::vector& active_grid) { - using namespace Eigen; - -#pragma omp parallel for - for(int i=0; iglobal_idx); + using namespace Eigen; - for(int j=0; jparticles.size(); ++j) - { - Particle *p = active_grid[i]->particles[j]; - Matrix3d Dp_inv = p->D.inverse(); + igl::parallel_for( + active_grid.size(), + [&](int i) + { + Eigen::Vector3i gridx = MPM::grid_pos(settings, active_grid[i]->global_idx); - Vector3d xixp(gridx[0]*cellsize[0]-p->x[0], gridx[1]*cellsize[1]-p->x[1], gridx[2]*cellsize[2]-p->x[2]); - Vector3d newv = p->v + p->B*Dp_inv * (xixp); - active_grid[i]->v += active_grid[i]->wip[j]*p->m*newv; // eq 173 course notes + for(size_t j = 0; j < active_grid[i]->particles.size(); ++j) + { + Particle* p = active_grid[i]->particles[j]; + Matrix3d Dp_inv = p->D.inverse(); - } // end loop neighborhood + Vector3d xixp(gridx[0] * settings.cellsize[0] - p->x[0], + gridx[1] * settings.cellsize[1] - p->x[1], + gridx[2] * settings.cellsize[2] - p->x[2]); + Vector3d newv = p->v + p->B * Dp_inv * (xixp); + active_grid[i]->v += active_grid[i]->wip[j] * p->m * newv; // eq 173 course notes - // Fix Velocity (remove mass) - assert(active_grid[i]->m > 0.0); - active_grid[i]->v /= active_grid[i]->m; + } // end loop neighborhood - } // end loop grid + // Fix Velocity (remove mass) + assert(active_grid[i]->m > 0.0); + active_grid[i]->v /= active_grid[i]->m; + }, + MPM_MIN_PARALLEL); -} // end p2g velocity +} // end p2g velocity // // Update the particle deformation gradient (step 6 course notes) // -void MPM::g2p_deformation(std::vector &particles, const std::vector &grid) +void MPM::g2p_deformation(const Settings& settings, + std::vector& particles, + const std::vector& grid) { - using namespace Eigen; + using namespace Eigen; -#pragma omp parallel for - for(int i=0; i grid_nodes; - std::vector wip; - std::vector dwip; - MPM::make_grid_loop(p, grid_nodes, wip, dwip); - for(int j=0; jv * dwip[j].transpose(); - } + std::vector grid_nodes; + std::vector wip; + std::vector dwip; + MPM::make_grid_loop(settings, p, grid_nodes, wip, dwip); + for(size_t j = 0; j < grid_nodes.size(); ++j) + { + velocity_grad += grid[grid_nodes[j]]->v * dwip[j].transpose(); + } - p->update_deform_grad(velocity_grad, timestep_s); // eq 181 course notes + p->update_deform_grad(velocity_grad, settings.timestep_s); // eq 181 course notes + }, + MPM_MIN_PARALLEL); - } // end loop particles - -} // end update deformation gradient +} // end update deformation gradient // // Grid to Particle velocity and affine matrix // -void MPM::g2p_velocity(std::vector &particles, const std::vector &grid) +void MPM::g2p_velocity(const Settings& settings, std::vector& particles, const std::vector& grid) { - using namespace Eigen; - - // - // Compute particle velocities and affine matrix - // -#pragma omp parallel for - for(int i=0; iv.fill(0.f); - p->B.fill(0.f); - - std::vector grid_nodes; - std::vector wip; - std::vector dwip; - MPM::make_grid_loop(p, grid_nodes, wip, dwip); - - for(int j=0; jv += wip[j]*grid[grid_nodes[j]]->v; // eq 175 course notes - Eigen::Vector3d xixp = xi_xp(grid_nodes[j], p); - p->B += wip[j]*grid[grid_nodes[j]]->v*xixp.transpose(); // eq 176 course notes - } - - } // end loop particles - -} // end grid to particle velocity + using namespace Eigen; + + // + // Compute particle velocities and affine matrix + // + igl::parallel_for( + particles.size(), + [&](int i) + { + Particle* p = particles[i]; + p->v.fill(0.f); + p->B.fill(0.f); + + std::vector grid_nodes; + std::vector wip; + std::vector dwip; + MPM::make_grid_loop(settings, p, grid_nodes, wip, dwip); + + for(size_t j = 0; j < grid_nodes.size(); ++j) + { + p->v += wip[j] * grid[grid_nodes[j]]->v; // eq 175 course notes + Eigen::Vector3d xixp = xi_xp(settings, grid_nodes[j], p); + p->B += wip[j] * grid[grid_nodes[j]]->v * xixp.transpose(); // eq 176 course notes + } + }, + MPM_MIN_PARALLEL); + +} // end grid to particle velocity // // Particle Advection // -void MPM::particle_advection(std::vector &particles) +void MPM::particle_advection(const Settings& settings, std::vector& particles) { -#pragma omp parallel for - for(int i=0; ix += timestep_s*particles[i]->v; // assumes APIC - } // end loop particles + igl::parallel_for( + particles.size(), + [&](int i) + { + particles[i]->x += settings.timestep_s * particles[i]->v; // assumes APIC + }, + MPM_MIN_PARALLEL); -} // end particle advection +} // end particle advection // // Grid Collision // -void MPM::grid_collision(std::vector &active_grid) +void MPM::grid_collision(const Settings& settings, std::vector& active_grid) { - using namespace Eigen; - -#pragma omp parallel for - for(int i=0; iglobal_idx); - - // Boundary: - for(int j=0; j<3; ++j) - { - if(gridx[j] < 2){ node->v[j] = 0.f; } - if(gridx[j] > gridsize[j]-2){ node->v[j] = 0.f; } - } - - } // end loop active grid - -} // end check grid collision - + using namespace Eigen; + + igl::parallel_for( + active_grid.size(), + [&](int i) + { + GridNode* node = active_grid[i]; + Eigen::Vector3i gridx = grid_pos(settings, node->global_idx); + + // Boundary: + for(int j = 0; j < 3; ++j) + { + if(gridx[j] < 2) + { + node->v[j] = 0.f; + } + if(gridx[j] > settings.gridsize[j] - 2) + { + node->v[j] = 0.f; + } + } + }, + MPM_MIN_PARALLEL); + +} // end check grid collision // @@ -295,65 +309,77 @@ void MPM::grid_collision(std::vector &active_grid) // -void MPM::make_grid_loop(const Particle *p, std::vector &grid_idx, std::vector &Wip, std::vector &dWip) +void MPM::make_grid_loop(const Settings& settings, + const Particle* p, + std::vector& grid_idx, + std::vector& Wip, + std::vector& dWip) { - grid_idx.clear(); Wip.clear(); dWip.clear(); - grid_idx.reserve(64); - Wip.reserve(64); - dWip.reserve(64); - - // Our smoothing kernel has a radius of 2: 64 nodes 3D - Eigen::Vector3i normed_x = Eigen::Vector3i(p->x[0]/cellsize[0], p->x[1]/cellsize[1], p->x[2]/cellsize[2]); - Eigen::Vector3i ll_idx = normed_x - Eigen::Vector3i(1,1,1); // get the lowest-left index - - for(int x=ll_idx[0]; x<(ll_idx[0]+4); ++x) - { - double wx_x = p->x[0]/cellsize[0]-double(x); - double wip_x = qspline(wx_x); - double dwip_x = d_qspline(wx_x); - - for(int y=ll_idx[1]; y<(ll_idx[1]+4); ++y) - { - double wy_x = p->x[1]/cellsize[1]-double(y); - double wip_y = qspline(wy_x); - double dwip_y = d_qspline(wy_x); - - for(int z=ll_idx[2]; z<(ll_idx[2]+4); ++z) - { - double wz_x = p->x[2]/cellsize[2]-double(z); - double wip_z = qspline(wz_x); - double dwip_z = d_qspline(wz_x); - - // Compute grid index and interpolation weights - int g_idx = grid_index(Eigen::Vector3i(x,y,z)); - double g_wip = wip_x*wip_y*wip_z; - Eigen::Vector3d g_dwip( - (1.0/cellsize[0])*dwip_x*wip_y*wip_z, - wip_x*(1.0/cellsize[1])*dwip_y*wip_z, - wip_x*wip_y*(1.0/cellsize[2])*dwip_z - ); // eq (after) 124 course notes - - // Check that the computed node is within the grid boundaries. - // Near the edges it may not be (if something went wrong with collision). - if(g_idx < 0 || g_idx > gridsize[0]*gridsize[1]*gridsize[2]-1 || g_wip < 0.0){ continue; } - - // All good, add it to the vectors - grid_idx.push_back(g_idx); - Wip.push_back(g_wip); - dWip.push_back(g_dwip); - - } // end loop z - } // end loop y - } // end loop x - -} // end make grid loop - - -Eigen::Vector3d MPM::xi_xp(const int grid_idx, const Particle *p) + grid_idx.clear(); + Wip.clear(); + dWip.clear(); + grid_idx.reserve(64); + Wip.reserve(64); + dWip.reserve(64); + + // Our smoothing kernel has a radius of 2: 64 nodes 3D + Eigen::Vector3i normed_x = Eigen::Vector3i(p->x[0] / settings.cellsize[0], + p->x[1] / settings.cellsize[1], + p->x[2] / settings.cellsize[2]); + Eigen::Vector3i ll_idx = normed_x - Eigen::Vector3i(1, 1, 1); // get the lowest-left index + + for(int x = ll_idx[0]; x < (ll_idx[0] + 4); ++x) + { + double wx_x = p->x[0] / settings.cellsize[0] - double(x); + double wip_x = qspline(wx_x); + double dwip_x = d_qspline(wx_x); + + for(int y = ll_idx[1]; y < (ll_idx[1] + 4); ++y) + { + double wy_x = p->x[1] / settings.cellsize[1] - double(y); + double wip_y = qspline(wy_x); + double dwip_y = d_qspline(wy_x); + + for(int z = ll_idx[2]; z < (ll_idx[2] + 4); ++z) + { + double wz_x = p->x[2] / settings.cellsize[2] - double(z); + double wip_z = qspline(wz_x); + double dwip_z = d_qspline(wz_x); + + // Compute grid index and interpolation weights + int g_idx = grid_index(settings, Eigen::Vector3i(x, y, z)); + double g_wip = wip_x * wip_y * wip_z; + Eigen::Vector3d g_dwip((1.0 / settings.cellsize[0]) * dwip_x * wip_y * wip_z, + wip_x * (1.0 / settings.cellsize[1]) * dwip_y * wip_z, + wip_x * wip_y * (1.0 / settings.cellsize[2]) * dwip_z); // eq (after) 124 course notes + + // Check that the computed node is within the grid boundaries. + // Near the edges it may not be (if something went wrong with collision). + if(g_idx < 0 || g_idx > settings.gridsize[0] * settings.gridsize[1] * settings.gridsize[2] - 1 + || g_wip < 0.0) + { + continue; + } + + // All good, add it to the vectors + grid_idx.push_back(g_idx); + Wip.push_back(g_wip); + dWip.push_back(g_dwip); + + } // end loop z + } // end loop y + } // end loop x + +} // end make grid loop + + +Eigen::Vector3d MPM::xi_xp(const Settings& settings, const int grid_idx, const Particle* p) { - Eigen::Vector3i gridx = grid_pos(grid_idx); - Eigen::Vector3d xixp(double(gridx[0])*cellsize[0]-p->x[0], double(gridx[1])*cellsize[1]-p->x[1], double(gridx[2])*cellsize[2]-p->x[2]); - return xixp; -} // end compute world coord diff - -} // ns mpm + Eigen::Vector3i gridx = grid_pos(settings, grid_idx); + Eigen::Vector3d xixp(double(gridx[0]) * settings.cellsize[0] - p->x[0], + double(gridx[1]) * settings.cellsize[1] - p->x[1], + double(gridx[2]) * settings.cellsize[2] - p->x[2]); + return xixp; +} // end compute world coord diff + +} // namespace mpm diff --git a/src/MPM.hpp b/src/MPM.hpp index 1a49305..6409b8b 100644 --- a/src/MPM.hpp +++ b/src/MPM.hpp @@ -1,16 +1,16 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER @@ -21,7 +21,10 @@ #define MPM_BASE_HPP 1 #include "Particle.hpp" -#include +#include +#include + +#define MPM_MIN_PARALLEL 100 // min count to parallelize namespace mpm { @@ -29,18 +32,25 @@ namespace mpm // Grid node base class class GridNode { -public: - GridNode(int global_index) : v(0,0,0), m(0), global_idx(global_index), active_idx(0) {} - virtual ~GridNode() {} - - double m; // mass - Eigen::Vector3d v; // velocity - unsigned int global_idx; // index into the global buffer - unsigned int active_idx; // index into the active buffer (set by MPM::p2g_mass) - - std::vector particles; // particles within interp radius (set by MPM::p2g_mass) - std::vector wip; // corresponding particle weight (set by MPM::p2g_mass) - std::vector dwip; // corresponding particle deriv. weight (set by MPM::p2g_mass) + public: + GridNode(int global_index) + : v(0, 0, 0) + , m(0) + , global_idx(global_index) + , active_idx(0) + { + } + virtual ~GridNode() {} + + double m; // mass + Eigen::Vector3d v; // velocity + unsigned int global_idx; // index into the global buffer + unsigned int active_idx; // index into the active buffer (set by MPM::p2g_mass) + std::mutex mutex; // for atomic updates + + std::vector particles; // particles within interp radius (set by MPM::p2g_mass) + std::vector wip; // corresponding particle weight (set by MPM::p2g_mass) + std::vector dwip; // corresponding particle deriv. weight (set by MPM::p2g_mass) }; @@ -48,53 +58,71 @@ class GridNode // for steps 1,2,7,8 of the SIGGRAPH course notes. class MPM { -public: - static Eigen::Vector3d cellsize; - static Eigen::Vector3i gridsize; - static double timestep_s; - - // - // Particle to Grid functions - // - - static void clear_grid(std::vector &active_grid); - static void p2g_mass(std::vector &particles, std::vector &grid, std::vector &active_grid); - static void init_volumes(std::vector &particles, const std::vector &grid); - static void p2g_velocity(const std::vector &particles, std::vector &active_grid); - - // - // Grid to Particle functions - // - - static void g2p_deformation(std::vector &particles, const std::vector &grid); - static void g2p_velocity(std::vector &particles, const std::vector &grid); - static void particle_advection(std::vector &particles); - - // - // Collision - // - - static void grid_collision(std::vector &active_grid); - - // - // Interpolation functions (see Interp.hpp) - // - - // Pass in (normalized) grid index and returns its index in the grid vector. - static int grid_index(const Eigen::Vector3i &pos); - - // Inverse of grid_index (position in grid, NOT world space!) - static Eigen::Vector3i grid_pos(int grid_idx); - - // Makes an easy-to-iterate list of grid nodes that are within interpolation range of a particle. - static void make_grid_loop(const Particle *p, std::vector &grid_idx, std::vector &Wip, std::vector &dWip); - - // Computes grid position minus particle position (world space) for APIC - static Eigen::Vector3d xi_xp(const int grid_idx, const Particle *p); - -}; // end class MPM - -} // end namespace mpm + public: + struct Settings + { + Eigen::Vector3d cellsize = Eigen::Vector3d(0.25, 0.25, 0.25); + Eigen::Vector3i gridsize = Eigen::Vector3i(32, 32, 32); + double timestep_s = 0.04; + }; + + // + // Particle to Grid functions + // + + static void clear_grid(std::vector& active_grid); + static void p2g_mass(const Settings& settings, + std::vector& particles, + std::vector& grid, + std::vector& active_grid); + static void init_volumes(const Settings& settings, + std::vector& particles, + const std::vector& grid); + static void p2g_velocity(const Settings& settings, + const std::vector& particles, + std::vector& active_grid); + + // + // Grid to Particle functions + // + + static void g2p_deformation(const Settings& settings, + std::vector& particles, + const std::vector& grid); + static void g2p_velocity(const Settings& settings, + std::vector& particles, + const std::vector& grid); + static void particle_advection(const Settings& settings, std::vector& particles); + + // + // Collision + // + + static void grid_collision(const Settings& settings, std::vector& active_grid); + + // + // Interpolation functions (see Interp.hpp) + // + + // Pass in (normalized) grid index and returns its index in the grid vector. + static int grid_index(const Settings& settings, const Eigen::Vector3i& pos); + + // Inverse of grid_index (position in grid, NOT world space!) + static Eigen::Vector3i grid_pos(const Settings& settings, int grid_idx); + + // Makes an easy-to-iterate list of grid nodes that are within interpolation range of a particle. + static void make_grid_loop(const Settings& settings, + const Particle* p, + std::vector& grid_idx, + std::vector& Wip, + std::vector& dWip); + + // Computes grid position minus particle position (world space) for APIC + static Eigen::Vector3d xi_xp(const Settings& settings, const int grid_idx, const Particle* p); + +}; // end class MPM + +} // end namespace mpm #endif diff --git a/src/Particle.hpp b/src/Particle.hpp index 6be8521..a559a07 100644 --- a/src/Particle.hpp +++ b/src/Particle.hpp @@ -1,16 +1,16 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER @@ -20,121 +20,142 @@ #ifndef MPM_PARTICLE_HPP #define MPM_PARTICLE_HPP 1 -#include -#include -#include #include "Interp.hpp" +#include +#include namespace mpm { // Projection, Singular Values, SVD's U, SVD's V transpose -static inline void projected_svd(const Eigen::Matrix3d &F, Eigen::Vector3d &S, Eigen::Matrix3d &U, Eigen::Matrix3d &Vt) +static inline void projected_svd(const Eigen::Matrix3d& F, Eigen::Vector3d& S, Eigen::Matrix3d& U, Eigen::Matrix3d& Vt) { - using namespace Eigen; - JacobiSVD svd(F, ComputeFullU | ComputeFullV); - S = svd.singularValues(); - U = svd.matrixU(); - Vt = svd.matrixV().transpose(); - Matrix3d J = Matrix3d::Identity(); J(2,2)=-1.0; - // Check for inversion - if(U.determinant() < 0.0) { U = U * J; S[2] *= -1.0; } - if(Vt.determinant() < 0.0) { Vt = J * Vt; S[2] *= -1.0; } - -} // end projected svd + using namespace Eigen; + JacobiSVD svd(F, ComputeFullU | ComputeFullV); + S = svd.singularValues(); + U = svd.matrixU(); + Vt = svd.matrixV().transpose(); + Matrix3d J = Matrix3d::Identity(); + J(2, 2) = -1.0; + // Check for inversion + if(U.determinant() < 0.0) + { + U = U * J; + S[2] *= -1.0; + } + if(Vt.determinant() < 0.0) + { + Vt = J * Vt; + S[2] *= -1.0; + } + +} // end projected svd // Particle (material point) class class Particle { -public: - Particle() : m(0.1), vol(0.0), v(0,0,0), x(0,0,0) - { - Fe = Eigen::Matrix3d::Identity(); - tempP = Eigen::Matrix3d::Identity(); - B.setZero(); - D.setZero(); - } - - virtual Eigen::Matrix3d get_piola_stress(Eigen::Matrix3d &currF) const = 0; - virtual double get_energy_density(Eigen::Matrix3d &newF) const = 0; - virtual Eigen::Matrix3d get_deform_grad(){ return Fe; } - virtual void update_deform_grad(Eigen::Matrix3d velocity_grad, double timestep_s) - { - Fe = (Eigen::Matrix3d::Identity() + timestep_s*velocity_grad) * Fe; - } - - double m; // mass - double vol; // rest volume, computed at first timestep - Eigen::Vector3d v; // velocity - Eigen::Vector3d x; // location - Eigen::Matrix3d B; // affine matrix (eq 176 course notes) - Eigen::Matrix3d D; // helper matrix, set in p_to_g mass (eq 174 course notes) - - Eigen::Matrix3d tempP; // temporary piola stress tensor computed by solver - Eigen::Matrix3d Fe; // elastic deformation gradient + public: + Particle() + : m(0.1) + , vol(0.0) + , v(0, 0, 0) + , x(0, 0, 0) + { + Fe = Eigen::Matrix3d::Identity(); + tempP = Eigen::Matrix3d::Identity(); + B.setZero(); + D.setZero(); + } + + virtual Eigen::Matrix3d get_piola_stress(Eigen::Matrix3d& currF) const = 0; + virtual double get_energy_density(Eigen::Matrix3d& newF) const = 0; + virtual Eigen::Matrix3d get_deform_grad() { return Fe; } + virtual void update_deform_grad(Eigen::Matrix3d velocity_grad, double timestep_s) + { + Fe = (Eigen::Matrix3d::Identity() + timestep_s * velocity_grad) * Fe; + } + + double m; // mass + double vol; // rest volume, computed at first timestep + Eigen::Vector3d v; // velocity + Eigen::Vector3d x; // location + Eigen::Matrix3d B; // affine matrix (eq 176 course notes) + Eigen::Matrix3d D; // helper matrix, set in p_to_g mass (eq 174 course notes) + + Eigen::Matrix3d tempP; // temporary piola stress tensor computed by solver + Eigen::Matrix3d Fe; // elastic deformation gradient }; // NeoHookean Particle class pNeoHookean : public Particle { -public: - - pNeoHookean() : mu(10), lambda(10) {} - double mu; - double lambda; - - inline Eigen::Matrix3d get_piola_stress( Eigen::Matrix3d &currF ) const - { - // Fix inversions: - Eigen::Matrix3d U, Vt, Ftemp; - Eigen::Vector3d S; - projected_svd(currF, S, U, Vt); - if (S[2] < 0.0) { S[2] *= -1.0; } - Ftemp = U * S.asDiagonal() * Vt; - - // Compute Piola stress tensor: - double J = Ftemp.determinant(); - assert(isreal(J)); - assert(J>0.0); - Eigen::Matrix3d Fit = (Ftemp.inverse()).transpose(); // F^(-T) - Eigen::Matrix3d P = mu*(Ftemp-Fit) + lambda*(log(J)*Fit); - for (int i=0; i<3; ++i) - { - for(int j=0; j<3; ++j) - { - assert(isreal(P(i,j))); - } - } - return P; - - } // end compute piola stress tensor - - inline double get_energy_density(Eigen::Matrix3d &newF) const - { - // Fix inversions: - Eigen::Matrix3d U, Vt, Ftemp; - Eigen::Vector3d S; - projected_svd(newF, S, U, Vt); - if (S[2] < 0.0) { S[2] *= -1.0; } - Ftemp = U * S.asDiagonal() * Vt; - - // Compute energy density: - double J = Ftemp.determinant(); - assert(isreal(J)); - assert(J>0.0); - double t1 = 0.5*mu*((Ftemp.transpose()*Ftemp).trace() - 3.0); - double t2 = -mu * log(J); - double t3 = 0.5*lambda*(log(J)*log(J)); - assert(isreal(t1)); - assert(isreal(t2)); - assert(isreal(t3)); - return (t1+t2+t3); - - } // end compute energy density - -}; // end class neohookean particle - -} // end namespace mpm + public: + pNeoHookean() + : mu(10) + , lambda(10) + { + } + double mu; + double lambda; + + Eigen::Matrix3d get_piola_stress(Eigen::Matrix3d& currF) const + { + // Fix inversions: + Eigen::Matrix3d U, Vt, Ftemp; + Eigen::Vector3d S; + projected_svd(currF, S, U, Vt); + if(S[2] < 0.0) + { + S[2] *= -1.0; + } + Ftemp = U * S.asDiagonal() * Vt; + + // Compute Piola stress tensor: + double J = Ftemp.determinant(); + assert(isreal(J)); + assert(J > 0.0); + Eigen::Matrix3d Fit = (Ftemp.inverse()).transpose(); // F^(-T) + Eigen::Matrix3d P = mu * (Ftemp - Fit) + lambda * (log(J) * Fit); + for(int i = 0; i < 3; ++i) + { + for(int j = 0; j < 3; ++j) + { + assert(isreal(P(i, j))); + } + } + return P; + + } // end compute piola stress tensor + + double get_energy_density(Eigen::Matrix3d& newF) const + { + // Fix inversions: + Eigen::Matrix3d U, Vt, Ftemp; + Eigen::Vector3d S; + projected_svd(newF, S, U, Vt); + if(S[2] < 0.0) + { + S[2] *= -1.0; + } + Ftemp = U * S.asDiagonal() * Vt; + + // Compute energy density: + double J = Ftemp.determinant(); + assert(isreal(J)); + assert(J > 0.0); + double t1 = 0.5 * mu * ((Ftemp.transpose() * Ftemp).trace() - 3.0); + double t2 = -mu * log(J); + double t3 = 0.5 * lambda * (log(J) * log(J)); + assert(isreal(t1)); + assert(isreal(t2)); + assert(isreal(t3)); + return (t1 + t2 + t3); + + } // end compute energy density + +}; // end class neohookean particle + +} // end namespace mpm #endif diff --git a/src/Solver.cpp b/src/Solver.cpp index 13be7ad..1b65356 100644 --- a/src/Solver.cpp +++ b/src/Solver.cpp @@ -1,16 +1,16 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER @@ -18,257 +18,301 @@ // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "Solver.hpp" +#include "LBFGS.hpp" +#include +#include namespace mpm { bool Solver::initialize() { - // Fill a sphere with particles - Eigen::Vector3d sphere_center(3.f,3.f,3.f); - double sphere_rad = 1.0; - double step = 0.2; - double p_mass = 0.01; - for(double x=sphere_center[0]-sphere_rad; xx = pos; - p->m = p_mass; - m_particles.push_back(p); - } - }}} - - // The boundaries are determined by these variables, and assumed positive. - // I.e. the boundary is 0 to gridsize*cellsize. - MPM::cellsize = Eigen::Vector3d(.5,.5,.5); - MPM::gridsize = Eigen::Vector3i(16,16,16); - - for(int i=0; ix = pos; + p->m = p_mass; + particles.push_back(p); + } + } + } + } + + // The boundaries are determined by these variables, and assumed positive. + // I.e. the boundary is 0 to gridsize*cellsize. + mpm_settings.cellsize = Eigen::Vector3d(.5, .5, .5); + mpm_settings.gridsize = Eigen::Vector3i(16, 16, 16); + + for(int i = 0; i < mpm_settings.gridsize[0] * mpm_settings.gridsize[1] * mpm_settings.gridsize[2]; ++i) + { + grid.push_back(new GridNode(i)); + } + + std::cout << "Gridsize: " << mpm_settings.gridsize[0] << ' ' << mpm_settings.gridsize[1] << ' ' + << mpm_settings.gridsize[2] << std::endl; + std::cout << "Cellsize: " << mpm_settings.cellsize[0] << ' ' << mpm_settings.cellsize[1] << ' ' + << mpm_settings.cellsize[2] << std::endl; + std::cout << "Num Particles: " << particles.size() << std::endl; + std::cout << "Num Nodes: " << grid.size() << std::endl; + + + return true; } Solver::~Solver() { - for(int i=0; ix; + int np = particles.size(); + if(np == 0) + { + return; + } + + if(X.rows() < np) + { + X.resize(np, 3); + } + + for(int i = 0; i < np; ++i) + { + X.row(i) = particles[i]->x; + } } -// -// Run a system step -// -bool Solver::step(float screen_dt) +void Solver::step(float screen_dt) { - // Compute timestep - MPM::timestep_s = compute_timestep(screen_dt); - MPM::clear_grid(active_grid); - - // Particle to grid - MPM::p2g_mass(m_particles, m_grid, active_grid); // step 1 - if(elapsed_s <= 0.f){ MPM::init_volumes(m_particles, m_grid); } - MPM::p2g_velocity(m_particles, active_grid); // step 1 + // Compute timestep + mpm_settings.timestep_s = compute_timestep(screen_dt); + MPM::clear_grid(active_grid); - // Steps 4, 5 in course notes -// explicit_solve(); - implicit_solve(); - MPM::grid_collision(active_grid); + // Particle to grid + MPM::p2g_mass(mpm_settings, particles, grid, active_grid); // step 1 + if(elapsed_s <= 0.f) + { + MPM::init_volumes(mpm_settings, particles, grid); + } + MPM::p2g_velocity(mpm_settings, particles, active_grid); // step 1 - // Grid to particle, and advection - // Steps 6, 7, 8 in course notes - MPM::g2p_deformation(m_particles, m_grid); - MPM::g2p_velocity(m_particles, m_grid); - MPM::particle_advection(m_particles); + // Steps 4, 5 in course notes + // explicit_solve(); + implicit_solve(); + MPM::grid_collision(mpm_settings, active_grid); - elapsed_s += MPM::timestep_s; + // Grid to particle, and advection + // Steps 6, 7, 8 in course notes + MPM::g2p_deformation(mpm_settings, particles, grid); + MPM::g2p_velocity(mpm_settings, particles, grid); + MPM::particle_advection(mpm_settings, particles); - // Done - return true; + elapsed_s += mpm_settings.timestep_s; -} // end step +} // end step void Solver::explicit_solve() { -#pragma omp parallel for - for(int i=0; iget_deform_grad(); - p->tempP = p->vol * p->get_piola_stress(p_F) * p->Fe.transpose(); - } // end loop particles - - // Velocity update -#pragma omp parallel for - for(int i=0; iparticles.size(); ++j){ - Particle *p = g->particles[j]; - Eigen::Matrix3d p_F = p->get_deform_grad(); - f -= p->tempP*g->dwip[j]; - } - - g->v += (MPM::timestep_s*f)/g->m + MPM::timestep_s*gravity; - - } // end for all active grid - -} // end compute forces - - -double Solver::gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad) + igl::parallel_for( + particles.size(), + [&](int i) + { + Particle* p = particles[i]; + Eigen::Matrix3d p_F = p->get_deform_grad(); + p->tempP = p->vol * p->get_piola_stress(p_F) * p->Fe.transpose(); + }, + MPM_MIN_PARALLEL); + + + // Velocity update + igl::parallel_for( + active_grid.size(), + [&](int i) + { + GridNode* g = active_grid[i]; + Eigen::Vector3d f(0, 0, 0); + for(size_t j = 0; j < g->particles.size(); ++j) + { + Particle* p = g->particles[j]; + Eigen::Matrix3d p_F = p->get_deform_grad(); + f -= p->tempP * g->dwip[j]; + } + + g->v += (mpm_settings.timestep_s * f) / g->m + mpm_settings.timestep_s * gravity; + }, + MPM_MIN_PARALLEL); + +} // end compute forces + + +double Solver::gradient(const Eigen::VectorXd& v, Eigen::VectorXd& grad) { - using namespace Eigen; - double tot_energy = 0.0; - - // - // Compute Particle Deformation Gradients for new grid velocities - // -#pragma omp parallel for reduction (+:tot_energy) - for(int i=0; iget_deform_grad(); - - Matrix3d vel_grad; vel_grad.fill(0.f); - - std::vector grid_nodes; - std::vector wip; - std::vector dwip; - MPM::make_grid_loop(p, grid_nodes, wip, dwip); - for(int j=0; jactive_idx; - Eigen::Vector3d curr_v(v[idx*3+0], v[idx*3+1], v[idx*3+2]); - vel_grad += (MPM::timestep_s*curr_v) * dwip[j].transpose(); - } - - Eigen::Matrix3d newF = (Matrix3d::Identity() + vel_grad) * p_F; // eq 193 course notes - p->tempP = p->vol * p->get_piola_stress(newF) * p_F.transpose(); - - // Update the total energy - double e = p->get_energy_density(newF)*p->vol; - tot_energy = tot_energy + e; - - } // end loop particles - - - // - // Compute energy gradient - // -#pragma omp parallel for reduction (+:tot_energy) - for(int i=0; im * (curr_v - node->v); - Eigen::Vector3d energy_grad(0,0,0); - - for(int j=0; jparticles.size(); ++j){ - energy_grad += node->particles[j]->tempP * node->dwip[j] * MPM::timestep_s; - } - - if (grad.rows() == v.rows()) - { - grad.segment(i*3, 3) = momentum_grad + energy_grad; - } - tot_energy = tot_energy + 0.5 * node->m * (curr_v - node->v).squaredNorm(); - - } // end loop grid - - return tot_energy; + using namespace Eigen; + VectorXd energy_sum = VectorXd::Zero(std::max(particles.size(), active_grid.size())); + + // + // Compute Particle Deformation Gradients for new grid velocities + // + igl::parallel_for( + particles.size(), + [&](int i) + { + Particle* p = particles[i]; + Matrix3d p_F = p->get_deform_grad(); + + Matrix3d vel_grad; + vel_grad.fill(0.f); + + std::vector grid_nodes; + std::vector wip; + std::vector dwip; + MPM::make_grid_loop(mpm_settings, p, grid_nodes, wip, dwip); + for(size_t j = 0; j < grid_nodes.size(); ++j) + { + int idx = grid[grid_nodes[j]]->active_idx; + Vector3d curr_v(v[idx * 3 + 0], v[idx * 3 + 1], v[idx * 3 + 2]); + vel_grad += (mpm_settings.timestep_s * curr_v) * dwip[j].transpose(); + } + + Matrix3d newF = (Matrix3d::Identity() + vel_grad) * p_F; // eq 193 course notes + p->tempP = p->vol * p->get_piola_stress(newF) * p_F.transpose(); + + // Update the total energy + energy_sum[i] = p->get_energy_density(newF) * p->vol; + }, + MPM_MIN_PARALLEL); + + + // + // Compute energy gradient + // + igl::parallel_for( + active_grid.size(), + [&](int i) + { + GridNode* node = active_grid[i]; + Vector3d curr_v(v[i * 3 + 0], v[i * 3 + 1], v[i * 3 + 2]); + + // Only compute gradient if requested (grad is sized) + if(grad.rows() == v.rows()) + { + Vector3d momentum_grad = node->m * (curr_v - node->v); + Vector3d energy_grad(0, 0, 0); + for(size_t j = 0; j < node->particles.size(); ++j) + { + energy_grad += node->particles[j]->tempP * node->dwip[j] * mpm_settings.timestep_s; + } + + grad.segment(i * 3, 3) = momentum_grad + energy_grad; + } + energy_sum[i] += 0.5 * node->m * (curr_v - node->v).squaredNorm(); + }, + MPM_MIN_PARALLEL); + + return energy_sum.sum(); } void Solver::implicit_solve() { - using namespace Eigen; - - // Initial guess of grid velocities - VectorXd v(active_grid.size()*3); -#pragma omp parallel for - for(int i=0; iv[j]; - } - } - - // Minimize - optimizer.gradient = [&](const VectorXd &x, VectorXd &g)->double - { - double objective = gradient(x, g); - return objective; - }; - - optimizer.minimize(v); - - // Copy solver results back to grid -#pragma omp parallel for - for(int i=0; iv[j] = v[i*3+j]; - } - active_grid[i]->v += MPM::timestep_s*gravity; // explicitly add gravity - } - -} // end implicit solve + using namespace Eigen; + + // Initial guess of grid velocities + VectorXd v = VectorXd::Zero(active_grid.size() * 3); + igl::parallel_for( + active_grid.size(), + [&](int i) + { + for(int j = 0; j < 3; ++j) + { + v[i * 3 + j] = active_grid[i]->v[j]; + } + }, + MPM_MIN_PARALLEL); + + // Minimize + mcl::LBFGS optimizer; + + optimizer.gradient = [&](const VectorXd& x, VectorXd& g) -> double + { + double objective = gradient(x, g); + return objective; + }; + + optimizer.minimize(v); + + // Copy solver results back to grid + igl::parallel_for( + active_grid.size(), + [&](int i) + { + for(int j = 0; j < 3; ++j) + { + active_grid[i]->v[j] = v[i * 3 + j]; + } + active_grid[i]->v += mpm_settings.timestep_s * gravity; // explicitly add gravity + }, + MPM_MIN_PARALLEL); + +} // end implicit solve double Solver::compute_timestep(double screen_dt) { - double max_v = 0.0; - double max_timestep = 0.01; - -#pragma omp parallel for reduction(max:max_v) - for(int i=0; iv.squaredNorm(); - if(curr_v > max_v){ max_v = curr_v; } - } max_v = sqrtf(max_v); - - if(max_v > 1e-8) - { - double min_cellsize = std::min(MPM::cellsize[0], std::min(MPM::cellsize[1], MPM::cellsize[2])); - double ts_adjust = 0.04; - double dt = ts_adjust * min_cellsize / sqrt(max_v); - if(dt > max_timestep){ return max_timestep; } - if(dt > screen_dt){ return screen_dt; } - return dt; - } - - return max_timestep; - -} // end adaptive timestep - -} // ns mpm + using namespace Eigen; + double max_timestep = 0.01; + VectorXd squared_vel = VectorXd::Zero(particles.size()); + + igl::parallel_for( + particles.size(), [&](int i) { squared_vel[i] = particles[i]->v.squaredNorm(); }, MPM_MIN_PARALLEL); + + double max_v = std::sqrt(squared_vel.maxCoeff()); + if(max_v > 1e-8) + { + double min_cellsize = std::min(mpm_settings.cellsize[0], + std::min(mpm_settings.cellsize[1], mpm_settings.cellsize[2])); + double ts_adjust = 0.04; + double dt = ts_adjust * min_cellsize / sqrt(max_v); + if(dt > max_timestep) + { + return max_timestep; + } + if(dt > screen_dt) + { + return screen_dt; + } + return dt; + } + + return max_timestep; + +} // end adaptive timestep + +} // namespace mpm diff --git a/src/Solver.hpp b/src/Solver.hpp index 53ebdbc..b63994f 100644 --- a/src/Solver.hpp +++ b/src/Solver.hpp @@ -1,16 +1,16 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER @@ -21,62 +21,46 @@ #define MPM_SOLVER_HPP 1 #include "MPM.hpp" -#include "LBFGS.hpp" namespace mpm { class Solver { -public: + public: + Solver() + : gravity(0.0, -9.8, 0.0) + , elapsed_s(0) + { + } + ~Solver(); - Solver() : elapsed_s(0) {} - ~Solver(); + std::vector grid; + std::vector particles; + std::vector active_grid; // resized each time step + MPM::Settings mpm_settings; + const Eigen::Vector3d gravity; + double elapsed_s; - // - // Data - // + // Initialize the system + bool initialize(); - std::vector m_grid; - std::vector active_grid; // resized each time step - std::vector m_particles; - void get_vertices(Eigen::MatrixXd &X) const; - mcl::LBFGS optimizer; + // Run a timestep. + void step(float screen_dt); - // - // Settings - // + // Gradient of the objective. + double gradient(const Eigen::VectorXd& v, Eigen::VectorXd& grad); - double elapsed_s; - const Eigen::Vector3d gravity = Eigen::Vector3d(0.0,-9.8,0.0); + // Retuns current particle positions + void get_vertices(Eigen::MatrixXd& X) const; - // - // Called by the gui - // + private: + double compute_timestep(double screen_dt); + void explicit_solve(); + void implicit_solve(); - // Initialize the system - bool initialize(); +}; // end class solver - // Run a timestep. - // Returns true on success. - bool step(float screen_dt); - - // Gradient of the objective. - double gradient(const Eigen::VectorXd &v, Eigen::VectorXd &grad); - -private: - - double compute_timestep(double screen_dt); - - // - // Integration - // - - void explicit_solve(); - void implicit_solve(); - -}; // end class solver - -} // end namespace mpm +} // end namespace mpm #endif diff --git a/test/solver.cpp b/test/solver.cpp index aa06e89..196c6fc 100644 --- a/test/solver.cpp +++ b/test/solver.cpp @@ -2,11 +2,11 @@ #include "Solver.hpp" // placeholder -int main(int argc, char *argv[]) +int main(int argc, char* argv[]) { - mpm::Solver solver; + mpm::Solver solver; - (void)(argc); - (void)(argv); - return EXIT_SUCCESS; + (void)(argc); + (void)(argv); + return EXIT_SUCCESS; } \ No newline at end of file diff --git a/test/sphere.cpp b/test/sphere.cpp index a2c001b..acc70ea 100644 --- a/test/sphere.cpp +++ b/test/sphere.cpp @@ -1,16 +1,16 @@ // Copyright (c) 2016 Matt Overby -// +// // MPM-OPTIMIZATION Uses the BSD 2-Clause License (http://www.opensource.org/licenses/BSD-2-Clause) // Redistribution and use in source and binary forms, with or without modification, are -// permitted provided that the following conditions are met: +// permitted provided that the following conditions are met: // 1. Redistributions of source code must retain the above copyright notice, this list of -// conditions and the following disclaimer. +// conditions and the following disclaimer. // 2. Redistributions in binary form must reproduce the above copyright notice, this list // of conditions and the following disclaimer in the documentation and/or other materials -// provided with the distribution. +// provided with the distribution. // THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA, DULUTH OR CONTRIBUTORS BE +// ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF MINNESOTA OR CONTRIBUTORS BE // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, // OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER @@ -20,48 +20,47 @@ #include "Solver.hpp" #include -int main(int argc, char *argv[]) +int main(int argc, char* argv[]) { - using namespace mpm; - using namespace Eigen; + using namespace mpm; + using namespace Eigen; - Solver solver; - solver.initialize(); + Solver solver; + solver.initialize(); - MatrixXd V; - solver.get_vertices(V); + MatrixXd V; + solver.get_vertices(V); - // Also add ground plane - MatrixXd gV(4,3); - gV.row(0) = RowVector3d(0,0.75,0); - gV.row(1) = RowVector3d(6,0.75,0); - gV.row(2) = RowVector3d(6,0.75,6); - gV.row(3) = RowVector3d(0,0.75,6); - MatrixXi gF(2,3); - gF.row(0) = RowVector3i(0,2,1); - gF.row(1) = RowVector3i(0,3,2); + // Also add ground plane + MatrixXd gV(4, 3); + gV.row(0) = RowVector3d(0, 0.75, 0); + gV.row(1) = RowVector3d(6, 0.75, 0); + gV.row(2) = RowVector3d(6, 0.75, 6); + gV.row(3) = RowVector3d(0, 0.75, 6); + MatrixXi gF(2, 3); + gF.row(0) = RowVector3i(0, 2, 1); + gF.row(1) = RowVector3i(0, 3, 2); - std::cout << "Press A to toggle animation" << std::endl; + std::cout << "Press A to toggle animation" << std::endl; - // Create viewer - igl::opengl::glfw::Viewer viewer; - viewer.core().is_animating = false; - viewer.data().set_mesh(gV, gF); - viewer.data().add_points(V, Eigen::RowVector3d(1,0,0)); - viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer&)->bool - { - if (viewer.core().is_animating) - { - solver.step(0.04); - solver.get_vertices(V); - viewer.data().clear(); - viewer.data().set_mesh(gV, gF); - viewer.data().add_points(V, Eigen::RowVector3d(1,0,0)); - } - return false; - }; - viewer.launch(); + // Create viewer + igl::opengl::glfw::Viewer viewer; + viewer.core().is_animating = false; + viewer.data().set_mesh(gV, gF); + viewer.data().add_points(V, Eigen::RowVector3d(1, 0, 0)); + viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer&) -> bool + { + if(viewer.core().is_animating) + { + solver.step(0.04); + solver.get_vertices(V); + viewer.data().clear(); + viewer.data().set_mesh(gV, gF); + viewer.data().add_points(V, Eigen::RowVector3d(1, 0, 0)); + } + return false; + }; + viewer.launch(); - return EXIT_SUCCESS; + return EXIT_SUCCESS; } -