From 4d05d6b9c041bd602cdeb5bbbf6d7b734bdec745 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 29 Oct 2025 13:14:48 -0400 Subject: [PATCH 1/2] Organize geometry and math utilities --- python/src/CMakeLists.txt | 2 + python/src/bindings.cpp | 12 +- python/src/bindings.hpp | 2 + python/src/geometry/CMakeLists.txt | 8 + python/src/geometry/angle.cpp | 60 ++++ .../area_gradient.cpp => geometry/area.cpp} | 4 +- python/src/geometry/bindings.hpp | 8 + .../src/{utils => geometry}/intersection.cpp | 2 +- python/src/geometry/normal.cpp | 306 ++++++++++++++++++ python/src/math/CMakeLists.txt | 5 + python/src/math/bindings.hpp | 5 + python/src/{utils => math}/interval.cpp | 2 +- python/src/utils/CMakeLists.txt | 3 - python/src/utils/bindings.hpp | 3 - src/ipc/CMakeLists.txt | 1 + src/ipc/ccd/nonlinear_ccd.hpp | 2 +- src/ipc/collision_mesh.cpp | 24 +- src/ipc/geometry/CMakeLists.txt | 6 + src/ipc/geometry/angle.cpp | 56 ++++ src/ipc/geometry/angle.hpp | 37 +++ .../area_gradient.cpp => geometry/area.cpp} | 2 +- .../area_gradient.hpp => geometry/area.hpp} | 23 ++ src/ipc/{utils => geometry}/intersection.cpp | 0 src/ipc/{utils => geometry}/intersection.hpp | 0 src/ipc/geometry/normal.cpp | 40 +++ src/ipc/geometry/normal.hpp | 20 +- src/ipc/ipc.cpp | 2 +- src/ipc/math/CMakeLists.txt | 9 + src/ipc/{utils => math}/interval.cpp | 0 src/ipc/{utils => math}/interval.hpp | 0 src/ipc/{utils => math}/math.cpp | 0 src/ipc/{utils => math}/math.hpp | 0 src/ipc/{utils => math}/math.tpp | 0 .../smooth_contact/distance/point_edge.hpp | 2 +- src/ipc/utils/CMakeLists.txt | 9 - tests/src/tests/CMakeLists.txt | 2 + tests/src/tests/barrier/test_barrier.cpp | 2 +- tests/src/tests/candidates/test_normals.cpp | 39 +++ tests/src/tests/distance/test_point_point.cpp | 2 +- tests/src/tests/geometry/CMakeLists.txt | 9 + tests/src/tests/geometry/test_angle.cpp | 59 ++++ tests/src/tests/math/CMakeLists.txt | 14 + .../tests/{utils => math}/test_interval.cpp | 2 +- tests/src/tests/utils/CMakeLists.txt | 1 - 44 files changed, 738 insertions(+), 47 deletions(-) create mode 100644 python/src/geometry/CMakeLists.txt create mode 100644 python/src/geometry/angle.cpp rename python/src/{utils/area_gradient.cpp => geometry/area.cpp} (91%) create mode 100644 python/src/geometry/bindings.hpp rename python/src/{utils => geometry}/intersection.cpp (96%) create mode 100644 python/src/geometry/normal.cpp create mode 100644 python/src/math/CMakeLists.txt create mode 100644 python/src/math/bindings.hpp rename python/src/{utils => math}/interval.cpp (99%) create mode 100644 src/ipc/geometry/angle.cpp create mode 100644 src/ipc/geometry/angle.hpp rename src/ipc/{utils/area_gradient.cpp => geometry/area.cpp} (98%) rename src/ipc/{utils/area_gradient.hpp => geometry/area.hpp} (61%) rename src/ipc/{utils => geometry}/intersection.cpp (100%) rename src/ipc/{utils => geometry}/intersection.hpp (100%) create mode 100644 src/ipc/math/CMakeLists.txt rename src/ipc/{utils => math}/interval.cpp (100%) rename src/ipc/{utils => math}/interval.hpp (100%) rename src/ipc/{utils => math}/math.cpp (100%) rename src/ipc/{utils => math}/math.hpp (100%) rename src/ipc/{utils => math}/math.tpp (100%) create mode 100644 tests/src/tests/geometry/CMakeLists.txt create mode 100644 tests/src/tests/geometry/test_angle.cpp create mode 100644 tests/src/tests/math/CMakeLists.txt rename tests/src/tests/{utils => math}/test_interval.cpp (99%) diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt index 5e8499e29..5a9b5f9d9 100644 --- a/python/src/CMakeLists.txt +++ b/python/src/CMakeLists.txt @@ -18,7 +18,9 @@ add_subdirectory(ccd) add_subdirectory(collisions) add_subdirectory(distance) add_subdirectory(friction) +add_subdirectory(geometry) add_subdirectory(implicits) +add_subdirectory(math) add_subdirectory(potentials) add_subdirectory(tangent) add_subdirectory(utils) \ No newline at end of file diff --git a/python/src/bindings.cpp b/python/src/bindings.cpp index 8bffaf376..479bc5255 100644 --- a/python/src/bindings.cpp +++ b/python/src/bindings.cpp @@ -90,9 +90,18 @@ PYBIND11_MODULE(ipctk, m) define_smooth_potential(m); + // geometry + define_angle(m); + define_area(m); + define_intersection(m); + define_normal(m); + // implicits define_plane_implicit(m); + // math + define_interval(m); + // potentials define_normal_potential(m); // define early because it is used next define_barrier_potential(m); @@ -102,9 +111,6 @@ PYBIND11_MODULE(ipctk, m) define_tangential_adhesion_potential(m); // utils - define_area_gradient(m); - define_interval(m); - define_intersection(m); define_logger(m); define_thread_limiter(m); define_vertex_to_min_edge(m); diff --git a/python/src/bindings.hpp b/python/src/bindings.hpp index a680b4f14..66db6b2b2 100644 --- a/python/src/bindings.hpp +++ b/python/src/bindings.hpp @@ -8,7 +8,9 @@ #include #include #include +#include #include +#include #include #include #include diff --git a/python/src/geometry/CMakeLists.txt b/python/src/geometry/CMakeLists.txt new file mode 100644 index 000000000..68c6738a5 --- /dev/null +++ b/python/src/geometry/CMakeLists.txt @@ -0,0 +1,8 @@ +set(SOURCES + area.cpp + angle.cpp + normal.cpp + intersection.cpp +) + +target_sources(ipctk PRIVATE ${SOURCES}) \ No newline at end of file diff --git a/python/src/geometry/angle.cpp b/python/src/geometry/angle.cpp new file mode 100644 index 000000000..39e017f97 --- /dev/null +++ b/python/src/geometry/angle.cpp @@ -0,0 +1,60 @@ +#include + +#include + +using namespace ipc; + +void define_angle(py::module_& m) +{ + m.def( + "dihedral_angle", &dihedral_angle, + R"ipc_Qu8mg5v7( + Compute the bending angle between two triangles sharing an edge. + x0---x2 + | \ | + x1---x3 + + Parameters + ---------- + x0 : Eigen::Vector3d + The first vertex of the edge. + x1 : Eigen::Vector3d + The second vertex of the edge. + x2 : Eigen::Vector3d + The opposite vertex of the first triangle. + x3 : Eigen::Vector3d + The opposite vertex of the second triangle. + + Returns + ------- + double + The bending angle between the two triangles. + )ipc_Qu8mg5v7", + py::arg("x0"), py::arg("x1"), py::arg("x2"), py::arg("x3")); + + m.def( + "dihedral_angle_gradient", &dihedral_angle_gradient, + R"ipc_Qu8mg5v7( + Compute the Jacobian of the bending angle between two triangles sharing an edge. + x0---x2 + | \ | + x1---x3 + + Parameters + ---------- + x0 : Eigen::Vector3d + The first vertex of the edge. + x1 : Eigen::Vector3d + The second vertex of the edge. + x2 : Eigen::Vector3d + The opposite vertex of the first triangle. + x3 : Eigen::Vector3d + The opposite vertex of the second triangle. + + Returns + ------- + Eigen::Vector + The Jacobian matrix of the bending angle with respect to the input vertices. + )ipc_Qu8mg5v7", + py::arg("x0"), py::arg("x1"), py::arg("x2"), py::arg("x3")); +} diff --git a/python/src/utils/area_gradient.cpp b/python/src/geometry/area.cpp similarity index 91% rename from python/src/utils/area_gradient.cpp rename to python/src/geometry/area.cpp index df6152676..ff917c2cb 100644 --- a/python/src/utils/area_gradient.cpp +++ b/python/src/geometry/area.cpp @@ -1,11 +1,11 @@ #include -#include +#include #include using namespace ipc; -void define_area_gradient(py::module_& m) +void define_area(py::module_& m) { m.def( "edge_length_gradient", &edge_length_gradient, diff --git a/python/src/geometry/bindings.hpp b/python/src/geometry/bindings.hpp new file mode 100644 index 000000000..7ee2455c9 --- /dev/null +++ b/python/src/geometry/bindings.hpp @@ -0,0 +1,8 @@ +#pragma once + +#include + +void define_angle(py::module_& m); +void define_area(py::module_& m); +void define_normal(py::module_& m); +void define_intersection(py::module_& m); \ No newline at end of file diff --git a/python/src/utils/intersection.cpp b/python/src/geometry/intersection.cpp similarity index 96% rename from python/src/utils/intersection.cpp rename to python/src/geometry/intersection.cpp index c30d25286..946568f81 100644 --- a/python/src/utils/intersection.cpp +++ b/python/src/geometry/intersection.cpp @@ -1,6 +1,6 @@ #include -#include +#include #include diff --git a/python/src/geometry/normal.cpp b/python/src/geometry/normal.cpp new file mode 100644 index 000000000..52bfe2885 --- /dev/null +++ b/python/src/geometry/normal.cpp @@ -0,0 +1,306 @@ +#include + +#include + +using namespace ipc; + +void define_normal(py::module_& m) +{ + m.def( + "normalization_and_jacobian", &normalization_and_jacobian, + R"ipc_qu8mg5v7( + Computes the normalization and Jacobian of a vector. + + Parameters + ---------- + x: The input vector. + + Returns + ------- + A tuple containing the normalized vector and its Jacobian. + )ipc_qu8mg5v7", + py::arg("x")); + + m.def( + "normalization_jacobian", &normalization_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the normalization operation. + + Parameters + ---------- + x: The input vector. + + Returns + ------- + The Jacobian of the normalization operation. + )ipc_qu8mg5v7", + py::arg("x")); + + m.def( + "normalization_and_jacobian_and_hessian", + &normalization_and_jacobian_and_hessian, + R"ipc_qu8mg5v7( + Computes the normalization, Jacobian, and Hessian of a vector. + + Parameters + ---------- + x: The input vector. + + Returns + ------- + A tuple containing the normalized vector, its Jacobian, and its Hessian. + )ipc_qu8mg5v7", + py::arg("x")); + + m.def( + "normalization_hessian", &normalization_hessian, + R"ipc_qu8mg5v7( + Computes the Hessian of the normalization operation. + + Parameters + ---------- + x: The input vector. + + Returns + ------- + The Hessian of the normalization operation. + )ipc_qu8mg5v7", + py::arg("x")); + + m.def( + "cross_product_matrix", &cross_product_matrix, + R"ipc_qu8mg5v7( + Cross product matrix for 3D vectors. + + Parameters + ---------- + v: Vector to create the cross product matrix for. + + Returns + ------- + The cross product matrix of the vector. + )ipc_qu8mg5v7", + py::arg("v")); + + m.def( + "cross_product_matrix_jacobian", &cross_product_matrix_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the cross product matrix. + Returns + ------- + The Jacobian of the cross product matrix. + )ipc_qu8mg5v7"); + + m.def( + "edge_vertex_unnormalized_normal", &edge_vertex_unnormalized_normal, + R"ipc_qu8mg5v7( + Computes the unnormalized normal vector of an edge-vertex pair. + + Parameters + ---------- + v: The vertex position. + e0: The start position of the edge. + e1: The end position of the edge. + + Returns + ------- + The unnormalized normal vector. + )ipc_qu8mg5v7", + py::arg("v"), py::arg("e0"), py::arg("e1")); + + m.def( + "edge_vertex_normal", &edge_vertex_normal, + R"ipc_qu8mg5v7( + Computes the normal vector of an edge-vertex pair. + + Parameters + ---------- + v: The vertex position. + e0: The start position of the edge. + e1: The end position of the edge. + + Returns + ------- + The normal vector. + The normal vector. + )ipc_qu8mg5v7", + py::arg("v"), py::arg("e0"), py::arg("e1")); + + m.def( + "edge_vertex_unnormalized_normal_jacobian", + &edge_vertex_unnormalized_normal_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the unnormalized normal vector of an edge-vertex pair. + + Parameters + ---------- + v: The vertex position. + e0: The start position of the edge. + e1: The end position of the edge. + + Returns + ------- + The Jacobian of the unnormalized normal vector. + )ipc_qu8mg5v7", + py::arg("v"), py::arg("e0"), py::arg("e1")); + + m.def( + "triangle_unnormalized_normal", &triangle_unnormalized_normal, + R"ipc_qu8mg5v7( + Computes the unnormalized normal vector of a triangle. + + Parameters + ---------- + a: The first vertex of the triangle. + b: The second vertex of the triangle. + c: The third vertex of the triangle. + + Returns + ------- + The unnormalized normal vector of the triangle. + )ipc_qu8mg5v7", + py::arg("a"), py::arg("b"), py::arg("c")); + + m.def( + "triangle_normal", &triangle_normal, + R"ipc_qu8mg5v7( + Computes the normal vector of a triangle. + + Parameters + ---------- + a: The first vertex of the triangle. + b: The second vertex of the triangle. + c: The third vertex of the triangle. + + Returns + ------- + The normal vector of the triangle. + )ipc_qu8mg5v7", + py::arg("a"), py::arg("b"), py::arg("c")); + + m.def( + "triangle_unnormalized_normal_jacobian", + &triangle_unnormalized_normal_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the unnormalized normal vector of a triangle. + + Parameters + ---------- + a: The first vertex of the triangle. + b: The second vertex of the triangle. + c: The third vertex of the triangle. + + Returns + ------- + The Jacobian of the unnormalized normal vector of the triangle. + )ipc_qu8mg5v7", + py::arg("a"), py::arg("b"), py::arg("c")); + + m.def( + "triangle_unnormalized_normal_hessian", + &triangle_unnormalized_normal_hessian, + R"ipc_qu8mg5v7( + Computes the Hessian of the unnormalized normal vector of a triangle. + + Parameters + ---------- + a: The first vertex of the triangle. + b: The second vertex of the triangle. + c: The third vertex of the triangle. + Returns + + ------- + The Hessian of the unnormalized normal vector of the triangle. + )ipc_qu8mg5v7", + py::arg("a"), py::arg("b"), py::arg("c")); + + m.def( + "triangle_normal_jacobian", &triangle_normal_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the normal vector of a triangle. + + Parameters + ---------- + a: The first vertex of the triangle. + b: The second vertex of the triangle. + c: The third vertex of the triangle. + Returns + + ------- + The Jacobian of the normal vector of the triangle. + )ipc_qu8mg5v7", + py::arg("a"), py::arg("b"), py::arg("c")); + + m.def( + "edge_edge_unnormalized_normal", &edge_edge_unnormalized_normal, + R"ipc_qu8mg5v7( + Computes the unnormalized normal vector of two edges. + + Parameters + ---------- + ea0: The first vertex of the first edge. + ea1: The second vertex of the first edge. + eb0: The first vertex of the second edge. + eb1: The second vertex of the second edge. + + Returns + ------- + The unnormalized normal vector of the two edges. + )ipc_qu8mg5v7", + py::arg("ea0"), py::arg("ea1"), py::arg("eb0"), py::arg("eb1")); + + m.def( + "edge_edge_normal", &edge_edge_normal, + R"ipc_qu8mg5v7( + Computes the normal vector of two edges. + + Parameters + ---------- + ea0: The first vertex of the first edge. + ea1: The second vertex of the first edge. + eb0: The first vertex of the second edge. + eb1: The second vertex of the second edge. + + Returns + ------- + The normal vector of the two edges. + )ipc_qu8mg5v7", + py::arg("ea0"), py::arg("ea1"), py::arg("eb0"), py::arg("eb1")); + + m.def( + "edge_edge_unnormalized_normal_jacobian", + &edge_edge_unnormalized_normal_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the unnormalized normal vector of two edges. + + Parameters + ---------- + ea0: The first vertex of the first edge. + ea1: The second vertex of the first edge. + eb0: The first vertex of the second edge. + eb1: The second vertex of the second edge. + + Returns + ------- + The Jacobian of the unnormalized normal vector of the two edges. + )ipc_qu8mg5v7", + py::arg("ea0"), py::arg("ea1"), py::arg("eb0"), py::arg("eb1")); + + m.def( + "edge_edge_normal_jacobian", &edge_edge_normal_jacobian, + R"ipc_qu8mg5v7( + Computes the Jacobian of the normal vector of two edges. + + Parameters + ---------- + ea0: The first vertex of the first edge. + ea1: The second vertex of the first edge. + eb0: The first vertex of the second edge. + eb1: The second vertex of the second edge. + + Returns + ------- + The Jacobian of the normal vector of the two edges. + )ipc_qu8mg5v7", + py::arg("ea0"), py::arg("ea1"), py::arg("eb0"), py::arg("eb1")); +} diff --git a/python/src/math/CMakeLists.txt b/python/src/math/CMakeLists.txt new file mode 100644 index 000000000..6f794c54d --- /dev/null +++ b/python/src/math/CMakeLists.txt @@ -0,0 +1,5 @@ +set(SOURCES + interval.cpp +) + +target_sources(ipctk PRIVATE ${SOURCES}) \ No newline at end of file diff --git a/python/src/math/bindings.hpp b/python/src/math/bindings.hpp new file mode 100644 index 000000000..ac30a52bd --- /dev/null +++ b/python/src/math/bindings.hpp @@ -0,0 +1,5 @@ +#pragma once + +#include + +void define_interval(py::module_& m); \ No newline at end of file diff --git a/python/src/utils/interval.cpp b/python/src/math/interval.cpp similarity index 99% rename from python/src/utils/interval.cpp rename to python/src/math/interval.cpp index c65f47476..ea3358896 100644 --- a/python/src/utils/interval.cpp +++ b/python/src/math/interval.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #ifdef IPC_TOOLKIT_WITH_FILIB using namespace filib; diff --git a/python/src/utils/CMakeLists.txt b/python/src/utils/CMakeLists.txt index 9032db4aa..3de7d7512 100644 --- a/python/src/utils/CMakeLists.txt +++ b/python/src/utils/CMakeLists.txt @@ -1,8 +1,5 @@ set(SOURCES - area_gradient.cpp eigen_ext.cpp - intersection.cpp - interval.cpp logger.cpp thread_limiter.cpp vertex_to_min_edge.cpp diff --git a/python/src/utils/bindings.hpp b/python/src/utils/bindings.hpp index 93794ebc4..9e3b26741 100644 --- a/python/src/utils/bindings.hpp +++ b/python/src/utils/bindings.hpp @@ -2,10 +2,7 @@ #include -void define_area_gradient(py::module_& m); void define_eigen_ext(py::module_& m); -void define_interval(py::module_& m); -void define_intersection(py::module_& m); void define_logger(py::module_& m); void define_thread_limiter(py::module_& m); void define_vertex_to_min_edge(py::module_& m); diff --git a/src/ipc/CMakeLists.txt b/src/ipc/CMakeLists.txt index 75397d903..85d6500d3 100644 --- a/src/ipc/CMakeLists.txt +++ b/src/ipc/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(distance) add_subdirectory(friction) add_subdirectory(geometry) add_subdirectory(implicits) +add_subdirectory(math) add_subdirectory(potentials) add_subdirectory(smooth_contact) add_subdirectory(tangent) diff --git a/src/ipc/ccd/nonlinear_ccd.hpp b/src/ipc/ccd/nonlinear_ccd.hpp index abd68ad9f..c0dcd0594 100644 --- a/src/ipc/ccd/nonlinear_ccd.hpp +++ b/src/ipc/ccd/nonlinear_ccd.hpp @@ -3,7 +3,7 @@ #include #include #ifdef IPC_TOOLKIT_WITH_FILIB -#include +#include #endif #include diff --git a/src/ipc/collision_mesh.cpp b/src/ipc/collision_mesh.cpp index 66aeac933..daef87f67 100644 --- a/src/ipc/collision_mesh.cpp +++ b/src/ipc/collision_mesh.cpp @@ -1,6 +1,6 @@ #include "collision_mesh.hpp" -#include +#include #include #include #include @@ -286,9 +286,7 @@ void CollisionMesh::init_areas() Eigen::VectorXd vertex_edge_areas = Eigen::VectorXd::Constant(num_vertices(), -1); for (int i = 0; i < m_edges.rows(); i++) { - const VectorMax3d e0 = m_rest_positions.row(m_edges(i, 0)); - const VectorMax3d e1 = m_rest_positions.row(m_edges(i, 1)); - double edge_len = (e1 - e0).norm(); + double edge_len = edge_length(i); for (int j = 0; j < m_edges.cols(); j++) { vertex_edge_areas[m_edges(i, j)] = @@ -303,10 +301,10 @@ void CollisionMesh::init_areas() m_edge_areas.setConstant(m_edges.rows(), -1); if (dim() == 3) { for (int i = 0; i < m_faces.rows(); i++) { - const Eigen::Vector3d f0 = m_rest_positions.row(m_faces(i, 0)); - const Eigen::Vector3d f1 = m_rest_positions.row(m_faces(i, 1)); - const Eigen::Vector3d f2 = m_rest_positions.row(m_faces(i, 2)); - double face_area = 0.5 * (f1 - f0).cross(f2 - f0).norm(); + double face_area = triangle_area( + m_rest_positions.row(m_faces(i, 0)), + m_rest_positions.row(m_faces(i, 1)), + m_rest_positions.row(m_faces(i, 2))); for (int j = 0; j < m_faces.cols(); ++j) { vertex_face_areas[m_faces(i, j)] = @@ -330,9 +328,7 @@ void CollisionMesh::init_areas() for (int i = 0; i < m_edge_areas.size(); i++) { if (m_edge_areas[i] < 0) { // Use the edge length for codim edges - const VectorMax3d e0 = m_rest_positions.row(m_edges(i, 0)); - const VectorMax3d e1 = m_rest_positions.row(m_edges(i, 1)); - m_edge_areas[i] = (e1 - e0).norm(); + m_edge_areas[i] = edge_length(i); } } } @@ -508,9 +504,9 @@ Eigen::MatrixXi CollisionMesh::construct_faces_to_edges( double CollisionMesh::edge_length(const index_t edge_id) const { - return (m_rest_positions.row(m_edges(edge_id, 0)) - - m_rest_positions.row(m_edges(edge_id, 1))) - .norm(); + return ::ipc::edge_length( + m_rest_positions.row(m_edges(edge_id, 0)), + m_rest_positions.row(m_edges(edge_id, 1))); } double CollisionMesh::max_edge_length() const diff --git a/src/ipc/geometry/CMakeLists.txt b/src/ipc/geometry/CMakeLists.txt index b3fccc83d..a771e9d6e 100644 --- a/src/ipc/geometry/CMakeLists.txt +++ b/src/ipc/geometry/CMakeLists.txt @@ -1,4 +1,10 @@ set(SOURCES + angle.cpp + angle.hpp + area.cpp + area.hpp + intersection.cpp + intersection.hpp normal.cpp normal.hpp ) diff --git a/src/ipc/geometry/angle.cpp b/src/ipc/geometry/angle.cpp new file mode 100644 index 000000000..935603c69 --- /dev/null +++ b/src/ipc/geometry/angle.cpp @@ -0,0 +1,56 @@ +#include "angle.hpp" + +#include + +namespace ipc { + +double dihedral_angle( + Eigen::ConstRef x0, + Eigen::ConstRef x1, + Eigen::ConstRef x2, + Eigen::ConstRef x3) +{ + double cos_theta = + triangle_normal(x0, x1, x2).dot(triangle_normal(x1, x0, x3)); + assert(std::abs(cos_theta) <= 1.0 + 1e-6); // Numerical tolerance + return std::acos(std::clamp(cos_theta, -1.0, 1.0)); +} + +Eigen::Vector dihedral_angle_gradient( + Eigen::ConstRef x0, + Eigen::ConstRef x1, + Eigen::ConstRef x2, + Eigen::ConstRef x3) +{ + const Eigen::Vector3d n0 = triangle_normal(x0, x1, x2); + const Eigen::Vector3d n1 = triangle_normal(x1, x0, x3); + + const Eigen::Matrix dn0_dx = + triangle_normal_jacobian(x0, x1, x2); + const Eigen::Matrix dn1_dx = + triangle_normal_jacobian(x1, x0, x3); + + const Eigen::Matrix3d dn0_dx0 = dn0_dx.block<3, 3>(0, 0); + const Eigen::Matrix3d dn0_dx1 = dn0_dx.block<3, 3>(0, 3); + const Eigen::Matrix3d dn0_dx2 = dn0_dx.block<3, 3>(0, 6); + const Eigen::Matrix3d dn1_dx0 = dn1_dx.block<3, 3>(0, 3); + const Eigen::Matrix3d dn1_dx1 = dn1_dx.block<3, 3>(0, 0); + const Eigen::Matrix3d dn1_dx3 = dn1_dx.block<3, 3>(0, 6); + + Eigen::Vector gradient; + // Product rule + gradient.segment<3>(0) = + dn0_dx0.transpose() * n1 + dn1_dx0.transpose() * n0; + gradient.segment<3>(3) = + dn0_dx1.transpose() * n1 + dn1_dx1.transpose() * n0; + gradient.segment<3>(6) = dn0_dx2.transpose() * n1; + gradient.segment<3>(9) = dn1_dx3.transpose() * n0; + + // Chain rule through acos + const double cos_theta = n0.dot(n1); + gradient *= -1 / sqrt(1.0 - cos_theta * cos_theta); + + return gradient; +} + +} // namespace ipc \ No newline at end of file diff --git a/src/ipc/geometry/angle.hpp b/src/ipc/geometry/angle.hpp new file mode 100644 index 000000000..f6f435722 --- /dev/null +++ b/src/ipc/geometry/angle.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include + +namespace ipc { + +/// @brief Compute the bending angle between two triangles sharing an edge. +/// x0---x2 +/// | \ | +/// x1---x3 +/// @param x0 The first vertex of the edge. +/// @param x1 The second vertex of the edge. +/// @param x2 The opposite vertex of the first triangle. +/// @param x3 The opposite vertex of the second triangle. +/// @return The bending angle between the two triangles. +double dihedral_angle( + Eigen::ConstRef x0, + Eigen::ConstRef x1, + Eigen::ConstRef x2, + Eigen::ConstRef x3); + +/// @brief Compute the Jacobian of the bending angle between two triangles sharing an edge. +/// x0---x2 +/// | \ | +/// x1---x3 +/// @param x0 The first vertex of the edge. +/// @param x1 The second vertex of the edge. +/// @param x2 The opposite vertex of the first triangle. +/// @param x3 The opposite vertex of the second triangle. +/// @return The Jacobian matrix of the bending angle with respect to the input vertices. +Eigen::Vector dihedral_angle_gradient( + Eigen::ConstRef x0, + Eigen::ConstRef x1, + Eigen::ConstRef x2, + Eigen::ConstRef x3); + +} // namespace ipc \ No newline at end of file diff --git a/src/ipc/utils/area_gradient.cpp b/src/ipc/geometry/area.cpp similarity index 98% rename from src/ipc/utils/area_gradient.cpp rename to src/ipc/geometry/area.cpp index c8f28f9b0..cbb6e223e 100644 --- a/src/ipc/utils/area_gradient.cpp +++ b/src/ipc/geometry/area.cpp @@ -1,4 +1,4 @@ -#include "area_gradient.hpp" +#include "area.hpp" #include diff --git a/src/ipc/utils/area_gradient.hpp b/src/ipc/geometry/area.hpp similarity index 61% rename from src/ipc/utils/area_gradient.hpp rename to src/ipc/geometry/area.hpp index 8cd4f4417..6700a7e48 100644 --- a/src/ipc/utils/area_gradient.hpp +++ b/src/ipc/geometry/area.hpp @@ -4,6 +4,16 @@ namespace ipc { +/// @brief Compute the length of an edge. +/// @param e0 The first vertex of the edge. +/// @param e1 The second vertex of the edge. +/// @return The length of the edge. +inline double +edge_length(Eigen::ConstRef e0, Eigen::ConstRef e1) +{ + return (e1 - e0).norm(); +} + /// @brief Compute the gradient of an edge's length. /// @param e0 The first vertex of the edge. /// @param e1 The second vertex of the edge. @@ -11,6 +21,19 @@ namespace ipc { VectorMax6d edge_length_gradient( Eigen::ConstRef e0, Eigen::ConstRef e1); +/// @brief Compute the area of a triangle. +/// @param t0 The first vertex of the triangle. +/// @param t1 The second vertex of the triangle. +/// @param t2 The third vertex of the triangle. +/// @return The area of the triangle. +inline double triangle_area( + Eigen::ConstRef t0, + Eigen::ConstRef t1, + Eigen::ConstRef t2) +{ + return 0.5 * (t1 - t0).cross(t2 - t0).norm(); +} + /// @brief Compute the gradient of the area of a triangle. /// @param t0 The first vertex of the triangle. /// @param t1 The second vertex of the triangle. diff --git a/src/ipc/utils/intersection.cpp b/src/ipc/geometry/intersection.cpp similarity index 100% rename from src/ipc/utils/intersection.cpp rename to src/ipc/geometry/intersection.cpp diff --git a/src/ipc/utils/intersection.hpp b/src/ipc/geometry/intersection.hpp similarity index 100% rename from src/ipc/utils/intersection.hpp rename to src/ipc/geometry/intersection.hpp diff --git a/src/ipc/geometry/normal.cpp b/src/ipc/geometry/normal.cpp index 088e0788e..877c6e716 100644 --- a/src/ipc/geometry/normal.cpp +++ b/src/ipc/geometry/normal.cpp @@ -87,4 +87,44 @@ MatrixMax edge_vertex_unnormalized_normal_jacobian( return dn; } +// --- triangle normal functions ---------------------------------------------- + +namespace { + void set_cross_product_matrix_jacobian( + Eigen::Ref> Jx, double chain_rule = 1.0) + { + Jx(2, 1) = Jx(3, 2) = Jx(7, 0) = -chain_rule; + Jx(1, 2) = Jx(5, 0) = Jx(6, 1) = chain_rule; + } +} // namespace + +Eigen::Matrix cross_product_matrix_jacobian() +{ + Eigen::Matrix J = Eigen::Matrix::Zero(); + set_cross_product_matrix_jacobian(J); + return J.reshaped(3, 9); +} + +Eigen::Matrix triangle_unnormalized_normal_hessian( + Eigen::ConstRef a, + Eigen::ConstRef b, + Eigen::ConstRef c) +{ + Eigen::Matrix H = Eigen::Matrix::Zero(); + + // ∂²n/∂a² = 0 + set_cross_product_matrix_jacobian(H.block<9, 3>(0, 3), -1.0); // ∂²n/∂a∂b + set_cross_product_matrix_jacobian(H.block<9, 3>(0, 6), 1.0); // ∂²n/∂a∂c + /**/ + set_cross_product_matrix_jacobian(H.block<9, 3>(9, 0), 1.0); // ∂²n/∂b∂a + // ∂²n/∂b² = 0 + set_cross_product_matrix_jacobian(H.block<9, 3>(9, 6), -1.0); // ∂²n/∂b∂c + /**/ + set_cross_product_matrix_jacobian(H.block<9, 3>(18, 0), -1.0); // ∂²n/∂c∂a + set_cross_product_matrix_jacobian(H.block<9, 3>(18, 3), 1.0); // ∂²n/∂c∂b + // ∂²n/∂c² = 0 + + return H.reshaped(3, 81); +} + } // namespace ipc \ No newline at end of file diff --git a/src/ipc/geometry/normal.hpp b/src/ipc/geometry/normal.hpp index d4867e5f3..782bc375a 100644 --- a/src/ipc/geometry/normal.hpp +++ b/src/ipc/geometry/normal.hpp @@ -47,6 +47,10 @@ inline Eigen::Matrix3d cross_product_matrix(Eigen::ConstRef v) return m; } +/// @brief Computes the Jacobian of the cross product matrix. +/// @return The Jacobian of the cross product matrix. +Eigen::Matrix cross_product_matrix_jacobian(); + /** * \defgroup geometry Edge-vertex normal * \brief Functions for computing an edge-vertex normal and resp. Jacobians. @@ -129,12 +133,22 @@ inline Eigen::Matrix triangle_unnormalized_normal_jacobian( Eigen::ConstRef c) { Eigen::Matrix J; - J.middleCols<3>(0) = cross_product_matrix(c - b); - J.middleCols<3>(3) = cross_product_matrix(a - c); - J.middleCols<3>(6) = cross_product_matrix(b - a); + J.middleCols<3>(0) = cross_product_matrix(c - b); // ∂n/∂a + J.middleCols<3>(3) = cross_product_matrix(a - c); // ∂n/∂b + J.middleCols<3>(6) = cross_product_matrix(b - a); // ∂n/∂c return J; } +/// @brief Computes the Hessian of the unnormalized normal vector of a triangle. +/// @param a The first vertex of the triangle. +/// @param b The second vertex of the triangle. +/// @param c The third vertex of the triangle. +/// @return The Hessian of the unnormalized normal vector of the triangle. +Eigen::Matrix triangle_unnormalized_normal_hessian( + Eigen::ConstRef a, + Eigen::ConstRef b, + Eigen::ConstRef c); + /// @brief Computes the Jacobian of the normal vector of a triangle. /// @param a The first vertex of the triangle. /// @param b The second vertex of the triangle. diff --git a/src/ipc/ipc.cpp b/src/ipc/ipc.cpp index d149b65b9..6431690ed 100644 --- a/src/ipc/ipc.cpp +++ b/src/ipc/ipc.cpp @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include #ifdef IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/math/CMakeLists.txt b/src/ipc/math/CMakeLists.txt new file mode 100644 index 000000000..f94256b8e --- /dev/null +++ b/src/ipc/math/CMakeLists.txt @@ -0,0 +1,9 @@ +set(SOURCES + interval.cpp + interval.hpp + math.cpp + math.hpp + math.tpp +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) \ No newline at end of file diff --git a/src/ipc/utils/interval.cpp b/src/ipc/math/interval.cpp similarity index 100% rename from src/ipc/utils/interval.cpp rename to src/ipc/math/interval.cpp diff --git a/src/ipc/utils/interval.hpp b/src/ipc/math/interval.hpp similarity index 100% rename from src/ipc/utils/interval.hpp rename to src/ipc/math/interval.hpp diff --git a/src/ipc/utils/math.cpp b/src/ipc/math/math.cpp similarity index 100% rename from src/ipc/utils/math.cpp rename to src/ipc/math/math.cpp diff --git a/src/ipc/utils/math.hpp b/src/ipc/math/math.hpp similarity index 100% rename from src/ipc/utils/math.hpp rename to src/ipc/math/math.hpp diff --git a/src/ipc/utils/math.tpp b/src/ipc/math/math.tpp similarity index 100% rename from src/ipc/utils/math.tpp rename to src/ipc/math/math.tpp diff --git a/src/ipc/smooth_contact/distance/point_edge.hpp b/src/ipc/smooth_contact/distance/point_edge.hpp index 24181bc7e..9749234a1 100644 --- a/src/ipc/smooth_contact/distance/point_edge.hpp +++ b/src/ipc/smooth_contact/distance/point_edge.hpp @@ -2,9 +2,9 @@ #include #include +#include #include #include -#include #include diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt index 1a551aa87..80786a29b 100644 --- a/src/ipc/utils/CMakeLists.txt +++ b/src/ipc/utils/CMakeLists.txt @@ -1,12 +1,6 @@ set(SOURCES - area_gradient.cpp - area_gradient.hpp eigen_ext.hpp eigen_ext.tpp - intersection.cpp - intersection.hpp - interval.cpp - interval.hpp local_to_global.hpp logger.cpp logger.hpp @@ -18,9 +12,6 @@ set(SOURCES vertex_to_min_edge.cpp vertex_to_min_edge.hpp world_bbox_diagonal_length.hpp - math.cpp - math.hpp - math.tpp autodiff_types.hpp maybe_parallel_for.hpp maybe_parallel_for.tpp diff --git a/tests/src/tests/CMakeLists.txt b/tests/src/tests/CMakeLists.txt index 2d49cfc53..5e9aac565 100644 --- a/tests/src/tests/CMakeLists.txt +++ b/tests/src/tests/CMakeLists.txt @@ -30,6 +30,8 @@ add_subdirectory(ccd) add_subdirectory(collisions) add_subdirectory(distance) add_subdirectory(friction) +add_subdirectory(geometry) +add_subdirectory(math) add_subdirectory(potential) add_subdirectory(tangent) add_subdirectory(utils) \ No newline at end of file diff --git a/tests/src/tests/barrier/test_barrier.cpp b/tests/src/tests/barrier/test_barrier.cpp index 3e7a2da9e..1416f647c 100644 --- a/tests/src/tests/barrier/test_barrier.cpp +++ b/tests/src/tests/barrier/test_barrier.cpp @@ -8,7 +8,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/src/tests/candidates/test_normals.cpp b/tests/src/tests/candidates/test_normals.cpp index e83dd7a01..cf032aa34 100644 --- a/tests/src/tests/candidates/test_normals.cpp +++ b/tests/src/tests/candidates/test_normals.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -190,6 +191,44 @@ TEST_CASE("Face-vertex collision normal", "[fv][normal]") } } +TEST_CASE("Triangle normal hessian", "[normal]") +{ + Eigen::Vector3d a(0, 0, 0); + Eigen::Vector3d b(1, 0, 0); + Eigen::Vector3d c(0, 1, 0); + Eigen::VectorXd x(9); + x << a, b, c; + + // Cross product matrix jacobian + Eigen::MatrixXd J_cross = cross_product_matrix_jacobian(); + Eigen::MatrixXd fd_J_cross; + fd::finite_jacobian( + a, + [](const Eigen::Vector3d& a_fd) { return cross_product_matrix(a_fd); }, + fd_J_cross); + CHECK(fd::compare_jacobian(J_cross, fd_J_cross, 1e-6)); + if (!fd::compare_jacobian(J_cross, fd_J_cross, 1e-6)) { + std::cout << "Hessian:\n" << J_cross << std::endl; + std::cout << "FD Hessian:\n" << fd_J_cross << std::endl; + } + + // Check hessian using finite differences + Eigen::MatrixXd hessian = triangle_unnormalized_normal_hessian(a, b, c); + Eigen::MatrixXd fd_hessian; + fd::finite_jacobian( + x, + [](const Eigen::VectorXd& x_fd) -> Eigen::MatrixXd { + return triangle_unnormalized_normal_jacobian( + x_fd.segment<3>(0), x_fd.segment<3>(3), x_fd.segment<3>(6)); + }, + fd_hessian); + CHECK(fd::compare_jacobian(hessian, fd_hessian, 1e-6)); + if (!fd::compare_jacobian(hessian, fd_hessian, 1e-6)) { + std::cout << "Hessian:\n" << hessian << std::endl; + std::cout << "FD Hessian:\n" << fd_hessian << std::endl; + } +} + TEST_CASE("Plane-vertex collision normal", "[pv][normal]") { Eigen::MatrixXd V(1, 3); diff --git a/tests/src/tests/distance/test_point_point.cpp b/tests/src/tests/distance/test_point_point.cpp index 0fa5c34b5..1a0dc11bd 100644 --- a/tests/src/tests/distance/test_point_point.cpp +++ b/tests/src/tests/distance/test_point_point.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include diff --git a/tests/src/tests/geometry/CMakeLists.txt b/tests/src/tests/geometry/CMakeLists.txt new file mode 100644 index 000000000..1c61142dd --- /dev/null +++ b/tests/src/tests/geometry/CMakeLists.txt @@ -0,0 +1,9 @@ +set(SOURCES + test_angle.cpp +) + +target_sources(ipc_toolkit_tests PRIVATE ${SOURCES}) + +################################################################################ +# Subfolders +################################################################################ \ No newline at end of file diff --git a/tests/src/tests/geometry/test_angle.cpp b/tests/src/tests/geometry/test_angle.cpp new file mode 100644 index 000000000..674d07e35 --- /dev/null +++ b/tests/src/tests/geometry/test_angle.cpp @@ -0,0 +1,59 @@ +#include +#include +#include + +#include + +#include +#include + +#include + +using namespace ipc; + +constexpr double deg2rad(double degrees) { return degrees * igl::PI / 180.0; } + +TEST_CASE("Dihedral angle and gradient", "[angle][dihedral]") +{ + Eigen::Vector3d x0(0, -1, 0); + Eigen::Vector3d x1(0, 1, 0); + Eigen::Vector3d x2(0.5, 0, 0); + Eigen::Vector3d x3(-0.5, 0, 0); + + double expected_angle = deg2rad(GENERATE(0, 15, 30, 45, 60, 75, 90)); + + Eigen::AngleAxisd rotation( + (igl::PI - expected_angle) / 2.0, Eigen::Vector3d::UnitY()); + x2 = rotation * x2; + x3 = rotation.inverse() * x3; + + double angle = dihedral_angle(x0, x1, x2, x3); + + if (expected_angle < igl::PI / 2) { + expected_angle = igl::PI - expected_angle; + } + + CHECK(angle == Catch::Approx(expected_angle)); + + // --- Check the gradient against finite differences --- + + Eigen::VectorXd x(12); + x << x0, x1, x2, x3; + + Eigen::VectorXd grad = dihedral_angle_gradient(x0, x1, x2, x3); + Eigen::VectorXd fd_grad(12); + fd::finite_gradient( + x, + [](const Eigen::VectorXd& x_fd) -> double { + return dihedral_angle( + x_fd.segment<3>(0), x_fd.segment<3>(3), x_fd.segment<3>(6), + x_fd.segment<3>(9)); + }, + fd_grad); + + CHECK(fd::compare_gradient(grad, fd_grad)); + if (!fd::compare_gradient(grad, fd_grad)) { + std::cout << " Gradient:\n" << grad.transpose() << std::endl; + std::cout << "FD Gradient:\n" << fd_grad.transpose() << std::endl; + } +} \ No newline at end of file diff --git a/tests/src/tests/math/CMakeLists.txt b/tests/src/tests/math/CMakeLists.txt new file mode 100644 index 000000000..2e01b5778 --- /dev/null +++ b/tests/src/tests/math/CMakeLists.txt @@ -0,0 +1,14 @@ +set(SOURCES + # Tests + test_interval.cpp + + # Benchmarks + + # Utilities +) + +target_sources(ipc_toolkit_tests PRIVATE ${SOURCES}) + +################################################################################ +# Subfolders +################################################################################ \ No newline at end of file diff --git a/tests/src/tests/utils/test_interval.cpp b/tests/src/tests/math/test_interval.cpp similarity index 99% rename from tests/src/tests/utils/test_interval.cpp rename to tests/src/tests/math/test_interval.cpp index 3c66232fa..22f6d4f82 100644 --- a/tests/src/tests/utils/test_interval.cpp +++ b/tests/src/tests/math/test_interval.cpp @@ -5,7 +5,7 @@ #ifdef IPC_TOOLKIT_WITH_FILIB -#include +#include #include #include diff --git a/tests/src/tests/utils/CMakeLists.txt b/tests/src/tests/utils/CMakeLists.txt index b15998743..2a547cbd7 100644 --- a/tests/src/tests/utils/CMakeLists.txt +++ b/tests/src/tests/utils/CMakeLists.txt @@ -1,6 +1,5 @@ set(SOURCES # Tests - test_interval.cpp test_utils.cpp test_matrixcache.cpp From dc650e52a4e5a19aa507b06d0b13152c8531cf45 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 29 Oct 2025 13:19:58 -0400 Subject: [PATCH 2/2] Fix formatting in documentation for normal vector functions --- python/src/geometry/normal.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/src/geometry/normal.cpp b/python/src/geometry/normal.cpp index 52bfe2885..8984fbbe8 100644 --- a/python/src/geometry/normal.cpp +++ b/python/src/geometry/normal.cpp @@ -122,7 +122,6 @@ void define_normal(py::module_& m) Returns ------- The normal vector. - The normal vector. )ipc_qu8mg5v7", py::arg("v"), py::arg("e0"), py::arg("e1")); @@ -207,8 +206,8 @@ void define_normal(py::module_& m) a: The first vertex of the triangle. b: The second vertex of the triangle. c: The third vertex of the triangle. - Returns + Returns ------- The Hessian of the unnormalized normal vector of the triangle. )ipc_qu8mg5v7", @@ -224,8 +223,8 @@ void define_normal(py::module_& m) a: The first vertex of the triangle. b: The second vertex of the triangle. c: The third vertex of the triangle. - Returns + Returns ------- The Jacobian of the normal vector of the triangle. )ipc_qu8mg5v7",