From dff890f60b92cab30c379b6693b4dc6c5f7387dc Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sat, 18 Mar 2023 13:59:21 -0400 Subject: [PATCH] Add tutorial (#31) --- docs/requirements.txt | 2 + docs/source/_static/css/custom.css | 44 ++ docs/source/conf.py | 28 +- docs/source/cpp-api/broad_phase.rst | 20 +- docs/source/cpp-api/candidates.rst | 34 ++ docs/source/index.md | 11 + docs/source/python-api/broad_phase.rst | 43 +- docs/source/python-api/candidates.rst | 59 ++- .../python-api/collision_constraints.rst | 27 +- docs/source/python-api/collision_mesh.rst | 3 +- docs/source/python-api/friction.rst | 24 +- docs/source/tutorial/getting_started.rst | 462 ++++++++++++++++++ docs/source/tutorial/misc.rst | 64 +++ docs/source/tutorial/simulation.rst | 89 ++++ python/src/collision_mesh.cpp | 40 +- src/ipc/collision_mesh.hpp | 54 +- src/ipc/collisions/collision_constraints.hpp | 3 +- 17 files changed, 906 insertions(+), 101 deletions(-) create mode 100644 docs/source/tutorial/getting_started.rst create mode 100644 docs/source/tutorial/misc.rst create mode 100644 docs/source/tutorial/simulation.rst diff --git a/docs/requirements.txt b/docs/requirements.txt index 08f2d54b..d352d5cb 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,6 +1,8 @@ sphinx sphinx-sitemap sphinx-immaterial +autoclasstoc +# sphinx-autodoc-toolbox breathe nbsphinx pandoc diff --git a/docs/source/_static/css/custom.css b/docs/source/_static/css/custom.css index 915d2e4b..0e6607e8 100644 --- a/docs/source/_static/css/custom.css +++ b/docs/source/_static/css/custom.css @@ -2,6 +2,50 @@ span.colon { margin-left: -1em; } +dl.cpp.objdesc, dl.py.objdesc { + /* border: 0.05rem solid var(--md-primary-fg-color); */ + border: 0.05rem solid rgb(68, 138, 255); + box-shadow: var(--md-shadow-z1); + border-radius: 0.4rem; + overflow: hidden; +} + +dl.cpp.objdesc > dt, dl.py.objdesc > dt { + padding-left: 0.5rem; + padding-right: 0.5rem; +} + +dl.cpp.objdesc > dd, dl.py.objdesc > dd { + margin: 0 1.875em; +} + +dl.py.objdesc > dd > details.toggle-details { + border: 0.05rem solid var(--md-primary-fg-color); + box-shadow: var(--md-shadow-z1); + border-radius: 0.4rem; + overflow: hidden; +} + +dl.py.objdesc > dd > details.toggle-details > summary { + border: 0; + font-size: 0.8rem; +} + +dl.py.objdesc > dd > details.toggle-details > summary::before, +dl.py.objdesc > dd > details.toggle-details > summary::after { + background-color: var(--md-primary-fg-color--light); + /* background-color: currentcolor; */ + top: auto; +} + +dl.py.objdesc > dd > details.toggle-details > summary > svg.tb-icon { + display: none; +} + +table.autosummary { + font-size: 0.75rem !important; +} + /* .breatheparameterlist li tt + p { display: inline; } diff --git a/docs/source/conf.py b/docs/source/conf.py index bf66c63d..9dd598bb 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -37,7 +37,9 @@ # extensions coming with Sphinx (named "sphinx.ext.*") or your custom # ones. extensions = [ + "autoclasstoc", "sphinx.ext.autodoc", + "sphinx.ext.autosummary", "sphinx.ext.napoleon", "sphinx.ext.intersphinx", "sphinx.ext.autosectionlabel", @@ -51,7 +53,13 @@ "myst_parser", "nbsphinx", "sphinx_immaterial", - "sphinx_immaterial.apidoc.python.apigen" + "sphinx_immaterial.apidoc.python.apigen", + "sphinx_immaterial.apidoc.format_signatures", + # 'sphinx_autodoc_toolbox.collapse', +] + +object_description_options = [ + ("cpp:.*", dict(clang_format_style={"BasedOnStyle": "WebKit"})), ] source_suffix = { @@ -69,10 +77,23 @@ project: "../build/doxyoutput/xml" } breathe_default_project = project -breathe_default_members = ("members", "undoc-members") +breathe_default_members = ( + "members", + "undoc-members", + "protected-members", + "private-members", +) breathe_show_define_initializer = True # breathe_show_include = True +autodoc_default_options = { + "members": True, + "undoc-members": True, + "private-members": True, + 'special-members': True, + 'show-inheritance': True, +} + # -- GraphViz configuration ---------------------------------- graphviz_output_format = 'svg' @@ -163,7 +184,8 @@ "navigation.tracking", "search.highlight", "search.share", - "toc.follow" + "toc.follow", + "content.tabs.link" ], "font": { diff --git a/docs/source/cpp-api/broad_phase.rst b/docs/source/cpp-api/broad_phase.rst index 6130e4c9..75fc17ad 100644 --- a/docs/source/cpp-api/broad_phase.rst +++ b/docs/source/cpp-api/broad_phase.rst @@ -5,16 +5,34 @@ Broad Phase .. doxygenvariable:: ipc::DEFAULT_BROAD_PHASE_METHOD +Broad Phase +----------- + .. doxygenclass:: ipc::BroadPhase -.. doxygenclass:: ipc::AABB +Brute Force +----------- .. doxygenclass:: ipc::BruteForce +Hash Grid +--------- + .. doxygenclass:: ipc::HashGrid +Spatial Hash +------------ + .. doxygenclass:: ipc::SpatialHash +Sweep and Tiniest Queue +----------------------- + .. doxygenclass:: ipc::SweepAndTiniestQueue .. .. doxygenclass:: ipc::SweepAndTiniestQueueGPU + +AABB +---- + +.. doxygenclass:: ipc::AABB \ No newline at end of file diff --git a/docs/source/cpp-api/candidates.rst b/docs/source/cpp-api/candidates.rst index 2a4a4045..f9f85e3f 100644 --- a/docs/source/cpp-api/candidates.rst +++ b/docs/source/cpp-api/candidates.rst @@ -1,13 +1,47 @@ Candidates ========== +Candidates +---------- + .. doxygenclass:: ipc::Candidates +Collision Stencil +----------------- + .. doxygenclass:: ipc::CollisionStencil + + +Continuous Collision Candidate +------------------------------ + .. doxygenclass:: ipc::ContinuousCollisionCandidate + +Vertex-Vertex Candidate +----------------------- + .. doxygenclass:: ipc::VertexVertexCandidate + + +Edge-Vertex Candidate +--------------------- + .. doxygenclass:: ipc::EdgeVertexCandidate + + +Edge-Edge Candidate +------------------- + .. doxygenclass:: ipc::EdgeEdgeCandidate + + +Edge-Face Candidate +------------------- + .. doxygenclass:: ipc::EdgeFaceCandidate + +Face-Vertex Candidate +--------------------- + .. doxygenclass:: ipc::FaceVertexCandidate \ No newline at end of file diff --git a/docs/source/index.md b/docs/source/index.md index 12d619e4..e26c20f4 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -7,6 +7,15 @@ changelog license ``` +```{toctree} +:caption: Tutorial +:hidden: + +tutorial/getting_started.rst +tutorial/simulation.rst +tutorial/misc.rst +``` + ```{toctree} :caption: C++ :hidden: @@ -48,5 +57,7 @@ style_guide Code of Conduct ``` + + ```{include} ../../README.md ``` diff --git a/docs/source/python-api/broad_phase.rst b/docs/source/python-api/broad_phase.rst index 9808dbc9..6ca877e4 100644 --- a/docs/source/python-api/broad_phase.rst +++ b/docs/source/python-api/broad_phase.rst @@ -2,27 +2,52 @@ Broad Phase =========== .. autoclass:: ipctk.BroadPhaseMethod - :members: + + .. autoclasstoc:: .. .. autovariable:: ipctk.DEFAULT_BROAD_PHASE_METHOD +Broad Phase +----------- + .. autoclass:: ipctk.BroadPhase - :members: -.. autoclass:: ipctk.AABB - :members: + .. autoclasstoc:: + +Brute Force +----------- .. autoclass:: ipctk.BruteForce - :members: + + .. autoclasstoc:: + +Hash Grid +--------- .. autoclass:: ipctk.HashGrid - :members: + + .. autoclasstoc:: + +Spatial Hash +------------ .. autoclass:: ipctk.SpatialHash - :members: + + .. autoclasstoc:: + +Sweep and Tiniest Queue +----------------------- .. autoclass:: ipctk.SweepAndTiniestQueue - :members: + + .. autoclasstoc:: .. .. autoclass:: ipctk.SweepAndTiniestQueueGPU -.. :members: \ No newline at end of file +.. :members: + +AABB +---- + +.. autoclass:: ipctk.AABB + + .. autoclasstoc:: \ No newline at end of file diff --git a/docs/source/python-api/candidates.rst b/docs/source/python-api/candidates.rst index f14f370c..d575f540 100644 --- a/docs/source/python-api/candidates.rst +++ b/docs/source/python-api/candidates.rst @@ -1,26 +1,59 @@ Candidates ========== +Candidates +---------- + .. autoclass:: ipctk.Candidates - :members: + + .. autoclasstoc:: + +Collision Stencil +----------------- + +.. error:: + This ``autoclass`` is not working. I don't know why. I have to investigate. .. .. autoclass:: ipctk.CollisionStencil -.. :members: + +Continuous Collision Candidate +------------------------------ + .. autoclass:: ipctk.ContinuousCollisionCandidate - :members: + + .. autoclasstoc:: + +Vertex-Vertex Candidate +----------------------- .. autoclass:: ipctk.VertexVertexCandidate - :members: - :show-inheritance: + + .. autoclasstoc:: + +Edge-Vertex Candidate +--------------------- + .. autoclass:: ipctk.EdgeVertexCandidate - :members: - :show-inheritance: + + .. autoclasstoc:: + +Edge-Edge Candidate +------------------- + .. autoclass:: ipctk.EdgeEdgeCandidate - :members: - :show-inheritance: + + .. autoclasstoc:: + +Edge-Face Candidate +------------------- + .. autoclass:: ipctk.EdgeFaceCandidate - :members: - :show-inheritance: + + .. autoclasstoc:: + +Face-Vertex Candidate +--------------------- + .. autoclass:: ipctk.FaceVertexCandidate - :members: - :show-inheritance: \ No newline at end of file + + .. autoclasstoc:: \ No newline at end of file diff --git a/docs/source/python-api/collision_constraints.rst b/docs/source/python-api/collision_constraints.rst index b3ce5166..77d0faf2 100644 --- a/docs/source/python-api/collision_constraints.rst +++ b/docs/source/python-api/collision_constraints.rst @@ -5,46 +5,47 @@ Collision Constraints --------------------- .. autoclass:: ipctk.CollisionConstraints - :members: + + .. autoclasstoc:: Collision Constraint -------------------- .. autoclass:: ipctk.CollisionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Vertex-Vertex Collision Constraint ---------------------------------- .. autoclass:: ipctk.VertexVertexConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Edge-Vertex Collision Constraint -------------------------------- .. autoclass:: ipctk.EdgeVertexConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Edge-Edge Collision Constraint ------------------------------ .. autoclass:: ipctk.EdgeEdgeConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Face-Vertex Collision Constraint -------------------------------- .. autoclass:: ipctk.FaceVertexConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Plane-Vertex Collision Constraint --------------------------------- .. autoclass:: ipctk.PlaneVertexConstraint - :members: - :show-inheritance: \ No newline at end of file + + .. autoclasstoc:: \ No newline at end of file diff --git a/docs/source/python-api/collision_mesh.rst b/docs/source/python-api/collision_mesh.rst index c53cbc44..8e51f5f1 100644 --- a/docs/source/python-api/collision_mesh.rst +++ b/docs/source/python-api/collision_mesh.rst @@ -2,4 +2,5 @@ Collision Mesh ============== .. autoclass:: ipctk.CollisionMesh - :members: \ No newline at end of file + + .. autoclasstoc:: \ No newline at end of file diff --git a/docs/source/python-api/friction.rst b/docs/source/python-api/friction.rst index 4a8ef273..ad2cefef 100644 --- a/docs/source/python-api/friction.rst +++ b/docs/source/python-api/friction.rst @@ -5,43 +5,43 @@ Friction Constraints -------------------- .. autoclass:: ipctk.FrictionConstraints - :members: - :show-inheritance: + + .. autoclasstoc:: Friction Constraint ^^^^^^^^^^^^^^^^^^^ .. autoclass:: ipctk.FrictionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Vertex-Vertex Friction Constraint ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ipctk.VertexVertexFrictionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Edge-Vertex Friction Constraint ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ipctk.EdgeVertexFrictionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Edge-Edge Friction Constraint ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ipctk.EdgeEdgeFrictionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Triangle-Vertex Friction Constraint ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: ipctk.FaceVertexFrictionConstraint - :members: - :show-inheritance: + + .. autoclasstoc:: Smooth Mollifier ---------------- diff --git a/docs/source/tutorial/getting_started.rst b/docs/source/tutorial/getting_started.rst new file mode 100644 index 00000000..cd2b65ee --- /dev/null +++ b/docs/source/tutorial/getting_started.rst @@ -0,0 +1,462 @@ +Getting Started +=============== + +This tutorial will walk you through the basics of using the IPC Toolkit. + +Collision Mesh +-------------- + +First, we need to create a collision mesh. The ``CollisionMesh`` data structure represents the surface geometry used for collision processing. + +We will start by creating a collision mesh from a ``bunny.obj`` mesh file (you can find the mesh `here `_): + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + #include + + #include + #include + #include + #include + + // ... + + Eigen::MatrixXd rest_positions; + Eigen::MatrixXi edges, faces; + igl::readOBJ("bunny.obj", rest_positions, faces); + igl::edges(faces, edges); + + ipc::CollisionMesh collision_mesh(rest_positions, edges, faces); + + .. md-tab-item:: Python + + .. code-block:: python + + import ipctk + + import meshio + + mesh = meshio.read("bunny.obj") + rest_positions = mesh.points + faces = mesh.cells[0].data + edges = ipctk.edges(faces) + + collision_mesh = ipctk.CollisionMesh(rest_positions, edges, faces) + +The matrix ``rest_positions`` contains the undeformed positions of the vertices. It should be :math:`|V| \times d` where :math:`|V|` is the number of vertices and :math:`d \in \{2, 3\}` is the dimension. +Each row of the ``edges`` and ``faces`` matrices contain the vertex IDs (row number in ``rest_positions``) of the edge or face's vertices. +The size of ``edges`` and ``faces`` are ``#E x 2`` and ``#F x 3`` respectively (``#E`` and ``#F`` are the number of edges and faces). + +.. NOTE:: + Only linear triangular faces are supported. If your mesh has nonlinear or non-triangular faces, you will need to triangulate them. + +.. NOTE:: + In 2D only the ``edges`` matrix is required. In 3D both ``edges`` and ``faces`` are required. + +Collision Constraints +--------------------- + +Now that we have a collision mesh, we can compute the contact potential. To do this we first need to build the set of active collision constraints (``CollisionConstraints``). + +To start we need the current positions of the ``vertices``. For this tutorial, let us use squash the bunny has to 1% of its original height. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::MatrixXd vertices = collision_mesh.rest_positions(); + vertices.col(1) *= 0.01; // Squash the bunny in the y-direction + + .. md-tab-item:: Python + + .. code-block:: python + + vertices = collision_mesh.rest_positions() + vertices[:, 1] *= 0.01 # Squash the bunny in the y-direction + +Using these deformed positions, we can build the set of active collision constraints. +For this, we need to specify :math:`\hat{d}` (``dhat``), the IPC barrier activation distance. +We will use a value of :math:`\hat{d} = 10^{-3}`. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + const double dhat = 1e-3; + + ipc::CollisionConstraints collision_constraints; + collision_constraints.build(collision_mesh, vertices, dhat); + + .. md-tab-item:: Python + + .. code-block:: python + + dhat = 1e-3 + + collision_constraints = ipctk.CollisionConstraints() + collision_constraints.build(collision_mesh, vertices, dhat) + +This will automatically use a spatial data structure to perform a broad-phase culling and then perform a narrow-phase culling by computing distances (discarding any constraint with a distance :math:`> \hat{d}`). + +Contact Potential +^^^^^^^^^^^^^^^^^ + +Now we can compute the contact potential using the ``CollisionConstraints``. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + double contact_potential = collision_constraints.compute_potential( + collision_mesh, vertices, dhat); + + .. md-tab-item:: Python + + .. code-block:: python + + contact_potential = collision_constraints.compute_potential( + collision_mesh, vertices, dhat) + +This returns a scalar value ``contact_potential`` which is the sum of the contact potential for each active constraint. + +Mathematically this is defined as + +.. math:: + B(x) = \sum_{k \in C} b\left(d_k(x)\right), + +where :math:`C` is the active collision constraints, :math:`d_k` is the distance (squared) of the :math:`k`-th active constraint, and :math:`b` is IPC's C2-clamped log-barrier function. + +.. NOTE:: + This is not premultiplied by the barrier stiffness :math:`\kappa`. + +Contact Potential Derivative +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +We can also compute the first and second derivatives of the contact potential with respect to the vertex positions. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::VectorXd contact_potential_grad = + collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat); + + Eigen::SparseMatrix contact_potential_hess = + collision_constraints.compute_potential_hessian( + collision_mesh, vertices, dhat); + + .. md-tab-item:: Python + + .. code-block:: python + + contact_potential_grad = collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat) + + contact_potential_hess = collision_constraints.compute_potential_hessian( + collision_mesh, vertices, dhat) + +These return the gradient and hessian of the contact potential as a dense vector and sparse matrix respectively. + +The derivatives are taken with respect to the row-wise flattened vertices. That is, for ``vertices`` + +.. math:: + \begin{bmatrix} + x_1 & y_1 & z_1 \\ + & \vdots & \\ + x_n & y_n & z_n \\ + \end{bmatrix}, + +you will get the gradient of size :math:`|V|d \times 1` with the order + +.. math:: + \begin{bmatrix} + \frac{\partial B}{\partial x_1} & + \frac{\partial B}{\partial y_1} & + \frac{\partial B}{\partial z_1} & + \cdots & + \frac{\partial B}{\partial x_n} & + \frac{\partial B}{\partial y_n} & + \frac{\partial B}{\partial z_n} + \end{bmatrix}^T, + +and the Hessian of size :math:`|V|d \times |V|d` with the order + +.. math:: + \begin{bmatrix} + \frac{\partial^2 B}{\partial x_1^2} & + \frac{\partial^2 B}{\partial x_1 \partial y_1} & + \frac{\partial^2 B}{\partial x_1 \partial z_1} & + \cdots & + \frac{\partial^2 B}{\partial x_1 \partial x_n} & + \frac{\partial^2 B}{\partial x_1 \partial y_n} & + \frac{\partial^2 B}{\partial x_1 \partial z_n} \\ + % + \frac{\partial^2 B}{\partial y_1 \partial x_1} & + \frac{\partial^2 B}{\partial y_1^2} & + \frac{\partial^2 B}{\partial y_1 \partial z_1} & + \cdots & + \frac{\partial^2 B}{\partial y_1 \partial x_n} & + \frac{\partial^2 B}{\partial y_1 \partial y_n} & + \frac{\partial^2 B}{\partial y_1 \partial z_n} \\ + % + \frac{\partial^2 B}{\partial z_1 \partial x_1} & + \frac{\partial^2 B}{\partial z_1 \partial y_1} & + \frac{\partial^2 B}{\partial z_1^2} & + \cdots & + \frac{\partial^2 B}{\partial z_1 \partial x_n} & + \frac{\partial^2 B}{\partial z_1 \partial y_n} & + \frac{\partial^2 B}{\partial z_1 \partial z_n} \\ + % + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + % + \frac{\partial^2 B}{\partial x_n \partial x_1} & + \frac{\partial^2 B}{\partial x_n \partial y_1} & + \frac{\partial^2 B}{\partial x_n \partial z_1} & + \cdots & + \frac{\partial^2 B}{\partial x_n^2} & + \frac{\partial^2 B}{\partial x_n \partial y_n} & + \frac{\partial^2 B}{\partial x_n \partial z_n} \\ + % + \frac{\partial^2 B}{\partial y_n \partial x_1} & + \frac{\partial^2 B}{\partial y_n \partial y_1} & + \frac{\partial^2 B}{\partial y_n \partial z_1} & + \cdots & + \frac{\partial^2 B}{\partial y_n \partial x_n} & + \frac{\partial^2 B}{\partial y_n^2} & + \frac{\partial^2 B}{\partial y_n \partial z_n} \\ + % + \frac{\partial^2 B}{\partial z_n \partial x_1} & + \frac{\partial^2 B}{\partial z_n \partial y_1} & + \frac{\partial^2 B}{\partial z_n \partial z_1} & + \cdots + & + \frac{\partial^2 B}{\partial z_n \partial x_n} & + \frac{\partial^2 B}{\partial z_n \partial y_n} & + \frac{\partial^2 B}{\partial z_n^2} + \end{bmatrix}. + +Adaptive Barrier Stiffness +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The last piece of the contact potential is the barrier stiffness. This is a weight that is multiplied by the barrier potential to better scale it relative to the energy potential. This can be a fixed value or adaptive. + +To compute the adaptive barrier stiffness, we can use two functions: ``initial_barrier_stiffness`` and ``update_barrier_stiffness``. The function ``initial_barrier_stiffness`` compute the initial value from the current energy and contact potential gradients. This function also provides a minimum and maximum value for the barrier stiffness. The function ``update_barrier_stiffness`` updates the barrier stiffness if the minimum distance has become too small. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + // (beginning of nonlinear solve) + + Eigen::VectorXd grad_energy = ...; // gradient of elastic energy potential + Eigen::VectorXd grad_contact = collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat); + + double bbox_diagonal = ipc::world_bbox_diagonal_length(vertices); + + double max_barrier_stiffness; // output of initial_barrier_stiffness + double barrier_stiffness = ipc::initial_barrier_stiffness( + bbox_diagonal, dhat, avg_mass, grad_energy, grad_contact, + max_barrier_stiffness); + + double prev_distance = + collision_constraints.compute_minimum_distance( + collision_mesh, vertices); + + // ... + + // (end of nonlinear iteration) + + double curr_distance = + collision_constraints.compute_minimum_distance(collision_mesh, vertices); + + barrier_stiffness = ipc::update_barrier_stiffness( + prev_distance, curr_distance, max_barrier_stiffness, barrier_stiffness, + bbox_diagonal); + + prev_distance = curr_distance; + + // (next iteration) + + .. md-tab-item:: Python + + .. code-block:: python + + # (beginning of nonlinear solve) + + grad_energy = ... # gradient of elastic energy potential + grad_contact = collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat) + + bbox_diagonal = ipctk.world_bbox_diagonal_length(vertices) + + barrier_stiffness, max_barrier_stiffness = ipctk.initial_barrier_stiffness( + bbox_diagonal, dhat, avg_mass, grad_energy, grad_contact, + max_barrier_stiffness) + + prev_distance = collision_constraints.compute_minimum_distance( + collision_mesh, vertices) + + # ... + + # (end of nonlinear iteration) + + curr_distance = collision_constraints.compute_minimum_distance( + collision_mesh, vertices) + + barrier_stiffness = ipctk.update_barrier_stiffness( + prev_distance, curr_distance, max_barrier_stiffness, barrier_stiffness, + bbox_diagonal) + + prev_distance = curr_distance + + # (next iteration) + +Friction +-------- + +Computing the friction dissipative potential is similar to the contact potential, but because it is a lagged model, we need to build it from a fixed set of constraints. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + ipc::FrictionConstraints friction_constraints; + collision_constraints.build( + collision_mesh, vertices, contact_constraints, dhat, barrier_stiffness, mu); + + .. md-tab-item:: Python + + .. code-block:: python + + friction_constraints = ipctk.FrictionConstraints() + friction_constraints.build( + collision_mesh, vertices, contact_constraints, dhat, barrier_stiffness, mu) + +Here ``mu`` (:math:`\mu`) is the (global) coefficient of friction, and ``barrier_stiffness`` (:math:`\kappa`) is the barrier stiffness. + +Friction Dissipative Potential +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Now we can compute the friction dissipative potential using the ``FrictionConstraints``. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + double friction_potential = friction_constraints.compute_potential( + collision_mesh, velocity, epsv_times_h); + + .. md-tab-item:: Python + + .. code-block:: python + + friction_potential = friction_constraints.compute_potential( + collision_mesh, velocity, epsv_times_h) + +Here ``epsv_times_h`` (:math:`\epsilon_v h`) is the static friction threshold (in units of velocity) used to smoothly transition from dynamic to static friction. + +Notice the friction potential is a function of the velocities rather than the positions. We can compute the velocities directly from the current and previous position(s) based on our time-integration scheme. For example, if we are using backward Euler integration, then the velocity is + +.. math:: + v = \frac{x - x_{n-1}}{h}, + +where :math:`x` is the current position, :math:`x_{n-1}` is the previous position, and :math:`h` is the time step size. + +This returns a scalar value ``friction_potential`` which is the sum of the individual friction potentials. + +Mathematically this is defined as + +.. math:: + D(x) = \sum_{k \in C} \mu\lambda_k^nf_0\left(\|T_k^Tv\|; \epsilon_v h\right), + +where :math:`C` is the lagged collision constraints, :math:`\lambda_k^n` is the normal force magnitude for the :math:`k`-th contact, :math:`T_k` is the tangential basis for the :math:`k`-th contact, and :math:`f_0` is the smooth friction function used to approximate the non-smooth transition from dynamic to static friction. + +Derivatives +^^^^^^^^^^^ + +We can also compute the first and second derivatives of the friction dissipative potential with respect to the velocities. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::VectorXd friction_potential_grad = + friction_constraints.compute_potential_gradient( + collision_mesh, velocity, epsv_times_h); + + Eigen::SparseMatrix friction_potential_hess = + friction_constraints.compute_potential_hessian( + collision_mesh, velocity, epsv_times_h); + + .. md-tab-item:: Python + + .. code-block:: python + + friction_potential_grad = friction_constraints.compute_potential_gradient( + collision_mesh, velocity, epsv_times_h) + + friction_potential_hess = friction_constraints.compute_potential_hessian( + collision_mesh, velocity, epsv_times_h) + +Continuous Collision Detection +------------------------------ + +The last high-level component of the IPC Toolkit library is continuous collision detection (CCD). This is a method for determining if and at what time two objects will collide. This can be incorporated in a simulation nonlinear solver's line-search to determine the maximum step size allowable before a collision occurs. + +There are two main functions for doing this: ``is_step_collision_free`` and ``compute_collision_free_stepsize``. The former returns a boolean value indicating if the step is collision-free, and the latter returns the maximum step size that is collision-free. Both functions take the same arguments, but ``compute_collision_free_stepsize`` is the more convenient function to use as it returns the maximum step size. + +The following example determines the maximum step size allowable between the rest_positions and the squashed bunny. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::MatrixXd vertices_t0 = collision_mesh.rest_positions(); + Eigen::MatrixXd vertices_t1 = vertices_t0; + vertices_t1.col(1) *= 0.01; + + double max_step_size = compute_collision_free_stepsize( + collision_mesh, vertices_t0, vertices_t1); + + Eigen::MatrixXd collision_free_vertices = + (vertices_t1 - vertices_t0) * max_step_size + vertices_t0; + assert(is_step_collision_free(mesh, vertices_t0, collision_free_vertices)); + + .. md-tab-item:: Python + + .. code-block:: python + + vertices_t0 = collision_mesh.rest_positions() + vertices_t1 = vertices_t0 + vertices_t1[:, 1] *= 0.01 + + max_step_size = compute_collision_free_stepsize( + collision_mesh, vertices_t0, vertices_t1) + + collision_free_vertices = + (vertices_t1 - vertices_t0) * max_step_size + vertices_t0 + assert(is_step_collision_free(mesh, vertices_t0, collision_free_vertices)) \ No newline at end of file diff --git a/docs/source/tutorial/misc.rst b/docs/source/tutorial/misc.rst new file mode 100644 index 00000000..dd820ab8 --- /dev/null +++ b/docs/source/tutorial/misc.rst @@ -0,0 +1,64 @@ +Miscellaneous +============= + +Static Intersection Checks +-------------------------- + +Static intersection checks are not a core part of the IPC algorithm, but they are useful for debugging and verifying your input is intersection free. The IPC Toolkit provides a function for checking if a mesh is intersection free: + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + bool is_intersecting = ipc::has_intersections(collision_mesh, displaced); + + .. md-tab-item:: Python + + .. code-block:: python + + is_intersecting = ipctk.has_intersections(collision_mesh, displaced) + + + +Logger +------ + +The IPC Toolkit uses the `spdlog `_ library for logging. By default, the IPC Toolkit will log to ``stdout``. To change this behavior, you can set the logger using the ``ipc::set_logger`` function. For example, to additionally log to a file, you can do the following: + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + std::string log_file = "log.txt"; + + std::vector sinks; + sinks.emplace_back(std::make_shared()); + sinks.emplace_back(std::make_shared(log_file, /*truncate=*/true)); + + ipc::set_logger(std::make_shared("ipctk", sinks.begin(), sinks.end())); + + .. md-tab-item:: Python + + .. code-block:: python + + # Unfourtnately, this is not yet supported in Python. + +You can also set the log level: + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + ipc::logger().set_level(spdlog::level::debug); + + .. md-tab-item:: Python + + .. code-block:: python + + ipctk.set_logger_level(ipctk.LoggerLevel.debug) diff --git a/docs/source/tutorial/simulation.rst b/docs/source/tutorial/simulation.rst new file mode 100644 index 00000000..e9732a09 --- /dev/null +++ b/docs/source/tutorial/simulation.rst @@ -0,0 +1,89 @@ +Physical Simulation +=================== + +While the IPC Toolkit provides all the principle components of the IPC algorithm, it does not provide a complete simulation framework. Instead, it provides building blocks that can be used to integrate the IPC algorithm into a physical simulation. If all you want is a complete simulation framework using the IPC algorithm, then you should check out our other project `PolyFEM `_ which uses the IPC Toolkit for its collision handling. + +We provide several helper functions to make your job easier. The following examples show how to use these functions. + +Volumetric Meshes +----------------- + +The IPC Toolkit only handles surface meshes (through the ``CollisionMesh``). However, the finite element method often relies on volumetric discretizations. In this case the computed gradients and Hessians need to be mapped back to the full volumetric mesh. The ``CollisionMesh`` class provides this functionality. + +From the full (volumetric) mesh vertices and surface edges/faces which index into the full mesh vertices, you can build a ``CollisionMesh`` using the function ``CollisionMesh::build_from_full_mesh``. This will internally build and store a selection matrix that goes from the full to surface vertices as well as map the edge/faces entries accordingly. + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::MatrixXd full_rest_positions; + Eigen::MatrixXi tets; + // TODO: Show how to load a volumetric mesh from a file (e.g., using MshIO) + + // Faces of the surface mesh with indices into full_rest_positions + Eigen::MatrixXd faces; + igl::boundary_facets(tets, faces); + + // Edges of the surface mesh with indices into full_rest_positions + Eigen::MatrixXi edges; + igl::edges(faces, edges); + + std::vector 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); + + .. md-tab-item:: Python + + .. code-block:: python + + mesh = meshio.read("bunny.msh") + full_rest_positions = mesh.points + tets = mesh.cells[0].data + + faces = igl.boundary_facets(tets) # pip install libigl + edges = ipctk.edges(faces) # same as igl.edges + + collision_mesh = ipctk.CollisionMesh.build_from_full_mesh( + 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``): + +.. md-tab-set:: + + .. md-tab-item:: C++ + + .. code-block:: c++ + + Eigen::VectorXd grad = collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat); + Eigen::VectorXd grad_full = collision_mesh.to_full_dof(grad); + + Eigen::SparseMatrix hess = collision_constraints.compute_potential_hessian( + collision_mesh, vertices, dhat); + Eigen::SparseMatrix hess_full = collision_mesh.to_full_dof(hess); + + .. md-tab-item:: Python + + .. code-block:: python + + grad = collision_constraints.compute_potential_gradient( + collision_mesh, vertices, dhat); + grad_full = collision_mesh.to_full_dof(grad); + + hess = collision_constraints.compute_potential_hessian( + collision_mesh, vertices, dhat); + hess_full = collision_mesh.to_full_dof(hess); + +Positive Semi-Definite Projection +--------------------------------- + +As described in [IPC]_ the Hessian of the potentials can be indefinite. This is problematic when using the Hessian in a Newton step [IPC]_. 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. + +------------ + +.. rubric:: References + +.. [IPC] Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, Danny M. Kaufman. 2020. Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics. *ACM Transactions on Graphics (SIGGRAPH).* \ No newline at end of file diff --git a/python/src/collision_mesh.cpp b/python/src/collision_mesh.cpp index 40097895..38f42d5a 100644 --- a/python/src/collision_mesh.cpp +++ b/python/src/collision_mesh.cpp @@ -16,9 +16,9 @@ void define_collision_mesh(py::module_& m) Construct a new Collision Mesh object directly from the collision mesh vertices. Parameters: - rest_positions: The vertices of the collision mesh at rest. - edges: The edges of the collision mesh. - faces: The faces of the collision mesh. + rest_positions: The vertices of the collision mesh at rest (#V × dim). + edges: The edges of the collision mesh (#E × 2). + faces: The faces of the collision mesh (#F × 3). displacement_map: The displacement mapping from displacements on the full mesh to the collision mesh. )ipc_Qu8mg5v7", py::arg("rest_positions"), py::arg("edges"), py::arg("faces"), @@ -33,9 +33,9 @@ void define_collision_mesh(py::module_& m) Parameters: include_vertex: Vector of bools indicating whether each vertex should be included in the collision mesh. - full_rest_positions: The vertices of the full mesh at rest. - edges: The edges of the collision mesh indexed into the full mesh vertices. - faces: The faces of the collision mesh indexed into the full mesh vertices. + full_rest_positions: The vertices of the full mesh at rest (#V × dim). + edges: The edges of the collision mesh indexed into the full mesh vertices (#E × 2). + faces: The faces of the collision mesh indexed into the full mesh vertices (#F × 3). displacement_map: The displacement mapping from displacements on the full mesh to the collision mesh. )ipc_Qu8mg5v7", py::arg("include_vertex"), py::arg("full_rest_positions"), @@ -47,9 +47,9 @@ void define_collision_mesh(py::module_& m) Helper function that automatically builds include_vertex using construct_is_on_surface. Parameters: - full_rest_positions: The full vertices at rest. - edges: The edge matrix of mesh. - faces: The face matrix of mesh. + full_rest_positions: The full vertices at rest (#V × dim). + edges: The edge matrix of mesh (#E × 2). + faces: The face matrix of mesh (#F × 3). Returns: Constructed CollisionMesh. @@ -82,35 +82,35 @@ void define_collision_mesh(py::module_& m) "Get the number of degrees of freedom in the full mesh.") .def_property_readonly( "rest_positions", &CollisionMesh::rest_positions, - "Get the vertices of the collision mesh at rest.") + "Get the vertices of the collision mesh at rest (#V × dim).") .def_property_readonly( "edges", &CollisionMesh::edges, - "Get the edges of the collision mesh.") + "Get the edges of the collision mesh (#E × 2).") .def_property_readonly( "faces", &CollisionMesh::faces, - "Get the faces of the collision mesh.") + "Get the faces of the collision mesh (#F × 3).") .def_property_readonly( "faces_to_edges", &CollisionMesh::faces_to_edges, - "Get the mapping from faces to edges of the collision mesh.") + "Get the mapping from faces to edges of the collision mesh (#F × 3).") .def( "vertices", &CollisionMesh::vertices, R"ipc_Qu8mg5v7( Compute the vertex positions from the positions of the full mesh. Parameters: - full_positions: The vertex positions of the full mesh. + full_positions: The vertex positions of the full mesh (#FV × dim). Returns: - The vertex positions of the collision mesh. + The vertex positions of the collision mesh (#V × dim). )ipc_Qu8mg5v7", py::arg("full_positions")) .def( "displace_vertices", &CollisionMesh::displace_vertices, R"ipc_Qu8mg5v7( - Compute the vertex positions from vertex displacements on the full mesh. + Compute the vertex positions from vertex displacements on the full mesh (#FV × dim). Parameters: - full_displacements: The vertex displacements on the full mesh. + full_displacements: The vertex displacements on the full mesh (#V × dim). Returns: The vertex positions of the collision mesh. @@ -268,7 +268,7 @@ void define_collision_mesh(py::module_& m) Parameters: num_vertices: The number of vertices in the mesh. - edges: The surface edges of the mesh. + edges: The surface edges of the mesh (#E × 2). Returns: A vector of bools indicating whether each vertex is on the surface. @@ -281,8 +281,8 @@ void define_collision_mesh(py::module_& m) Construct a matrix that maps from the faces' edges to rows in the edges matrix. Parameters: - faces: The face matrix of mesh. - edges: The edge matrix of mesh. + faces: The face matrix of mesh (#F × 3). + edges: The edge matrix of mesh (#E × 2). Returns: Matrix that maps from the faces' edges to rows in the edges matrix. diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp index 2acc9201..738bf4fa 100644 --- a/src/ipc/collision_mesh.hpp +++ b/src/ipc/collision_mesh.hpp @@ -15,9 +15,9 @@ class CollisionMesh { CollisionMesh() { } /// @brief Construct a new Collision Mesh object directly from the collision mesh vertices. - /// @param rest_positions The vertices of the collision mesh at rest. - /// @param edges The edges of the collision mesh. - /// @param faces The faces of the collision mesh. + /// @param rest_positions The vertices of the collision mesh at rest (#V × dim). + /// @param edges The edges of the collision mesh (#E × 2). + /// @param faces The faces of the collision mesh (#F × 3). /// @param displacement_map The displacement mapping from displacements on the full mesh to the collision mesh. CollisionMesh( const Eigen::MatrixXd& rest_positions, @@ -28,9 +28,9 @@ class CollisionMesh { /// @brief Construct a new Collision Mesh object from a full mesh vertices. /// @param include_vertex Vector of bools indicating whether each vertex should be included in the collision mesh. - /// @param full_rest_positions The vertices of the full mesh at rest. - /// @param edges The edges of the collision mesh indexed into the full mesh vertices. - /// @param faces The faces of the collision mesh indexed into the full mesh vertices. + /// @param full_rest_positions The vertices of the full mesh at rest (#V × dim). + /// @param edges The edges of the collision mesh indexed into the full mesh vertices (#E × 2). + /// @param faces The faces of the collision mesh indexed into the full mesh vertices (#F × 3). /// @param displacement_map The displacement mapping from displacements on the full mesh to the collision mesh. CollisionMesh( const std::vector& include_vertex, @@ -41,9 +41,9 @@ class CollisionMesh { Eigen::SparseMatrix()); /// @brief Helper function that automatically builds include_vertex using construct_is_on_surface. - /// @param full_rest_positions The full vertices at rest. - /// @param edges The edge matrix of mesh. - /// @param faces The face matrix of mesh. + /// @param full_rest_positions The full vertices at rest (#FV × dim). + /// @param edges The edge matrix of mesh (#E × 2). + /// @param faces The face matrix of mesh (#F × 3). /// @return Constructed CollisionMesh. static CollisionMesh build_from_full_mesh( const Eigen::MatrixXd& full_rest_positions, @@ -87,16 +87,16 @@ class CollisionMesh { /// @brief Get the number of degrees of freedom in the full mesh. size_t full_ndof() const { return full_num_vertices() * dim(); } - /// @brief Get the vertices of the collision mesh at rest. + /// @brief Get the vertices of the collision mesh at rest (#V × dim). const Eigen::MatrixXd& rest_positions() const { return m_rest_positions; } - /// @brief Get the edges of the collision mesh. + /// @brief Get the edges of the collision mesh (#E × 2). const Eigen::MatrixXi& edges() const { return m_edges; } - /// @brief Get the faces of the collision mesh. + /// @brief Get the faces of the collision mesh (#F × 3). const Eigen::MatrixXi& faces() const { return m_faces; } - /// @brief Get the mapping from faces to edges of the collision mesh. + /// @brief Get the mapping from faces to edges of the collision mesh (#F × 3). const Eigen::MatrixXi& faces_to_edges() const { return m_faces_to_edges; } // const std::vector>& vertices_to_edges() const @@ -112,19 +112,19 @@ class CollisionMesh { // ----------------------------------------------------------------------- /// @brief Compute the vertex positions from the positions of the full mesh. - /// @param full_positions The vertex positions of the full mesh. - /// @return The vertex positions of the collision mesh. + /// @param full_positions The vertex positions of the full mesh (#FV × dim). + /// @return The vertex positions of the collision mesh (#V × dim). Eigen::MatrixXd vertices(const Eigen::MatrixXd& full_positions) const; /// @brief Compute the vertex positions from vertex displacements on the full mesh. - /// @param full_displacements The vertex displacements on the full mesh. - /// @return The vertex positions of the collision mesh. + /// @param full_displacements The vertex displacements on the full mesh (#FV × dim). + /// @return The vertex positions of the collision mesh (#V × dim). Eigen::MatrixXd displace_vertices(const Eigen::MatrixXd& full_displacements) const; /// @brief Map vertex displacements on the full mesh to vertex displacements on the collision mesh. - /// @param full_displacements The vertex displacements on the full mesh. - /// @return The vertex displacements on the collision mesh. + /// @param full_displacements The vertex displacements on the full mesh (#FV × dim). + /// @return The vertex displacements on the collision mesh (#V × dim). Eigen::MatrixXd map_displacements(const Eigen::MatrixXd& full_displacements) const; @@ -241,14 +241,14 @@ class CollisionMesh { /// @brief Construct a vector of bools indicating whether each vertex is on the surface. /// @param num_vertices The number of vertices in the mesh. - /// @param edges The surface edges of the mesh. + /// @param edges The surface edges of the mesh (#E × 2). /// @return A vector of bools indicating whether each vertex is on the surface. static std::vector construct_is_on_surface( const int num_vertices, const Eigen::MatrixXi& edges); /// @brief Construct a matrix that maps from the faces' edges to rows in the edges matrix. - /// @param faces The face matrix of mesh. - /// @param edges The edge matrix of mesh. + /// @param faces The face matrix of mesh (#F × 3). + /// @param edges The edge matrix of mesh (#E × 2). /// @return Matrix that maps from the faces' edges to rows in the edges matrix. static Eigen::MatrixXi construct_faces_to_edges( const Eigen::MatrixXi& faces, const Eigen::MatrixXi& edges); @@ -274,15 +274,15 @@ class CollisionMesh { // ----------------------------------------------------------------------- - /// @brief The full vertex positions at rest. + /// @brief The full vertex positions at rest (#FV × dim). Eigen::MatrixXd m_full_rest_positions; - /// @brief The vertex positions at rest. + /// @brief The vertex positions at rest (#V × dim). Eigen::MatrixXd m_rest_positions; - /// @brief Edges as rows of indicies into vertices. + /// @brief Edges as rows of indicies into vertices (#E × 2). Eigen::MatrixXi m_edges; - /// @brief Triangular faces as rows of indicies into vertices. + /// @brief Triangular faces as rows of indicies into vertices (#F × 3). Eigen::MatrixXi m_faces; - /// @brief Map from faces edges to rows of edges. + /// @brief Map from faces edges to rows of edges (#F × 3). Eigen::MatrixXi m_faces_to_edges; /// @brief Map from full vertices to collision vertices. diff --git a/src/ipc/collisions/collision_constraints.hpp b/src/ipc/collisions/collision_constraints.hpp index 2fcb5a6b..f815c576 100644 --- a/src/ipc/collisions/collision_constraints.hpp +++ b/src/ipc/collisions/collision_constraints.hpp @@ -33,8 +33,7 @@ class CollisionConstraints { const Eigen::MatrixXd& vertices, const double dhat, const double dmin = 0, - const BroadPhaseMethod broad_phase_method = - BroadPhaseMethod::HASH_GRID); + const BroadPhaseMethod broad_phase_method = DEFAULT_BROAD_PHASE_METHOD); /// @brief Initialize the set of constraints used to compute the barrier potential. /// @param candidates Distance candidates from which the constraint set is built.