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 2d619d7..7203552 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -12,7 +12,7 @@ jobs: steps: - name: clone - uses: actions/checkout@v4 # Action to check out your repository + uses: actions/checkout@v4 - name: dependencies run: | 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 fd36516..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 @@ -15,8 +22,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. -It is not efficient nor bug free. +My ideas didn't work out, but I put this code up on github instead of letting it collect dust on my hard drive. 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; } -