ipc-sim/ipc-toolkit

A set of reusable functions to integrate IPC into an existing simulation.

Description

IPC Toolkit is a set of reusable functions to integrate Incremental Potential Contact (IPC) into a simulation.

Features

• IPC barrier function and its derivatives and adaptive barrier stiffness algorithm
• Broad- and narrow-phase continuous collision detection (CCD) of linear and nonlinear trajectories
• Distance computation and derivatives between edges in 2D and triangles in 3D
• Distance barrier potential and its derivatives
• Smooth and lagged dissipative friction potential and its derivatives

Limitations

This is not a full simulation library. As such it does not include any physics or solvers. For a full simulation implementation, we recommend PolyFEM (a finite element library) or Rigid IPC (rigid-body dynamics) both of which utilize the IPC Toolkit.

Build

The easiest way to add the toolkit to an existing CMake project is to download it through CMake. CMake provides functionality for doing this called FetchContent (requires CMake ≥ 3.14). We use a very similar process to download all external dependencies (using CPM).

For example,

include(FetchContent)
FetchContent_Declare(
    ipc_toolkit
    GIT_REPOSITORY https://github.com/ipc-sim/ipc-toolkit.git
    GIT_TAG ${IPC_TOOLKIT_GIT_TAG}
)
FetchContent_MakeAvailable(ipc_toolkit)

where IPC_TOOLKIT_GIT_TAG is set to the version of the toolkit you want to use. This will download and add the toolkit to CMake. The toolkit can then be linked against using

# Link against the IPC Toolkit
target_link_libraries(${PROJECT_NAME} PUBLIC ipc::toolkit)

where PROJECT_NAME is the name of your library/binary.

Dependencies

All required dependencies are downloaded through CMake depending on the build options.

The following libraries are used in this project:

• Eigen: linear algebra
• libigl: basic geometry functions and predicates
• oneTBB: parallelism
• Tight-Inclusion: provably conservative CCD of [Wang and Ferguson et al. 2021]
• SimpleBVH: a simple bounding volume hierarchy data structure
• Scalable-CCD: scalable (GPU) CCD of [Belgrod et al. 2023]
• spdlog: logging information

Optional

The following dependencies are optionally used based on CMake options:

• robin-map: faster hash set/map than std::unordered_set/std::unordered_map
  • Enable by using the CMake option IPC_TOOLKIT_WITH_ROBIN_MAP
  • Enabled by default
• Abseil: hashing utilities
  • Enable by using the CMake option IPC_TOOLKIT_WITH_ABSEIL
  • Enabled by default
• filib: interval arithmetic for nonlinear trajectories/CCD
  • Enable by using the CMake option IPC_TOOLKIT_WITH_FILIB
  • Enabled by default
• rational-cpp: rational arithmetic used for exact intersection checks
  • Enable by using the CMake option IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION
  • Requires GMP to be installed at a system level
• Etienne Vouga's Collision Detection Library: inexact CCD
  • Included for comparison with the original IPC library
  • Enable by using the CMake option IPC_TOOLKIT_WITH_INEXACT_CCD
  • Replaces the default Tight-Inclusion CCD

Usage

See the tutorial for a quick introduction to the toolkit, or the documentation for a full reference.

Unit Tests

We provide unit tests to ensure the correctness of our algorithmic pieces. To enable the unit tests use the CMake option IPC_TOOLKIT_BUILD_TESTS.

Dependencies

The following are downloaded when unit tests are enabled:

• Catch2: testing framework
• finite-diff: finite-difference comparisons
• Nlohman's JSON library: loading test data from JSON files

Python Bindings

We provide Python bindings for functions in the toolkit using pybind11.

For more information see the Python documentation.

Contributing

This project is open to contributors! Contributions can come in the form of feature requests, bug fixes, documentation, tutorials, and the like. We highly recommend filing an Issue first before submitting a Pull Request.

Simply fork this repository and make a Pull Request! We would appreciate:

• Implementation of new features
• Bug Reports
• Documentation
• Testing

Citation

If you use the IPC Toolkit in your project, please consider citing our work:

@software{ipc_toolkit,
  author = {Zachary Ferguson and others},
  title = {{IPC Toolkit}},
  url = {https://github.com/ipc-sim/ipc-toolkit},
  year = {2020},
}

Additionally, you can cite the original IPC paper:

@article{Li2020IPC,
    author = {Minchen Li and Zachary Ferguson and Teseo Schneider and Timothy Langlois and
        Denis Zorin and Daniele Panozzo and Chenfanfu Jiang and Danny M. Kaufman},
    title = {Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics},
    journal = {{ACM} Trans. Graph. (SIGGRAPH)},
    year = {2020},
    volume = {39},
    number = {4},
    articleno = {49}
}

License

MIT License © 2020, the IPC-Sim organization (See LICENSE for details).

ipc-toolkit's Issues

Slightly confusing documentation for CollisionMesh

Following #27, I realized that there is one thing that has caused me more confusion than anything else. I understand the purpose of CollisionMesh, however some terminology is not quite clear to me.

For example, the CollisionConstraints::compute_potential_hessian method says:

...
@param[in] vertices Vertices of the collision mesh.
...
@returns The hessian of all barrier potentials (not scaled by the barrier stiffness). This will have a size of |vertices|x|vertices|.

This refers to vertices of the collision mesh. It's not quite clear to me if this generally refers to the "full mesh vertices" or the "surface mesh vertices". From some preliminary testing with CollisionMesh::build_from_full_mesh, it looks like it computes the Hessian with respect to the "full" mesh, is this correct? (Also, should that perhaps be 3|vertices| x 3 |vertices|, or something like that?)

However, in one of the CollisionMesh constructor, it also says:

@brief Construct a new Collision Mesh object directly from the collision mesh vertices.

So this sounds as if the collision mesh vertices here is something entirely different from the full mesh vertices, i.e. as opposed to full mesh vertices. Or perhaps, in this case, the full mesh vertices and the collision mesh vertices are anyway the same?

In short, the exact meaning of "vertices of the collision mesh" / "collision mesh vertices", as used throughout the library (perhaps in slightly different forms), is not entirely clear to me.

As I spend more time playing around with the library and through experimentation, what is expected and returned from various methods starts to become clearer, but for new users in particular I can imagine that this might be a source of confusion.

Possible Bug in the Definition of the Hash Structure

I'm learning the awesome work of IPC, and I really appreciate your generous contribution to this toolkit!

But when compiling the unordered_map_and_set.hpp (I don't set the abseil option, so the Hash structure is defined separately), some warning from GCC occurs, which says that the Hash::combine function is "returning reference to temporary". If I ignore the warning from compiler, segmentation fault happens when running into this line.

I'm not quite familiar with std::move or Rvalue reference, but I think there might be something wrong in this line, where you seem to return the reference of a temporary Hash variable.

Clarify dimensions of matrix inputs

Hi Zach and others! I thought I'd give the ipc-toolkit a spin for integrating contact with some deformable sims. Thanks for releasing the toolkit to the public!

I'm still getting used to the API, but one immediate point of confusion is that the matrix conventions are not so obvious to me. A number of functions take, for example, MatrixXd or MatrixXi, for example for describing a collision mesh. Of course, given the context, the options are limited: Edges for example are either 2xN or Nx2. However, since this doesn't appear to be documented (or is it? Maybe I missed it), working with the otherwise intuitive-looking APIs become a bit of guess work and/or detective work to reverse-engineer the conventions from the code. It would be great if this could be more clearly documented!

EDIT: In retrospect, I suppose that the library uses libigl conventions? I guess if you're not coming from the geometry processing community, these conventions may not be so familiar...

Control threading behavior

I am experiencing some non-determinism problems. This appears to be due to threading behavior. When running without contacts, it looks as if I am able to get deterministic (enough) output by disabling threading on my end. However, when running simulations with contacts with ipc-toolkit, I'm getting very non-deterministic behavior (seen, for example, through very varying number of solver iterations across runs with identical parameters).

I'd like to disable threading in ipc-toolkit to see if I get deterministic results this way. However, I am not so familiar with TBB and I'm not even sure if TBB is the only thing that's being used in ipc-toolkit, or if there are other sources of parallelism. Is there a way to disable parallelism or otherwise control the number of threads in ipc-toolkit?

Note: The problems I'm looking at are extremely sensitive (presumably due to near-singularity of matrices), I don't believe there's anything wrong with ipc-toolkit.

issue in distance_type.tpp

The return value of point_triangle_distance_type is incorrect.

For example: p =[0.488166 0.0132623 0.289055]; v0 = [0.456476 0.0526442 0.260834]; v1 = [0.609111 0.0595969 0.275928]; v2 = [0.431262 0.0508414 0.255831].

Questions regarding RigidIPC and GPU CCD

Hi there!

Thank you for making this great library open-source with such high-quality code!

I hope you don't mind me reaching out for some assistance. I was really excited to try out the GPU-based CCD implementation, and I was hoping to use it to speed up RigidIPC. However, I ran into some issues while updating the ipc-toolkit in RigidIPC, as there were a lot of changes since the last update. After resolving CMake errors, updating ipc-toolkit in RigidIPC gave me all sort of compilation errors (as shown in the screenshots).

I would really appreciate it if you could give me some pointers on how to resolve these issues. I understand that you might not have the time to fix them yourself, but any rough summary of the code changes or insights would be incredibly helpful for me.

Thank you!

Consider changing naming convention to indicate squared distances

Just spent some time scratching my head with CollisionConstraints::compute_minimum_distance. The documentation says it returns the minimum distance, but I found this in the implementation:

// NOTE: Actually distance squared

I realize that you're basically working with squared distances throughout (although for an outside user it is not clear if this holds in all cases), but it's still very confusing that "distance", which is a very unambiguously defined term in my opinion, actually means "squared distance". In my opinion, it would reduce confusion a lot for the new users and possible contributors if squared distances were always marked as such, e.g. by using naming like compute_squared_distance, compute_distance2 or similar.

(I labeled this as "bug" as, if you take the documentation and function names at their face value, what is returned is not what the naming suggests)

Co-dimensional simulations

Quick question: Does this toolkit include the ability to handle co-dimensional simulations similar to Co-IPC?

Shallow checkout with FetchContent may fail if tag is not HEAD

The suggested CMakeLists.txt configuration to fetch ipc-toolkit automatically is:

include(FetchContent)
FetchContent_Declare(
    ipc_toolkit
    GIT_REPOSITORY https://github.com/ipc-sim/ipc-toolkit.git
    GIT_TAG ${IPC_TOOLKIT_GIT_TAG}
    GIT_SHALLOW TRUE
)
FetchContent_MakeAvailable(ipc_toolkit)

However, this appears to only work reliably if GIT_TAG is set to a commit that corresponds to HEAD or possibly some other branch. For example, right now HEAD is 51cf348b50939913fa8f101845e7c55ac98c62f4, for which it works. However, if I set ${IPC_TOOLKIT_GIT_TAG} to c8bc2a1c13fa2c4627fc9a20c8e4bf02e5215164, which is one commit behind HEAD, then I get the following error:

-- The CXX compiler identification is GNU 11.3.0
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
[ 11%] Creating directories for 'ipc_toolkit-populate'
[ 22%] Performing download step (git clone) for 'ipc_toolkit-populate'
Cloning into 'ipc_toolkit-src'...
fatal: reference is not a tree: c8bc2a1c13fa2c4627fc9a20c8e4bf02e5215164
CMake Error at ipc_toolkit-subbuild/ipc_toolkit-populate-prefix/tmp/ipc_toolkit-populate-gitclone.cmake:40 (message):
  Failed to checkout tag: 'c8bc2a1c13fa2c4627fc9a20c8e4bf02e5215164'


gmake[2]: *** [CMakeFiles/ipc_toolkit-populate.dir/build.make:102: ipc_toolkit-populate-prefix/src/ipc_toolkit-populate-stamp/ipc_toolkit-populate-download] Error 1
gmake[1]: *** [CMakeFiles/Makefile2:83: CMakeFiles/ipc_toolkit-populate.dir/all] Error 2
gmake: *** [Makefile:91: all] Error 2

CMake Error at /usr/share/cmake-3.22/Modules/FetchContent.cmake:1087 (message):
  Build step for ipc_toolkit failed: 2
Call Stack (most recent call first):
  /usr/share/cmake-3.22/Modules/FetchContent.cmake:1216:EVAL:2 (__FetchContent_directPopulate)
  /usr/share/cmake-3.22/Modules/FetchContent.cmake:1216 (cmake_language)
  /usr/share/cmake-3.22/Modules/FetchContent.cmake:1259 (FetchContent_Populate)
  CMakeLists.txt:16 (FetchContent_MakeAvailable)


-- Configuring incomplete, errors occurred!

Similarly with older commits. It looks as if this might be a bug in cmake, related to the shallow clone functionality. Indeed, if I set GIT_SHALLOW FALSE it works also with other commits.

I'm not sure how to resolve this issue - since it appears to be a bug in cmake, perhaps just add a note in the documentation that it may be necessary to disable shallow cloning if you want to pin the toolkit to a particular commit? (which I guess you usually want to do in order to ensure that your build doesn't suddenly fail due to upstream changes)

Unexpectedly slow CCD

Hi, it's me again 👋

I noticed that CCD sometimes takes a very long time relative to other components of my simulator. For example, for a small impact problem, with a block dropping (with some initial velocity) on top of another, CCD can take ~20 seconds or more where factorization takes milliseconds. I can reproduce similar behavior for a very simple example, see this Gist here, where a single tetrahedron is squeezed so that the top vertex gets within dhat of the lower face. Here I'm getting timings like this:

Potential duration: 6.1011e-05 seconds.
CCD duration: 2.25247 seconds.

Here's the most crucial part of the code (see the above gist for the full setup):

    // Run this several times just to demonstrate that it's not due to some initialization or similar
    for (int i = 0; i < 5; ++i) {
        const auto potential_begin = std::chrono::steady_clock::now();
        const double potential = constraints.compute_potential(mesh, deformed_vertices, dhat);
        const Eigen::VectorXd grad = constraints.compute_potential_gradient(mesh, deformed_vertices, dhat);
        const Eigen::SparseMatrix<double> hessian = constraints.compute_potential_hessian(mesh, deformed_vertices, dhat,project_spd);
        const auto potential_end = std::chrono::steady_clock::now();
        const auto potential_duration = std::chrono::duration<double>(potential_end - potential_begin).count();

        const auto ccd_begin = std::chrono::steady_clock::now();
        const double alpha = ipc::compute_collision_free_stepsize(mesh, rest_vertices, deformed_vertices);
        const auto ccd_end = std::chrono::steady_clock::now();
        const auto ccd_duration = std::chrono::duration<double>(ccd_end - ccd_begin).count();

        std::cout << "Potential duration: " << potential_duration << " seconds." << std::endl;
        std::cout << "CCD duration: " << ccd_duration << " seconds." << std::endl;
    }

Build type is Release. I don't have all that much experience with CCD, but this seems slower than what I was expecting. Is there something I should configure that I haven't? I've basically just followed the docs/tutorial here. Any hints as to what I can do here, if anything, would be much appreciated :-)

std::array issue

I use your codes as suggested. But encountered with error: implicit instantiation of undefined template 'std::array<long, 4>' for collision_constraint.hpp.

Fix: add #include inside the collision_constraint.hpp.

Basic example of initialization and simulation loop

Hi, first of all, thank you for sharing this implementation of IPC.
Would it be possible to have a complete example showing the initialization of deformable bodies and the simulation loop?
Is there any documentation on how to use this toolkit?

Possibly inconsistent results for convergent formulation with respect to dhat?

Describe the bug
I'm using CollisionConstraints::set_use_convergent_formulation(true) to obtain quantities that are less tessellation-dependent and that - to my understanding - vary less with the value of dhat.

I am not sure if perhaps I've misunderstood the true purpose of this formulation, but I'm observing inconsistent values for different values of dhat in a setting in which I would expect the values to be more similar.

Consider the following example. Take a squished unit tetrahedron with vertices (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 0.5 * dhat). Since the z-coordinate of the "top" vertex is proportional to dhat, I would expect a similar potential energy value for different values of dhat (noting of course that changing dhat for this example also changes the area of the triangles).

Here's the output I get for different values of dhat for this example:

dhat = 0.0001
  Potential: 5.26234e-09
  Scaled potential: 0.526234
dhat = 0.001
  Potential: 5.26234e-07
  Scaled potential: 0.526234
dhat = 0.01
  Potential: 5.26251e-05
  Scaled potential: 0.526251

Since the inputs to the barrier function are squared, I believe the "real" barrier contribution works out to -log(d^2/dhat^2) (d^2 - dhat^2)^2. To make this unitless, we'd have to divide by dhat^4, however my understanding is that the weight constant currently used to scale the barrier function is divided by dhat^2 (when dmin = 0). The Scaled potential in the above output refers to dividing the potential by dhat^2, to make up for the missing powers.

To Reproduce
Here's a self-contained minimal example that I used to produce the output above:

#include <ipc/ipc.hpp>
#include <iostream>
#include <array>

int main() {
    const auto dhat_values = std::array<double, 3> { 1e-4, 1e-3, 1e-2 };

    for (const auto dhat : dhat_values) {
        Eigen::MatrixXd rest_vertices(4, 3);
        rest_vertices <<
                      0.0, 0.0, 0.0,
                1.0, 0.0, 0.0,
                0.0, 1.0, 0.0,
                0.0, 0.0, 1.0;

        Eigen::MatrixXd deformed_vertices = rest_vertices;
        deformed_vertices.row(3)(2) = 0.5 * dhat;

        Eigen::MatrixXi edges(6, 2);
        edges << 0, 1,
                0, 2,
                0, 3,
                1, 2,
                1, 3,
                2, 3;
        Eigen::MatrixXi faces(4, 3);
        faces << 0, 2, 1,
                0, 1, 3,
                0, 3, 2,
                1, 2, 3;

        ipc::CollisionMesh mesh = ipc::CollisionMesh::build_from_full_mesh(rest_vertices, edges, faces);
        ipc::CollisionConstraints constraints;
        constraints.set_use_convergent_formulation(true);
        constraints.build(mesh, deformed_vertices, dhat);

        const double potential = constraints.compute_potential(mesh, deformed_vertices, dhat);
        std::cout << "dhat = " << dhat << std::endl;
        std::cout << "  Potential: " << potential << std::endl;
        std::cout << "  Scaled potential: " << potential / (dhat * dhat) << std::endl;
    }

    return 0;
}

Expected behavior
I had expected to get the value of the "scaled potential" above when using the "convergent formulation". This might however be due to a misunderstanding on my part.

Additional context
This is not a big deal on my part, because I can work around it by dividing by dhat^2, but I thought I ought to let you know, in case this is an issue that has gone undiscovered.

EdgeEdgeConstraint distance gradient crash

The method EdgeEdgeConstraint::compute_distance_gradient calls edge_edge_distance_gradient with VectorMax9d as its output, and I'm hitting an Eigen runtime size assert. Shouldn't this vector be VectorMax12d? Switching to VectorMax12d fixes it on my end, but I wanted to make sure.

If this is an unexpected bug I can post more details, thanks :)

Closed-form expressions for the Distance Gradients / Hessians

Hi,
I'm implementing IPC in my advisor's physics library. However, the paper does not appear to provide closed-form expressions for the {gradients, Hessians} of the {point-plane, point-line, line-line} distance functions WRT all inputs. Instead, the toolkit only has mangled, auto-generated functions. Since my advisor wants to make sure that I fully understand the math, I need to know the un-mangled expressions for these functions. So far, I've managed to reverse-engineer the following functions below. Could you please provide me the un-mangled forms for all the mangled functions? I'm relatively new to this math (e.g. I have never tried deriving something like the gradient of point-plane distance WRT the plane itself).

Thank you for your time.

VECTOR<T, 12>
Point_Plane_Distance_Gradient(const TV& p, const TV& a, const TV& b, const TV& c)
{
    TV pa = p - a;
    TV ba = b - a;
    TV ca = c - a;
    TV cb = c - b;
    TV normal = ba.Cross(ca);
    T one_over_normal_sq = 1 / normal.Magnitude_Squared();
    T one_over_normal_sq_sq = one_over_normal_sq * one_over_normal_sq;
    T p_distance = pa.Dot(normal);
    T p_distance_sq = p_distance * p_distance;

    return 2 * p_distance * one_over_normal_sq * VECTOR<T, 12>(
        normal,                                                                     //g_p
        VECTOR<T, 9>(
            p_distance * one_over_normal_sq * cb.Cross(normal) - (normal + pa.Cross(cb)), //g_a
            VECTOR<T, 6>(
                ca.Cross(normal) + pa.Cross(ca) * p_distance * one_over_normal_sq,  //g_b
                pa.Cross(ba) - ba.Cross(normal) * p_distance * one_over_normal_sq   //g_c
            )
        )
    );
}

VECTOR<T, 9>
Point_Line_Distance_Gradient(const TV& p, const TV& a, const TV& b) const
{
    TV pa = p - a;
    TV pb = p - b;
    TV ab = a - b;
    T over_ab_sq = 1.0 / ab.Magnitude_Squared();
    TV paXpb = (p - a).Cross(p - b);
    T over_ab_sq_sq = over_ab_sq * over_ab_sq;
    T paXpb_sq = paXpb.Magnitude_Squared();
    TV f = 2 * ab * paXpb_sq * over_ab_sq_sq;
    return 2 * over_ab_sq * VECTOR<T, 9>(
        ab.Cross(paXpb),
        VECTOR<T, 6>(
            f + pb.Cross(paXpb),
            f + pa.Cross(paXpb)
        )
    );

Average body mass in IPC-toolkit vs average lumped nodal mass in IPC

I have a question about the arguments in adaptive barrier stiffness function called ipc::initial_barrier_stiffness, it has one argument called const double average_mass which indicates as average mass of all bodies while in IPC function called suggestKappa(kappa) and according to Supplement of IPC paper, it should be average of lumped nodal mass. Does it the same as IPC?
And what types of CCD IPC-toolkit does use for ipc::compute_collision_free_stepsize, it seems that it is TIGHT_INCLUSION while it seems that IPC uses FLOATING_POINT_ROOT_FINDER as default. Can I use this with IPC-toolkit?