First Commit

.github/workflows/build_and_test.yml

@@ -20,9 +20,11 @@ jobs:
         sudo apt-get install -y build-essential cmake
 
     - name: build
+      working-directory: ${{ github.workspace }} 
       run: |
-        cmake -DCLOTH_DEBUG=ON .
+        cmake -DCLOTH_TESTING_ONLY=ON .
         make -j
 
     - name: test
+      working-directory: ${{ github.workspace }} 
       run: make test

.gitignore

@@ -5,7 +5,6 @@
 *.slo
 *.lo
 *.o
-*.obj
 
 # Precompiled Headers
 *.gch

CMakeLists.txt

@@ -0,0 +1,46 @@
+# Copyright Matt Overby 2021.
+# Distributed under the MIT License.
+cmake_minimum_required(VERSION 3.5)
+
+project(cloth_sim_tutorial LANGUAGES CXX)
+
+option(CLOTH_DEBUG "run with debug options" OFF)
+option(CLOTH_TESTING_ONLY "only compile unit tests" OFF)
+
+set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+add_definitions(-DCLOTH_ROOT_DIR="${CMAKE_CURRENT_SOURCE_DIR}")
+
+# Change compile flags for debug
+if (CLOTH_DEBUG)
+	add_definitions(-DCLOTH_DEBUG)
+	set(CMAKE_BUILD_TYPE Debug)
+	if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
+		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wundef -Wextra -Wno-long-long -Wlogical-op -Wshadow -Wmaybe-uninitialized")
+	endif()
+else()
+	set(CMAKE_BUILD_TYPE Release)
+endif()
+
+# Libigl also includes Eigen, etc
+include(libigl)
+
+# Gather src
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src)
+file(GLOB SRC
+	${CMAKE_CURRENT_SOURCE_DIR}/src/*.hpp
+	${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp)
+
+# Create binary
+if (NOT CLOTH_TESTING_ONLY)
+	igl_include(glfw)
+	add_executable(${PROJECT_NAME} ${SRC})
+	target_link_libraries(${PROJECT_NAME} PUBLIC igl::glfw)
+endif()
+
+# Unit testing (work in progress)
+enable_testing()
+add_executable(cloth_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test.cpp)
+target_link_libraries(cloth_test PUBLIC igl::core)
+add_test(NAME TestClothSolver COMMAND cloth_test)

README.md

@@ -1,2 +1,91 @@
-# optimization-cloth-tutorial
-A tutorial to implicit cloth simulation with numerical optimization 
+# Implicit Cloth Simulation with Gradient Descent
+
+A gentle introduction to optimization-based implicit time integration for [CSCI 5611](https://umtc.catalog.prod.coursedog.com/courses/8103481).
+
+_This is a work in progress. I released this around 2021 and it hasn't seen many updates. Please let me know if you have suggestions: mattoverby@gmail.com_
+
+## Introduction
+
+Cloth simulation is standard coursework for undergraduate computer graphics.
+Many assignments focus on explicit methods like [forward Euler](https://en.wikipedia.org/wiki/Euler_method).
+While efficient and simple, explicit integrators are notoriously unstable for large time steps and high stiffness.
+Implicit methods, e.g., [backward Euler](https://en.wikipedia.org/wiki/Backward_Euler_method), offer a solution to that problem.
+The standard [Baraff & Witkin approach](https://doi.org/10.1145/280814.280821) (B&W) forms a Taylor approximation to the Eqs. of motion, i.e., Newton's method.
+This results in a block-sparse linear system for which the velocity deltas are solved via [conjugate gradient](https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf).
+
+The B&W method is conceptually straight forward. However, there are numerous gotchas not apparent in the original references.
+In particular, working out the Jacobians of some models can be nontrivial. Even when correct, the derivatives can still produce hard-to-debug instabilities.
+The [dynamic deformables course notes](https://doi.org/10.1145/3388769.3407490) describes this well, and I highly recommend reviewing the later chapters on cloth simulation for a thorough review of the topic.
+Nonetheless, implicit time integration remains a valuable and sometimes necessary tool for modern graphics applications.
+
+Recently, optimization-based methods have grown in popularity. This approach formulates the Eqs. of motion as an [iteratively minimized objective function](https://en.wikipedia.org/wiki/Optimization_problem).
+These methods provide unparalleled robustness and extensibility, e.g., [state of the art solvers](https://doi.org/10.1145/3450626.3459767) for cloth simulation.
+The implementation is often easier to construct piece-by-piece than B&W-style integrators.
+
+Below we describe a simple implicit solver for mass-spring systems using gradient descent.
+This is meant solely as a tutorial to optimization-based time stepping, and not for practical applications.
+For that, you'll want to consider higher-order optimization algorithms or [local-global type solvers](https://github.com/alecjacobson/computer-graphics-mass-spring-systems).
+
+## Implementation
+
+The objective we aim to minimize is the sum of internal elastic energy potentials plus a quadratic penalty of linear momentum:
+
+<img src="https://render.githubusercontent.com/render/math?math=\bar{x} = x %2B hv %2B h^2 M^{-1} f_{ext}">, <br/>
+<img src="https://render.githubusercontent.com/render/math?math=g(x) = \frac{1}{2h^2}\|M^{1/2}(x-\bar{x})\|^2 %2B \sum E(x)">,
+
+for vertex locations *x*, time step (sec) *h*, diagonal mass matrix *M*, velocities *v*, and external forces (gravity, wind, etc...)
+<img src="https://render.githubusercontent.com/render/math?math=f_{ext}">. The term <img src="https://render.githubusercontent.com/render/math?math=\bar{x}"> is the *explicit predictor*
+that is computed at the beginning of the time step. Computing a frame of animation amounts to iteratively minimizing the above objective until some convergence criteria is met
+(e.g., the norm of the gradient is below some threshold).
+
+### Interface
+
+Below we lay out the steps to produce an optimization-based solver for implicit cloth simulation.
+The dependencies will be fetched through a [DownloadProject](https://github.com/Crascit/DownloadProject) CMake script.
+[Eigen](https://eigen.tuxfamily.org) is used for linear algebra, which is an industry standard library for vector math.
+Rendering and mesh I/O is handled with [libigl](https://libigl.github.io) which also contains many excellent routines for mesh processing.
+In libigl, meshes are stored as two dense matrices, an nx3 matrix for vertices and an mx3 matrix for faces.
+Animating the mesh is done through a pre-draw callback in which we call our solver to update the vertex positions.
+It takes a little finesse to improve the rendering, so we'll save that for another day.
+
+To compile, execute the following commands:
+1. mkdir build
+2. cd build
+3. cmake ..
+4. make -j
+
+### Energies
+
+Whenever I write a new simulator, I start with the most basic elastic energy as the deformable primitives: a Hookean spring with stiffness *k* and rest length <img src="https://render.githubusercontent.com/render/math?math=\ell">. The potential energy of a spring is <img src="https://render.githubusercontent.com/render/math?math=(k/2)(\|x_a-x_b\|-\ell)^2">.
+These springs are created from the unique edges of tCe triangle mesh, with bending springs created between adjacent triangles.
+The [ClothMesh class](src/ClothMesh.hpp) generates this data given vertex and face buffers.
+We can also evaluate temporary springs to deal with collisions. One end point is the vertex, and the other is the projection onto the surface it's penetrating.
+
+### Optimization
+
+[Gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) is a simple first-order method that steps along the negative gradient w.r.t. *x* at each iteration: <img src="https://render.githubusercontent.com/render/math?math=p = -\nabla g(x^i)">. We use *i* to denote the iteration and *p* as the descent direction. By all accounts, gradient descent is a poor choice due to its low rate of convergence. You can certainly improve performance with other higher-order optimizers. However, I personally advocate starting with gradient descent when writing new optimization code. The optimizer is so simple it's hard to make a mistake - if there are bugs, they are surely in your energy or gradient calculations. This can be a helpful debugging strategy!
+
+Each iteration, we have to make sure we aren't overshooting the objective and increase the energy. This is accomplished by scaling the descent direction with a scalar *s* using [line search](https://en.wikipedia.org/wiki/Line_search). All together, an iteration of gradient descent involves:
+
+<img src="https://render.githubusercontent.com/render/math?math=p = -\nabla g(x^i)">,<br>
+<img src="https://render.githubusercontent.com/render/math?math=s = linesearch(x^i, p)">,<br>
+<img src="https://render.githubusercontent.com/render/math?math=x^{i %2B 1} = x^i %2B s * p">.
+
+Both energy and gradient use similar calculations so you can save resources (and implementation efforts) by computing them both at same time. I like to use one function for both and skip gradient calculation with a conditional during line search. You can see the evaluation of the objective and gradient in the [Objective class](src/Objective.hpp), with gradient descent in [Solver](src/Solver.hpp).
+
+### Common Problems
+With these pieces you have enough to implement an implicit cloth solver. But there are a number of common problems to be aware of if your solver starts behaving erratically:
+- Using masses that aren't scaled by the area of the mesh
+- The cloth mesh is unrealistically large or small: Each unit of the implementation should correspond to one meter or centimeter. Whatever you choose, be consistent.
+- The time step is too large or small (it should be somewhere between 1/1000 and 1/20)
+- Forgetting to resize or initialize matrix values
+- Dividing by zero and error propogation: A single zero-length spring can produce instabilities, as its result is passed all the way up the food chain to the solver. Using asserts, guards, and std::isfinite any time you do a division will save you a lot of trouble (see the [Objective class](src/Objective.hpp) for an example).
+
+## Next Steps
+
+There are many ways to improve the solver illustrated above:
+- Parallelizing the objective and gradient calculation (I recommend [TBB](https://software.intel.com/en-us/oneapi/onetbb))
+- Collision against triangle meshes with a [signed distance field](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance)
+- Better energy models with triangular elements ([Dynamic Deformables appendix D](https://doi.org/10.1145/3388769.3407490))
+- Preconditioning or accelerated versions of gradient descent
+- Optimizers with better convergence like [L-BFGS](https://en.wikipedia.org/wiki/Limited-memory_BFGS) or [projected Newton](https://doi.org/10.1145/1073368.1073394)

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)