Geometric Contact Potential

CMakeLists.txt

@@ -77,6 +77,8 @@ else()
   set(IPC_TOOLKIT_BUILD_PYTHON OFF CACHE BOOL "Build Python bindings" FORCE)
 endif()
 
+set(IPC_TOOLKIT_TESTS_GCP_DATA_DIR "${PROJECT_SOURCE_DIR}/tests/gcp-data/" CACHE PATH "")
+
 if(CMAKE_CUDA_COMPILER)
   option(IPC_TOOLKIT_WITH_CUDA "Enable CUDA CCD" ON)
 else()
@@ -91,6 +93,7 @@ option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD meth
 # Advanced options
 option(IPC_TOOLKIT_WITH_SIMD                  "Enable SIMD"                                   OFF)
 option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
+option(IPC_TOOLKIT_WITH_AUTODIFF              "Enable AutoDiff for computing derivatives"     OFF)
 
 mark_as_advanced(IPC_TOOLKIT_WITH_SIMD)          # This does not work reliably
 mark_as_advanced(IPC_TOOLKIT_WITH_CODE_COVERAGE) # This is used in GitHub Actions
@@ -235,6 +238,11 @@ endif()
 # Use C++17
 target_compile_features(ipc_toolkit PUBLIC cxx_std_17)
 
+if(IPC_TOOLKIT_WITH_AUTODIFF)
+  target_compile_definitions(ipc_toolkit PUBLIC DERIVATIVES_WITH_AUTODIFF)
+endif()
+target_compile_definitions(ipc_toolkit PUBLIC IPC_TOOLKIT_WITH_TBB)
+
 ################################################################################
 # CUDA
 ################################################################################
@@ -306,4 +314,4 @@ if (CMAKE_GENERATOR STREQUAL "Xcode")
   if(IPC_TOOLKIT_BUILD_PYTHON)
     set_target_properties(ipctk PROPERTIES XCODE_GENERATE_SCHEME ON)
   endif()
-endif()
+endif()

cmake/recipes/autodiff.cmake

@@ -0,0 +1,16 @@
+# autodiff (https://github.com/autodiff/autodiff.git)
+# BSD license
+if(TARGET autodiff::autodiff)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'autodiff::autodiff'")
+
+option(AUTODIFF_BUILD_TESTS "Enable the compilation of the test files." OFF)
+option(AUTODIFF_BUILD_PYTHON "Enable the compilation of the python bindings." OFF)
+option(AUTODIFF_BUILD_EXAMPLES "Enable the compilation of the example files." OFF)
+option(AUTODIFF_BUILD_DOCS "Enable the build of the documentation and website." OFF)
+
+include(CPM)
+CPMAddPackage("gh:autodiff/[email protected]")
+target_link_libraries(ipc_toolkit PRIVATE autodiff::autodiff)

python/src/bindings.cpp

@@ -88,6 +88,8 @@ PYBIND11_MODULE(ipctk, m)
     define_smooth_friction_mollifier(m);
     define_smooth_mu(m);
 
+    define_smooth_potential(m);
+
     // implicits
     define_plane_implicit(m);
 

python/src/collision_mesh.cpp

@@ -151,7 +151,8 @@ void define_collision_mesh(py::module_& m)
             "displacement_map"_a = Eigen::SparseMatrix<double>())
         .def(
             py::init<
-                const std::vector<bool>&, Eigen::ConstRef<Eigen::MatrixXd>,
+                const std::vector<bool>&, const std::vector<bool>&,
+                Eigen::ConstRef<Eigen::MatrixXd>,
                 Eigen::ConstRef<Eigen::MatrixXi>,
                 Eigen::ConstRef<Eigen::MatrixXi>,
                 const Eigen::SparseMatrix<double>&>(),
@@ -160,12 +161,13 @@ void define_collision_mesh(py::module_& m)
 
             Parameters:
                 include_vertex: Vector of bools indicating whether each vertex should be included in the collision mesh.
+                orient_vertex: Vector of bools indicating whether each vertex is orientable.
                 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",
-            "include_vertex"_a, "full_rest_positions"_a,
+            "include_vertex"_a, "orient_vertex"_a, "full_rest_positions"_a,
             "edges"_a = Eigen::MatrixXi(), "faces"_a = Eigen::MatrixXi(),
             "displacement_map"_a = Eigen::SparseMatrix<double>())
         .def_static(

python/src/collisions/normal/normal_collisions.cpp

@@ -1,9 +1,86 @@
 #include <common.hpp>
 
 #include <ipc/collisions/normal/normal_collisions.hpp>
+#include <ipc/smooth_contact/smooth_collisions.hpp>
 
 using namespace ipc;
 
+template <typename collision_type, typename parent_type>
+void define_SmoothCollisionTemplate(py::module_& m, std::string name)
+{
+    py::class_<collision_type, parent_type>(m, name.c_str())
+        .def("name", &collision_type::name, "Get the type name of collision")
+        .def(
+            "num_vertices", &collision_type::num_vertices,
+            "Get the number of vertices");
+}
+
+void define_SmoothCollisions(py::module_& m, std::string name)
+{
+    py::class_<SmoothCollisions>(m, name.c_str())
+        .def(py::init())
+        .def(
+            "build",
+            py::overload_cast<
+                const CollisionMesh&, Eigen::ConstRef<Eigen::MatrixXd>,
+                const ParameterType, const bool, std::shared_ptr<BroadPhase>>(
+                &SmoothCollisions::build),
+            R"ipc_Qu8mg5v7(
+            Initialize the set of collisions used to compute the barrier potential.
+
+            Parameters:
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+                param: ParameterType.
+                use_adaptive_dhat: If the adaptive dhat should be used.
+                broad_phase: Broad phase method.
+            )ipc_Qu8mg5v7",
+            py::arg("mesh"), py::arg("vertices"), py::arg("param"),
+            py::arg("use_adaptive_dhat") = false,
+            py::arg("broad_phase") = make_default_broad_phase())
+        .def(
+            "compute_minimum_distance",
+            &SmoothCollisions::compute_minimum_distance,
+            R"ipc_Qu8mg5v7(
+            Computes the minimum distance between any non-adjacent elements.
+
+            Parameters:
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+
+            Returns:
+                The minimum distance between any non-adjacent elements.
+            )ipc_Qu8mg5v7",
+            py::arg("mesh"), py::arg("vertices"))
+        .def(
+            "__len__", &SmoothCollisions::size, "Get the number of collisions.")
+        .def(
+            "empty", &SmoothCollisions::empty,
+            "Get if the collision set is empty.")
+        .def("clear", &SmoothCollisions::clear, "Clear the collision set.")
+        .def(
+            "__getitem__",
+            [](SmoothCollisions& self, size_t i) ->
+            typename SmoothCollisions::value_type& { return self[i]; },
+            py::return_value_policy::reference,
+            R"ipc_Qu8mg5v7(
+            Get a reference to collision at index i.
+
+            Parameters:
+                i: The index of the collision.
+
+            Returns:
+                A reference to the collision.
+            )ipc_Qu8mg5v7",
+            py::arg("i"))
+        .def(
+            "to_string", &SmoothCollisions::to_string, py::arg("mesh"),
+            py::arg("vertices"), py::arg("param"))
+        .def(
+            "n_candidates", &SmoothCollisions::n_candidates,
+            "Get the number of candidates.");
+}
+
 void define_normal_collisions(py::module_& m)
 {
     py::class_<NormalCollisions>(m, "NormalCollisions")
@@ -159,4 +236,42 @@ void define_normal_collisions(py::module_& m)
         .def_readwrite("ee_collisions", &NormalCollisions::ee_collisions)
         .def_readwrite("fv_collisions", &NormalCollisions::fv_collisions)
         .def_readwrite("pv_collisions", &NormalCollisions::pv_collisions);
+
+    py::class_<SmoothCollision>(m, "SmoothCollision2")
+        .def("n_dofs", &SmoothCollision::n_dofs, "Get the degree of freedom")
+        .def(
+            "__call__", &SmoothCollision::operator(),
+            R"ipc_Qu8mg5v7(
+            Compute the potential.
+
+            Parameters:
+                positions: The vertex positions.
+                params: The parameters.
+
+            Returns:
+                The potential (not scaled by the barrier stiffness) of this collision pair.
+            )ipc_Qu8mg5v7",
+            py::arg("positions"), py::arg("params"))
+        .def(
+            "__getitem__",
+            [](SmoothCollision& self, size_t i) -> long { return self[i]; },
+            R"ipc_Qu8mg5v7(
+            Get primitive id.
+
+            Parameters:
+                i: 0 or 1.
+
+            Returns:
+                The index of the primitive.
+            )ipc_Qu8mg5v7",
+            py::arg("i"));
+
+    define_SmoothCollisionTemplate<
+        SmoothCollisionTemplate<Edge2, Point2>, SmoothCollision>(
+        m, "Edge2Point2Collision");
+    define_SmoothCollisionTemplate<
+        SmoothCollisionTemplate<Point2, Point2>, SmoothCollision>(
+        m, "Point2Point2Collision");
+
+    define_SmoothCollisions(m, "SmoothCollisions");
 }

python/src/potentials/barrier_potential.cpp

@@ -1,6 +1,7 @@
 #include <common.hpp>
 
 #include <ipc/potentials/barrier_potential.hpp>
+#include <ipc/smooth_contact/smooth_contact_potential.hpp>
 
 using namespace ipc;
 
@@ -37,3 +38,153 @@ void define_barrier_potential(py::module_& m)
             &BarrierPotential::set_barrier,
             "Barrier function used to compute the potential.");
 }
+
+void define_smooth_potential(py::module_& m)
+{
+    py::class_<ParameterType>(m, "ParameterType")
+        .def(
+            py::init<
+                const double&, const double&, const double&, const double&,
+                const double&, const int&>(),
+            R"ipc_Qu8mg5v7(
+            Construct parameter set for smooth contact.
+
+            Parameters:
+                dhat, alpha_t, beta_t, alpha_n, beta_n, r
+            )ipc_Qu8mg5v7",
+            py::arg("dhat"), py::arg("alpha_t"), py::arg("beta_t"),
+            py::arg("alpha_n"), py::arg("beta_n"), py::arg("r"))
+        .def(
+            py::init<const double&, const double&, const double&, const int&>(),
+            R"ipc_Qu8mg5v7(
+            Construct parameter set for smooth contact.
+
+            Parameters:
+                dhat, alpha_t, beta_t, r
+            )ipc_Qu8mg5v7",
+            py::arg("dhat"), py::arg("alpha_t"), py::arg("beta_t"),
+            py::arg("r"))
+        .def_readonly("dhat", &ParameterType::dhat)
+        .def_readonly("alpha_t", &ParameterType::alpha_t)
+        .def_readonly("beta_t", &ParameterType::beta_t)
+        .def_readonly("alpha_n", &ParameterType::alpha_n)
+        .def_readonly("beta_n", &ParameterType::beta_n)
+        .def_readonly("r", &ParameterType::r);
+
+    py::class_<SmoothContactPotential>(m, "SmoothPotential")
+        .def(
+            py::init<const ParameterType&>(),
+            R"ipc_Qu8mg5v7(
+            Construct a smooth barrier potential.
+
+            Parameters:
+                param: A set of parameters.
+            )ipc_Qu8mg5v7",
+            py::arg("param"))
+        .def(
+            "__call__",
+            py::overload_cast<
+                const SmoothCollisions&, const CollisionMesh&,
+                Eigen::ConstRef<Eigen::MatrixXd>>(
+                &ipc::SmoothContactPotential::operator(), py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the barrier potential for a set of collisions.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+
+            Returns:
+                The sum of all barrier potentials (not scaled by the barrier stiffness).
+            )ipc_Qu8mg5v7",
+            py::arg("collisions"), py::arg("mesh"), py::arg("vertices"))
+        .def(
+            "gradient",
+            py::overload_cast<
+                const SmoothCollisions&, const CollisionMesh&,
+                Eigen::ConstRef<Eigen::MatrixXd>>(
+                &ipc::SmoothContactPotential::gradient, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the gradient of the barrier potential.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+
+            Returns:
+                The gradient of all barrier potentials (not scaled by the barrier stiffness). This will have a size of |vertices|.
+            )ipc_Qu8mg5v7",
+            py::arg("collisions"), py::arg("mesh"), py::arg("vertices"))
+        .def(
+            "hessian",
+            py::overload_cast<
+                const SmoothCollisions&, const CollisionMesh&,
+                Eigen::ConstRef<Eigen::MatrixXd>, const PSDProjectionMethod>(
+                &ipc::SmoothContactPotential::hessian, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the hessian of the barrier potential.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+                project_hessian_to_psd: Make sure the hessian is positive semi-definite.
+
+            Returns:
+                The hessian of all barrier potentials (not scaled by the barrier stiffness). This will have a size of |vertices|x|vertices|.
+            )ipc_Qu8mg5v7",
+            py::arg("collisions"), py::arg("mesh"), py::arg("vertices"),
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE)
+        .def(
+            "__call__",
+            py::overload_cast<
+                const SmoothCollision&, Eigen::ConstRef<Eigen::VectorXd>>(
+                &ipc::SmoothContactPotential::operator(), py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the potential for a single collision.
+
+            Parameters:
+                collision: The collision.
+                x: The collision stencil's degrees of freedom.
+
+            Returns:
+                The potential.
+            )ipc_Qu8mg5v7",
+            py::arg("collision"), py::arg("x"))
+        .def(
+            "gradient",
+            py::overload_cast<
+                const SmoothCollision&, Eigen::ConstRef<Eigen::VectorXd>>(
+                &SmoothContactPotential::gradient, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the gradient of the potential for a single collision.
+
+            Parameters:
+                collision: The collision.
+                x: The collision stencil's degrees of freedom.
+
+            Returns:
+                The gradient of the potential.
+            )ipc_Qu8mg5v7",
+            py::arg("collision"), py::arg("x"))
+        .def(
+            "hessian",
+            py::overload_cast<
+                const SmoothCollision&, Eigen::ConstRef<Eigen::VectorXd>,
+                const PSDProjectionMethod>(
+                &SmoothContactPotential::hessian, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the hessian of the potential for a single collision.
+
+            Parameters:
+                collision: The collision.
+                x: The collision stencil's degrees of freedom.
+
+            Returns:
+                The hessian of the potential.
+            )ipc_Qu8mg5v7",
+            py::arg("collision"), py::arg("x"),
+            py::arg("project_hessian_to_psd") = PSDProjectionMethod::NONE);
+}

python/src/potentials/bindings.hpp

@@ -3,6 +3,7 @@
 #include <pybind11/pybind11.h>
 
 void define_barrier_potential(py::module_& m);
+void define_smooth_potential(py::module_& m);
 void define_friction_potential(py::module_& m);
 void define_normal_adhesion_potential(py::module_& m);
 void define_normal_potential(py::module_& m);

src/ipc/CMakeLists.txt

@@ -22,4 +22,5 @@ add_subdirectory(friction)
 add_subdirectory(implicits)
 add_subdirectory(potentials)
 add_subdirectory(tangent)
-add_subdirectory(utils)
+add_subdirectory(utils)
+add_subdirectory(smooth_contact)