Update tutorial documentation (#101)

zfergus · web-flow · commit 3af9200c3b01 · 2024-07-12T00:52:38.000-04:00
* Update project to PSD docs

* Update simulation.rst
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -20,10 +20,11 @@
 
 # -- Project information -----------------------------------------------------
 
+from datetime import datetime
 import ipctk
 
 project = "IPC Toolkit"
-copyright = '2020-2023, IPC-Sim Organization; MIT License'
+copyright = f'2020-{datetime.now().year}, IPC-Sim Organization; MIT License'
 author = "Zachary Ferguson"
 version = ipctk.__version__
 
diff --git a/docs/source/cpp-api/utils.rst b/docs/source/cpp-api/utils.rst
@@ -11,4 +11,6 @@ Positive Semi-Definite Projection
 ---------------------------------
 
 .. doxygenfunction:: ipc::project_to_psd
-.. doxygenfunction:: ipc::project_to_pd
+.. doxygenfunction:: ipc::project_to_pd
+
+.. doxygenenum:: ipc::PSDProjectionMethod
diff --git a/docs/source/python-api/utils.rst b/docs/source/python-api/utils.rst
@@ -17,4 +17,5 @@ Positive Semi-Definite Projection
 ---------------------------------
 
 .. autofunction:: ipctk.project_to_psd
-.. autofunction:: ipctk.project_to_pd
+.. autofunction:: ipctk.project_to_pd
+.. autoclass:: ipctk.PSDProjectionMethod
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
@@ -67,4 +67,11 @@ @article{Ferguson2021RigidIPC
 	number    = 4,
 	articleno = 183,
 	note      = {\url{https://ipc-sim.github.io/rigid-ipc}}
-}
+}
+@inproceedings{Chen2024Stabler,
+	title        = {Stabler Neo-Hookean Simulation: Absolute Eigenvalue Filtering for Projected Newton},
+	author       = {Honglin Chen and Hsueh-Ti Derek Liu and David I.W. Levin and Changxi Zheng and Alec Jacobson},
+	year         = 2024,
+	booktitle    = {{ACM} {SIGGRAPH} 2024 Conference Proceedings},
+	note         = {\url{https://www.cs.columbia.edu/cg/abs-psd/}}
+}
diff --git a/docs/source/tutorial/simulation.rst b/docs/source/tutorial/simulation.rst
@@ -30,10 +30,8 @@ From the full (volumetric) mesh vertices and surface edges/faces which index int
             Eigen::MatrixXi edges;
             igl::edges(faces, edges);
 
-            std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(node_positions.rows(), boundary_edges);
-
             ipc::CollisionMesh collision_mesh =
-                CollisionMesh::build_from_full_mesh(full_rest_positions, edges, faces);
+                ipc::CollisionMesh::build_from_full_mesh(full_rest_positions, edges, faces);
 
     .. md-tab-item:: Python
 
@@ -49,7 +47,69 @@ From the full (volumetric) mesh vertices and surface edges/faces which index int
             collision_mesh = ipctk.CollisionMesh.build_from_full_mesh(
                 full_rest_positions, edges, faces)
 
-This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. However, when computing the gradient and Hessian of the potentials, the derivatives will be with respect to the surface DOF. If you want the derivatives with respect to the full mesh DOF, then we need to apply the chain rule. Fortunately, the ``CollisionMesh`` class provides a function to do this (``CollisionMesh::to_full_dof``):
+This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. However, when passing the collision mesh to toolkit functions, the vertices have to be the surface vertices. The ``CollisionMesh`` class provides a function to map the full vertices and velocities to the surface vertices and velocities. These are ``CollisionMesh::vertices(full_vertices)`` and ``CollisionMesh::map_displacements(full_displacements)``, respectively.
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // Convert full vertices to surface vertices
+            Eigen::VectorXd vertices = collision_mesh.vertices(full_vertices);
+
+            // Construct the set of collisions
+            ipc::Collisions collisions;
+            collisions.build(collision_mesh, vertices, dhat);
+
+            // Construct a barrier potential
+            ipc::BarrierPotential B(dhat);
+
+            // Evaluate the potential
+            double b = B(collisions, collision_mesh, vertices);
+
+            // Convert full velocities to surface velocities
+            Eigen::VectorXd velocities = collision_mesh.map_displacements(full_velocities);
+
+            // Construct the set of friction collisions
+            ipc::FrictionCollisions friction_collisions;
+            friction_collisions.build(collision_mesh, vertices, collisions, B, barrier_stiffness, mu);
+
+            // Construct a friction dissipative potential
+            ipc::FrictionPotential D(epsv);
+
+            double d = D(friction_collisions, collision_mesh, velocities);
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            # Convert full vertices to surface vertices
+            vertices = collision_mesh.vertices(full_vertices)
+
+            # Construct the set of collisions
+            collisions = ipctk.Collisions()
+            collisions.build(collision_mesh, vertices, dhat)
+
+            # Construct a barrier potential
+            B = ipctk.BarrierPotential(dhat)
+
+            # Evaluate the potential
+            b = B(collisions, collision_mesh, vertices)
+
+            # Convert full velocities to surface velocities
+            velocities = collision_mesh.map_displacements(full_velocities)
+
+            # Construct the set of friction collisions
+            friction_collisions = ipctk.FrictionCollisions()
+            friction_collisions.build(collision_mesh, vertices, collisions, B, barrier_stiffness, mu)
+
+            # Construct a friction dissipative potential
+            D = ipctk.FrictionPotential(epsv)
+
+            d = D(friction_collisions, collision_mesh, velocities)
+
+When computing the gradient and Hessian of the potentials, the derivatives will be with respect to the surface DOF. If you want the derivatives with respect to the full mesh DOF, then we need to apply the chain rule. Fortunately, the ``CollisionMesh`` class provides a function to do this (``CollisionMesh::to_full_dof``):
 
 .. md-tab-set::
 
@@ -77,6 +137,43 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How
             hess = B.hessian(collision, collision_mesh, vertices)
             hess_full = collision_mesh.to_full_dof(hess)
 
+Codimensional Vertices
+^^^^^^^^^^^^^^^^^^^^^^
+
+In some cases, the collision mesh vertices are not the same as the surface vertices of the volumetric mesh vertices. One such case is when simulating codimensional vertices in conjunction with shell or volumetric meshes. In this case, simply calling ``build_from_full_mesh`` will not work as it will ignore the vertices that are not connected to any boundary edge. Instead, you can build a vector of booleans that indicate which vertices are on the surface and pass it to the ``CollisionMesh`` constructor.
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // codim_vertices is a vector of indices of the codimensional vertices
+            Eigen::VectorXi codim_vertices = ...;
+
+            // is_on_surface is a vector of booleans indicating which vertices are on the surface
+            std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(
+                full_rest_positions.rows(), boundary_edges, codim_vertices);
+
+            // Construct the collision mesh from the is_on_surface vector and full mesh data
+            ipc::CollisionMesh collision_mesh(
+                is_on_surface, full_rest_positions, edges, faces);
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            # codim_vertices is an array of indices of the codimensional vertices
+            codim_vertices = ...
+
+            # is_on_surface is a list of booleans indicating which vertices are on the surface
+            is_on_surface = ipctk.CollisionMesh.construct_is_on_surface(
+                len(full_rest_positions), boundary_edges, codim_vertices)
+
+            # Construct the collision mesh from the is_on_surface vector and full mesh data
+            collision_mesh = ipctk.CollisionMesh(
+                is_on_surface, full_rest_positions, edges, faces)
+
 Nonlinear Bases and Curved Meshes
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -118,24 +215,37 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me
 
             # Collision proxy mesh
             # Load the proxy mesh from a file
-            proxy_mesh = meshio.read("proxy.msh")
+            proxy_mesh = meshio.read("proxy.obj")
             proxy_rest_positions = proxy_mesh.points
             proxy_faces = proxy_mesh.cells_dict["triangle"]
             proxy_edges = igl.edges(proxy_faces)
-            # Or build it from the volumetric mesh
+            # or build it from the volumetric mesh ...
 
             # Linear map from the finite element mesh to the collision proxy
             displacement_map = ... # build or load the displacement map
 
             collision_mesh = CollisionMesh(
                 proxy_rest_positions, proxy_edges, proxy_faces, displacement_map)
 
-We can then map the displacements using ``collision_mesh.map_displacement(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``.
+We can then map the displacements using ``collision_mesh.map_displacement(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/potential gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``.
 
 .. warning::
     The function ``CollisionMesh::vertices(full_positions)`` should not be used in this case because the rest positions used to construct the ``CollisionMesh`` are not the same as the finite element mesh's rest positions. Instead, use ``CollisionMesh::displace_vertices(fe_displacements)`` where ``fe_displacements`` is already the solution of the PDE or can be computed as ``fe_displacements = fe_positions - fe_rest_positions`` from deformed and rest positions.
 
 Positive Semi-Definite Projection
 ---------------------------------
 
-As described by :cite:t:`Li2020IPC`, the Hessian of the potentials can be indefinite. This is problematic when using the Hessian in a Newton step :cite:p:`Li2020IPC`. To remedy this, we can project the Hessian onto the positive semidefinite (PSD) cone. To do this set the optional parameter ``project_hessian_to_psd`` of ``compute_potential_hessian`` to true.
+As described by :cite:t:`Li2020IPC`, the Hessian of the potentials can be indefinite. This is problematic when using the Hessian in a Newton step :cite:p:`Li2020IPC`.
+To remedy this, we can project the Hessian onto the positive semidefinite (PSD) cone. To do this set the optional parameter ``project_hessian_to_psd`` in ``Potential::hessian`` to one of the following.
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        - ``ProjectToPSD::CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
+        - ``ProjectToPSD::ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
+
+    .. md-tab-item:: Python
+
+        - ``ProjectToPSD.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
+        - ``ProjectToPSD.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
diff --git a/python/src/utils/eigen_ext.cpp b/python/src/utils/eigen_ext.cpp
@@ -28,7 +28,8 @@ void define_eigen_ext(py::module_& m)
         "Enumeration of implemented PSD projection methods.")
         .value("NONE", PSDProjectionMethod::NONE, "No PSD projection")
         .value(
-            "CLAMP", PSDProjectionMethod::CLAMP, "Clamp negative eigenvalues")
+            "CLAMP", PSDProjectionMethod::CLAMP,
+            "Clamp negative eigenvalues to zero")
         .value("ABS", PSDProjectionMethod::ABS, "Flip negative eigenvalues")
         .export_values();
 
@@ -47,5 +48,5 @@ void define_eigen_ext(py::module_& m)
         Returns:
             Projected matrix
         )ipc_Qu8mg5v7",
-        py::arg("A"), py::arg("method"));
+        py::arg("A"), py::arg("method") = PSDProjectionMethod::CLAMP);
 }
diff --git a/python/src/utils/logger.cpp b/python/src/utils/logger.cpp
@@ -7,14 +7,16 @@ using namespace ipc;
 
 void define_logger(py::module_& m)
 {
-    py::enum_<spdlog::level::level_enum>(m, "LoggerLevel")
-        .value("trace", spdlog::level::level_enum::trace)
-        .value("debug", spdlog::level::level_enum::debug)
-        .value("info", spdlog::level::level_enum::info)
-        .value("warn", spdlog::level::level_enum::warn)
-        .value("error", spdlog::level::level_enum::err)
-        .value("critical", spdlog::level::level_enum::critical)
-        .value("off", spdlog::level::level_enum::off)
+    py::enum_<spdlog::level::level_enum>(
+        m, "LoggerLevel", "Enumeration of log levels")
+        .value("trace", spdlog::level::level_enum::trace, "Trace level")
+        .value("debug", spdlog::level::level_enum::debug, "Debug level")
+        .value("info", spdlog::level::level_enum::info, "Info level")
+        .value("warn", spdlog::level::level_enum::warn, "Warning level")
+        .value("error", spdlog::level::level_enum::err, "Error level")
+        .value(
+            "critical", spdlog::level::level_enum::critical, "Critical level")
+        .value("off", spdlog::level::level_enum::off, "Off level")
         .export_values();
 
     m.def(
diff --git a/src/ipc/collision_mesh.cpp b/src/ipc/collision_mesh.cpp
@@ -60,7 +60,7 @@ CollisionMesh::CollisionMesh(
 
     const int dim = full_rest_positions.cols();
 
-    // Selection matrix S ∈ ℝ^{collision×full}
+    // Initializes m_select_vertices and m_select_dof
     init_selection_matrices(dim);
 
     if (displacement_map.size() == 0) {
diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp
@@ -326,13 +326,13 @@ class CollisionMesh {
 
     /// @brief Selection matrix S ∈ ℝ^{collision×full} for vertices
     Eigen::SparseMatrix<double> m_select_vertices;
-    /// @brief Selection matrix S ∈ ℝ^{collision×full} for DOF
+    /// @brief Selection matrix S ∈ ℝ^{(dim*collision)×(dim*full)} for DOF
     Eigen::SparseMatrix<double> m_select_dof;
 
     /// @brief Mapping from full displacements to collision displacements
     /// @note this is premultiplied by m_select_vertices
     Eigen::SparseMatrix<double> m_displacement_map;
-    /// @brief Mapping from full displacements to collision displacements
+    /// @brief Mapping from full displacements DOF to collision displacements DOF
     /// @note this is premultiplied by m_select_dof
     Eigen::SparseMatrix<double> m_displacement_dof_map;
 
diff --git a/src/ipc/utils/eigen_ext.hpp b/src/ipc/utils/eigen_ext.hpp
@@ -138,10 +138,16 @@ project_to_pd(
     const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& A,
     double eps = 1e-8);
 
-enum class PSDProjectionMethod { NONE, CLAMP, ABS };
+/// @brief Enumeration of implemented PSD projection methods
+enum class PSDProjectionMethod {
+    NONE,  ///< No PSD projection
+    CLAMP, ///< Clamp negative eigenvalues to zero
+    ABS    ///< Flip negative eigenvalues to positive
+};
 
 /// @brief Matrix projection onto positive semi-definite cone
 /// @param A Symmetric matrix to project
+/// @param method PSD projection method
 /// @return Projected matrix
 template <
     typename _Scalar,
