diff --git a/python/src/barrier/adaptive_stiffness.cpp b/python/src/barrier/adaptive_stiffness.cpp index c308afae1..0f5799df6 100644 --- a/python/src/barrier/adaptive_stiffness.cpp +++ b/python/src/barrier/adaptive_stiffness.cpp @@ -95,7 +95,7 @@ void define_adaptive_stiffness(py::module_& m) static_cast, const NormalCollisions&, Eigen::ConstRef, - const Eigen::SparseMatrix&, const double)>( + const Eigen::SparseMatrix&, const std::optional)>( &semi_implicit_stiffness), R"ipc_Qu8mg5v7( Compute the semi-implicit stiffness's for all collisions. @@ -108,20 +108,20 @@ void define_adaptive_stiffness(py::module_& m) collisions: Normal collisions. vertex_masses: Lumped vertex masses. hess: Hessian of the elasticity energy function. - dmin: Minimum distance between elements. + dmin: Deprecated minimum distance between elements. If given, it overrides each stencil's own (a collision's dmin, or a candidate's thickness). Returns: The semi-implicit stiffness's. )ipc_Qu8mg5v7", "mesh"_a, "vertices"_a, "collisions"_a, "vertex_masses"_a, "hess"_a, - "dmin"_a); + "dmin"_a = py::none()); m.def( "semi_implicit_stiffness", static_cast, const Candidates&, Eigen::ConstRef, - const Eigen::SparseMatrix&, const double)>( + const Eigen::SparseMatrix&, const std::optional)>( &semi_implicit_stiffness), R"ipc_Qu8mg5v7( Compute the semi-implicit stiffness's for all collisions. @@ -134,11 +134,11 @@ void define_adaptive_stiffness(py::module_& m) collisions: Collisions candidates. vertex_masses: Lumped vertex masses. hess: Hessian of the elasticity energy function. - dmin: Minimum distance between elements. + dmin: Deprecated minimum distance between elements. If given, it overrides each stencil's own (a collision's dmin, or a candidate's thickness). Returns: The semi-implicit stiffness's. )ipc_Qu8mg5v7", "mesh"_a, "vertices"_a, "collisions"_a, "vertex_masses"_a, "hess"_a, - "dmin"_a); + "dmin"_a = py::none()); } diff --git a/python/src/candidates/candidates.cpp b/python/src/candidates/candidates.cpp index d902a294b..40a61ab9a 100644 --- a/python/src/candidates/candidates.cpp +++ b/python/src/candidates/candidates.cpp @@ -76,13 +76,14 @@ void define_candidates(py::module_& m) mesh: The collision mesh. vertices_t0: Surface vertex starting positions (rowwise). vertices_t1: Surface vertex ending positions (rowwise). - min_distance: The minimum distance allowable between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). narrow_phase_ccd: The narrow phase CCD algorithm to use. Returns: True if any collisions occur. )ipc_Qu8mg5v7", - "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, "min_distance"_a = 0.0, + "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, + "min_distance"_a = py::none(), "narrow_phase_ccd"_a = DEFAULT_NARROW_PHASE_CCD) .def( "compute_collision_free_stepsize", @@ -97,13 +98,14 @@ void define_candidates(py::module_& m) mesh: The collision mesh. vertices_t0: Surface vertex starting positions (rowwise). Assumed to be intersection free. vertices_t1: Surface vertex ending positions (rowwise). - min_distance: The minimum distance allowable between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). narrow_phase_ccd: The narrow phase CCD algorithm to use. Returns: A step-size :math:`\in [0, 1]` that is collision free. A value of 1.0 if a full step and 0.0 is no step. )ipc_Qu8mg5v7", - "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, "min_distance"_a = 0.0, + "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, + "min_distance"_a = py::none(), "narrow_phase_ccd"_a = DEFAULT_NARROW_PHASE_CCD) .def( "compute_noncandidate_conservative_stepsize", @@ -130,7 +132,7 @@ void define_candidates(py::module_& m) vertices_t0: Surface vertex starting positions (rowwise). vertices_t1: Surface vertex ending positions (rowwise). dhat: Barrier activation distance. - min_distance: Minimum distance allowable between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). broad_phase: Broad phase algorithm to use. narrow_phase_ccd: Narrow phase CCD algorithm to use. @@ -138,7 +140,7 @@ void define_candidates(py::module_& m) A step-size :math:`\in [0, 1]` that is collision free. )ipc_Qu8mg5v7", "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, "dhat"_a, - "min_distance"_a = 0.0, "broad_phase"_a = nullptr, + "min_distance"_a = py::none(), "broad_phase"_a = nullptr, "narrow_phase_ccd"_a = DEFAULT_NARROW_PHASE_CCD) .def( "compute_per_vertex_safe_distances", @@ -153,13 +155,13 @@ void define_candidates(py::module_& m) mesh: The collision mesh. vertices: Collision mesh vertex positions (rowwise). inflation_radius: The inflation radius used to build the candidates. - min_distance: The minimum allowable distance between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). Returns: A vector of minimum distances, one for each vertex. )ipc_Qu8mg5v7", "mesh"_a, "vertices"_a, "inflation_radius"_a, - "min_distance"_a = 0.0) + "min_distance"_a = py::none()) .def( "edge_vertex_to_vertex_vertex", [](const Candidates& C, const CollisionMesh& mesh, diff --git a/python/src/candidates/stencil_methods.hpp b/python/src/candidates/stencil_methods.hpp index 25214b935..29bb98a6e 100644 --- a/python/src/candidates/stencil_methods.hpp +++ b/python/src/candidates/stencil_methods.hpp @@ -3,8 +3,14 @@ #include #include +#include +#include +#include +#include +#include #include +#include #include namespace ipc { @@ -55,6 +61,45 @@ template void define_stencil_methods(PyClass& cls) }, "Select this stencil's DOF from the full matrix of DOF.", "X"_a, "edges"_a, "faces"_a) + .def( + "thickness", + [](const C& self, Eigen::ConstRef vertex_thickness, + Eigen::ConstRef edge_thickness, + Eigen::ConstRef face_thickness) { + // C++ indexes the vectors unchecked, so an empty or short one + // (e.g., a mesh without a thickness) must not reach it. + const auto check = [](Eigen::ConstRef t, + const index_t i, const char* name) { + if (i < 0 || i >= t.size()) { + throw py::index_error( + fmt::format( + "{} has no entry {} (size {})", name, i, + t.size())); + } + }; + if constexpr (std::is_same_v) { + check( + vertex_thickness, self.vertex0_id, "vertex_thickness"); + check( + vertex_thickness, self.vertex1_id, "vertex_thickness"); + } else if constexpr (std::is_same_v) { + check(edge_thickness, self.edge_id, "edge_thickness"); + check(vertex_thickness, self.vertex_id, "vertex_thickness"); + } else if constexpr (std::is_same_v) { + check(edge_thickness, self.edge0_id, "edge_thickness"); + check(edge_thickness, self.edge1_id, "edge_thickness"); + } else if constexpr (std::is_same_v) { + check(face_thickness, self.face_id, "face_thickness"); + check(vertex_thickness, self.vertex_id, "vertex_thickness"); + } else { + static_assert(std::is_same_v); + check(vertex_thickness, self.vertex_id, "vertex_thickness"); + } + return self.thickness( + vertex_thickness, edge_thickness, face_thickness); + }, + "Get the sum of the two primitives' thicknesses (see CollisionMesh.set_thickness).", + "vertex_thickness"_a, "edge_thickness"_a, "face_thickness"_a) // -- distance ------------------------------------------------------ .def( "compute_distance", diff --git a/python/src/collision_mesh.cpp b/python/src/collision_mesh.cpp index ab2dc7e93..327e48a14 100644 --- a/python/src/collision_mesh.cpp +++ b/python/src/collision_mesh.cpp @@ -519,6 +519,69 @@ void define_collision_mesh(py::module_& m) "are_area_jacobians_initialized", &CollisionMesh::are_area_jacobians_initialized, "Determine if the area Jacobians have been initialized by calling init_area_jacobians().") + .def( + "has_thickness", &CollisionMesh::has_thickness, + "Determine if the primitives have a thickness (see set_thickness).") + .def( + "vertex_thickness", &CollisionMesh::vertex_thickness, + "Get the thickness of a vertex (zero without a thickness).", "vi"_a) + .def( + "edge_thickness", &CollisionMesh::edge_thickness, + "Get the thickness of an edge (zero without a thickness).", "ei"_a) + .def( + "face_thickness", &CollisionMesh::face_thickness, + "Get the thickness of a face (zero without a thickness).", "fi"_a) + .def_property_readonly( + "vertex_thicknesses", &CollisionMesh::vertex_thicknesses, + "The thickness of the vertices (empty without a thickness).") + .def_property_readonly( + "edge_thicknesses", &CollisionMesh::edge_thicknesses, + "The thickness of the edges (empty without a thickness).") + .def_property_readonly( + "face_thicknesses", &CollisionMesh::face_thicknesses, + "The thickness of the faces (empty without a thickness).") + .def_property_readonly( + "max_thickness", &CollisionMesh::max_thickness, + "The largest thickness of any primitive (zero without a thickness).") + .def( + "set_thickness", &CollisionMesh::set_thickness, + R"ipc_Qu8mg5v7( + Give every vertex, edge, and face a thickness. + + Each primitive is offset outward by its thickness, so two primitives touch when they are the sum of their thicknesses apart. Collisions built from this mesh add that sum to their minimum distance (NormalCollision.dmin), the broad phase inflates each primitive's box by its thickness, and CCD stops each pair at that sum plus the minimum distance. A sheet of total thickness h, for example, gives each of its primitives h / 2. + + Parameters: + vertex_thickness: Thickness of each vertex (size num_vertices). + edge_thickness: Thickness of each edge (size num_edges). + face_thickness: Thickness of each face (size num_faces). + )ipc_Qu8mg5v7", + "vertex_thickness"_a, "edge_thickness"_a, "face_thickness"_a) + .def( + "set_vertex_thickness", &CollisionMesh::set_vertex_thickness, + R"ipc_Qu8mg5v7( + Give each vertex a thickness, which its edges and faces inherit. + + An edge or face takes the largest thickness of its vertices. + + Parameters: + vertex_thickness: Thickness of each vertex (size num_vertices). + )ipc_Qu8mg5v7", + "vertex_thickness"_a) + .def( + "set_face_thickness", &CollisionMesh::set_face_thickness, + R"ipc_Qu8mg5v7( + Give each face a thickness, which its edges and vertices inherit. + + An edge or vertex takes the largest thickness of the faces it bounds. An edge on no face takes the larger thickness of its two vertices. + + Parameters: + face_thickness: Thickness of each face (size num_faces). + vertex_thickness: Least thickness of each vertex (size num_vertices), or empty for zero. It is the whole thickness of a vertex on no face. + )ipc_Qu8mg5v7", + "face_thickness"_a, "vertex_thickness"_a = Eigen::VectorXd()) + .def( + "clear_thickness", &CollisionMesh::clear_thickness, + "Remove the primitives' thickness.") .def_static( "construct_is_on_surface", &CollisionMesh::construct_is_on_surface, R"ipc_Qu8mg5v7( diff --git a/python/src/collisions/normal/normal_collisions.cpp b/python/src/collisions/normal/normal_collisions.cpp index d590d2b9e..cdc6df4fb 100644 --- a/python/src/collisions/normal/normal_collisions.cpp +++ b/python/src/collisions/normal/normal_collisions.cpp @@ -166,7 +166,7 @@ void define_normal_collisions(py::module_& m) "build", py::overload_cast< const CollisionMesh&, Eigen::ConstRef, - const double, const double, BroadPhase*>( + const double, const std::optional, BroadPhase*>( &NormalCollisions::build), R"ipc_Qu8mg5v7( Initialize the set of collisions used to compute the barrier potential. @@ -175,17 +175,17 @@ void define_normal_collisions(py::module_& m) mesh: The collision mesh. vertices: Vertices of the collision mesh. dhat: The activation distance of the barrier. - dmin: Minimum distance. + dmin: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every collision keeps this minimum distance and the mesh's thickness is ignored; otherwise each keeps its two primitives' thicknesses. broad_phase: Broad-phase to use. )ipc_Qu8mg5v7", - "mesh"_a, "vertices"_a, "dhat"_a, "dmin"_a = 0, + "mesh"_a, "vertices"_a, "dhat"_a, "dmin"_a = py::none(), "broad_phase"_a = nullptr) .def( "build", py::overload_cast< const Candidates&, const CollisionMesh&, - Eigen::ConstRef, const double, const double>( - &NormalCollisions::build), + Eigen::ConstRef, const double, + const std::optional>(&NormalCollisions::build), R"ipc_Qu8mg5v7( Initialize the set of collisions used to compute the barrier potential. @@ -194,12 +194,13 @@ void define_normal_collisions(py::module_& m) mesh: The collision mesh. vertices: Vertices of the collision mesh. dhat: The activation distance of the barrier. - dmin: Minimum distance. + dmin: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every collision keeps this minimum distance and the mesh's thickness is ignored; otherwise each keeps its two primitives' thicknesses. Raises: ValueError: In 3D, if the candidates were built with all_types=True. )ipc_Qu8mg5v7", - "candidates"_a, "mesh"_a, "vertices"_a, "dhat"_a, "dmin"_a = 0) + "candidates"_a, "mesh"_a, "vertices"_a, "dhat"_a, + "dmin"_a = py::none()) .def( "compute_minimum_distance", &NormalCollisions::compute_minimum_distance, diff --git a/python/src/collisions/tangential/tangential_collision.cpp b/python/src/collisions/tangential/tangential_collision.cpp index 2d84fed6e..8baf9c779 100644 --- a/python/src/collisions/tangential/tangential_collision.cpp +++ b/python/src/collisions/tangential/tangential_collision.cpp @@ -121,6 +121,9 @@ void define_tangential_collision(py::module_& m) "normal_force_magnitude", &TangentialCollision::normal_force_magnitude, "Normal force magnitude") + .def_readwrite( + "dmin", &TangentialCollision::dmin, + "Minimum distance of the normal collision this was built from") .def_readwrite( "mu_s", &TangentialCollision::mu_s, "Ratio between normal and static tangential forces (e.g., friction coefficient)") diff --git a/python/src/ipc.cpp b/python/src/ipc.cpp index 9cb5f967e..8d897281c 100644 --- a/python/src/ipc.cpp +++ b/python/src/ipc.cpp @@ -23,15 +23,15 @@ void define_ipc(py::module_& m) mesh: The collision mesh. vertices_t0: Surface vertex vertices at start as rows of a matrix. vertices_t1: Surface vertex vertices at end as rows of a matrix. - min_distance: The minimum distance allowable between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). broad_phase: Broad phase to use. narrow_phase_ccd: The narrow phase CCD algorithm to use. Returns: True if any collisions occur. )ipc_Qu8mg5v7", - "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, "min_distance"_a = 0.0, - "broad_phase"_a = nullptr, + "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, + "min_distance"_a = py::none(), "broad_phase"_a = nullptr, py::arg_v( "narrow_phase_ccd", DEFAULT_NARROW_PHASE_CCD, "TightInclusionCCD()")); @@ -51,15 +51,15 @@ void define_ipc(py::module_& m) mesh: The collision mesh. vertices_t0: Vertex vertices at start as rows of a matrix. Assumes vertices_t0 is intersection free. vertices_t1: Surface vertex vertices at end as rows of a matrix. - min_distance: The minimum distance allowable between any two elements. + min_distance: Deprecated: give the mesh a thickness (CollisionMesh.set_thickness) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). broad_phase: Broad phase to use. narrow_phase_ccd: The narrow phase CCD algorithm to use. Returns: A step-size :math:`\in [0, 1]` that is collision free. A value of 1.0 if a full step and 0.0 is no step. )ipc_Qu8mg5v7", - "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, "min_distance"_a = 0.0, - "broad_phase"_a = nullptr, + "mesh"_a, "vertices_t0"_a, "vertices_t1"_a, + "min_distance"_a = py::none(), "broad_phase"_a = nullptr, py::arg_v( "narrow_phase_ccd", DEFAULT_NARROW_PHASE_CCD, "TightInclusionCCD()")); diff --git a/python/src/ogc/trust_region.cpp b/python/src/ogc/trust_region.cpp index 9b6eb2dbd..772db71f6 100644 --- a/python/src/ogc/trust_region.cpp +++ b/python/src/ogc/trust_region.cpp @@ -34,11 +34,11 @@ void define_trust_region(py::module_& m) pred_x: Predicted vertex positions. collisions: Collisions to be initialized. dhat: The offset distance for contact. - min_distance: Minimum distance between elements. + min_distance: Deprecated minimum distance between elements, which overrides the mesh's thickness when given. broad_phase: Broad phase collision detection. )ipc_Qu8mg5v7", "mesh"_a, "x"_a, "pred_x"_a, "collisions"_a, "dhat"_a, - "min_distance"_a = 0.0, "broad_phase"_a = nullptr) + "min_distance"_a = py::none(), "broad_phase"_a = nullptr) .def( "update", &ogc::TrustRegion::update, R"ipc_Qu8mg5v7( @@ -48,10 +48,10 @@ void define_trust_region(py::module_& m) mesh: The collision mesh. x: Current vertex positions. collisions: Collisions to be updated. - min_distance: Minimum distance between elements. + min_distance: Deprecated minimum distance between elements, which overrides the mesh's thickness when given. broad_phase: Broad phase collision detection. )ipc_Qu8mg5v7", - "mesh"_a, "x"_a, "collisions"_a, "min_distance"_a = 0.0, + "mesh"_a, "x"_a, "collisions"_a, "min_distance"_a = py::none(), "broad_phase"_a = nullptr) .def( "update_if_needed", &ogc::TrustRegion::update_if_needed, @@ -62,10 +62,10 @@ void define_trust_region(py::module_& m) mesh: The collision mesh. x: Current vertex positions. collisions: Collisions to be updated. - min_distance: Minimum distance between elements. + min_distance: Deprecated minimum distance between elements, which overrides the mesh's thickness when given. broad_phase: Broad phase collision detection. )ipc_Qu8mg5v7", - "mesh"_a, "x"_a, "collisions"_a, "min_distance"_a = 0.0, + "mesh"_a, "x"_a, "collisions"_a, "min_distance"_a = py::none(), "broad_phase"_a = nullptr) .def( "filter_step", &ogc::TrustRegion::filter_step, diff --git a/python/src/potentials/tangential_potential.cpp b/python/src/potentials/tangential_potential.cpp index 630f84be6..5b1a4625a 100644 --- a/python/src/potentials/tangential_potential.cpp +++ b/python/src/potentials/tangential_potential.cpp @@ -33,7 +33,7 @@ void define_tangential_potential(py::module_& m) Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, const NormalPotential&, - const double, const bool>( + const std::optional, const bool>( &TangentialPotential::force, py::const_), R"ipc_Qu8mg5v7( Compute the friction force from the given velocities. @@ -45,7 +45,7 @@ void define_tangential_potential(py::module_& m) lagged_displacements: Previous displacements of the vertices (rowwise). velocities: Current displacements of the vertices (rowwise). normal_potential: Normal potential (used for normal force magnitude). - dmin: Minimum distance (used for normal force magnitude). + dmin: Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin. no_mu: whether to not multiply by mu Returns: @@ -53,7 +53,7 @@ void define_tangential_potential(py::module_& m) )ipc_Qu8mg5v7", "collisions"_a, "mesh"_a, "rest_positions"_a, "lagged_displacements"_a, "velocities"_a, "normal_potential"_a, - "dmin"_a = 0, "no_mu"_a = false) + "dmin"_a = py::none(), "no_mu"_a = false) .def( "force_jacobian", py::overload_cast< @@ -61,7 +61,8 @@ void define_tangential_potential(py::module_& m) Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, const NormalPotential&, - const TangentialPotential::DiffWRT, const double>( + const TangentialPotential::DiffWRT, + const std::optional>( &TangentialPotential::force_jacobian, py::const_), R"ipc_Qu8mg5v7( Compute the Jacobian of the friction force wrt the velocities. @@ -74,21 +75,21 @@ void define_tangential_potential(py::module_& m) velocities: Current displacements of the vertices (rowwise). normal_potential: Normal potential (used for normal force magnitude). wrt: The variable to take the derivative with respect to. - dmin: Minimum distance (used for normal force magnitude). + dmin: Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin. Returns: The Jacobian of the friction force wrt the velocities. )ipc_Qu8mg5v7", "collisions"_a, "mesh"_a, "rest_positions"_a, "lagged_displacements"_a, "velocities"_a, "normal_potential"_a, - "wrt"_a, "dmin"_a = 0) + "wrt"_a, "dmin"_a = py::none()) .def( "force", py::overload_cast< const TangentialCollision&, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, - const NormalPotential&, const double, const bool>( - &TangentialPotential::force, py::const_), + const NormalPotential&, const std::optional, + const bool>(&TangentialPotential::force, py::const_), R"ipc_Qu8mg5v7( Compute the friction force. @@ -98,14 +99,14 @@ void define_tangential_potential(py::module_& m) lagged_displacements: Previous displacements of the vertices (rowwise). velocities: Current displacements of the vertices (rowwise). normal_potential: Normal potential (used for normal force magnitude). - dmin: Minimum distance (used for normal force magnitude). + dmin: Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin. no_mu: Whether to not multiply by mu Returns: Friction force )ipc_Qu8mg5v7", "collision"_a, "rest_positions"_a, "lagged_displacements"_a, - "velocities"_a, "normal_potential"_a, "dmin"_a = 0, + "velocities"_a, "normal_potential"_a, "dmin"_a = py::none(), "no_mu"_a = false) .def( "force_jacobian", @@ -113,7 +114,8 @@ void define_tangential_potential(py::module_& m) const TangentialCollision&, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, const NormalPotential&, const TangentialPotential::DiffWRT, - const double>(&TangentialPotential::force_jacobian, py::const_), + const std::optional>( + &TangentialPotential::force_jacobian, py::const_), R"ipc_Qu8mg5v7( Compute the friction force Jacobian. @@ -124,11 +126,12 @@ void define_tangential_potential(py::module_& m) velocities: Current displacements of the vertices (rowwise). normal_potential: Normal potential (used for normal force magnitude). wrt: Variable to differentiate the friction force with respect to. - dmin: Minimum distance (used for normal force magnitude). + dmin: Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin. Returns: Friction force Jacobian )ipc_Qu8mg5v7", "collision"_a, "rest_positions"_a, "lagged_displacements"_a, - "velocities"_a, "normal_potential"_a, "wrt"_a, "dmin"_a = 0); + "velocities"_a, "normal_potential"_a, "wrt"_a, + "dmin"_a = py::none()); } diff --git a/python/tests/test_collision_mesh.py b/python/tests/test_collision_mesh.py index f050ca385..c08d9e734 100644 --- a/python/tests/test_collision_mesh.py +++ b/python/tests/test_collision_mesh.py @@ -5,7 +5,9 @@ import scipy from ipctk import ( CollisionMesh, + EdgeEdgeCandidate, NormalCollisions, + compute_collision_free_stepsize, make_sparse_filter, make_vertex_patches_filter, ) @@ -263,3 +265,55 @@ def test_zero_custom_areas_drop_their_collisions(): # Zero-weight collisions are dropped, so no contact force is produced. assert len(collisions) == 0 + + +def test_thickness(): + V, E = _two_segments() # 0.1 apart + mesh = CollisionMesh(V, E) + assert not mesh.has_thickness() + + collisions = NormalCollisions() + collisions.build(mesh, V, dhat=0.06) + assert len(collisions) == 0 + + # Each segment is 0.025 thick, so the pair rests 0.05 apart and is + # within dhat of contact. + mesh.set_vertex_thickness(np.full(4, 0.025)) + assert mesh.has_thickness() + assert mesh.max_thickness == 0.025 + assert np.allclose(mesh.edge_thicknesses, 0.025) + assert ( + EdgeEdgeCandidate(0, 1).thickness( + mesh.vertex_thicknesses, + mesh.edge_thicknesses, + mesh.face_thicknesses, + ) + == 0.05 + ) + + collisions.build(mesh, V, dhat=0.06) + assert len(collisions) > 0 + for i in range(len(collisions)): + assert np.isclose(collisions[i].dmin, 0.05) + + # Moving the top segment down stops 0.05 above the bottom one. + V1 = V.copy() + V1[2:, 1] = -0.1 + step = compute_collision_free_stepsize(mesh, V, V1) + assert 0 < step <= 0.25 + # A min_distance (deprecated) overrides the thickness. + assert compute_collision_free_stepsize(mesh, V, V1, min_distance=0) > step + + mesh.clear_thickness() + assert not mesh.has_thickness() + # Without a thickness the arrays are empty, which raises instead of + # reading out of bounds. + try: + EdgeEdgeCandidate(0, 1).thickness( + mesh.vertex_thicknesses, + mesh.edge_thicknesses, + mesh.face_thicknesses, + ) + assert False + except IndexError: + pass diff --git a/src/ipc/barrier/adaptive_stiffness.cpp b/src/ipc/barrier/adaptive_stiffness.cpp index 4aa697fb3..7766748b8 100644 --- a/src/ipc/barrier/adaptive_stiffness.cpp +++ b/src/ipc/barrier/adaptive_stiffness.cpp @@ -4,10 +4,12 @@ #include #include +#include #include #include // std::min/max #include +#include namespace ipc { @@ -92,8 +94,12 @@ Eigen::VectorXd semi_implicit_stiffness( const StencilsT& collisions, Eigen::ConstRef vertex_masses, const Eigen::SparseMatrix& hess, - const double dmin) + const std::optional dmin) { + if (dmin.has_value()) { + warn_dmin_deprecated(); + } + const int dim = mesh.dim(); assert(vertices.cols() == dim); // Vertex positions must be 3D assert(hess.rows() == hess.cols()); // Hessian must be square @@ -138,8 +144,19 @@ Eigen::VectorXd semi_implicit_stiffness( } } + // A collision keeps the dmin it was built with; a candidate its + // two primitives' thicknesses. + double stencil_dmin; + if constexpr (std::is_base_of_v< + NormalCollision, + std::decay_t>) { + stencil_dmin = dmin.value_or(collision.dmin); + } else { + stencil_dmin = dmin.value_or(mesh.pair_thickness(collision)); + } + return semi_implicit_stiffness( - collision, positions, local_mass, local_hess, dmin); + collision, positions, local_mass, local_hess, stencil_dmin); }); } @@ -152,7 +169,7 @@ template Eigen::VectorXd semi_implicit_stiffness( const NormalCollisions& collisions, Eigen::ConstRef vertex_masses, const Eigen::SparseMatrix& hess, - const double dmin); + const std::optional dmin); template Eigen::VectorXd semi_implicit_stiffness( const CollisionMesh& mesh, @@ -160,7 +177,7 @@ template Eigen::VectorXd semi_implicit_stiffness( const Candidates& collisions, Eigen::ConstRef vertex_masses, const Eigen::SparseMatrix& hess, - const double dmin); + const std::optional dmin); // ----------------------------------------------------------------------------- diff --git a/src/ipc/barrier/adaptive_stiffness.hpp b/src/ipc/barrier/adaptive_stiffness.hpp index 70a72e4db..8d5893e6d 100644 --- a/src/ipc/barrier/adaptive_stiffness.hpp +++ b/src/ipc/barrier/adaptive_stiffness.hpp @@ -8,6 +8,8 @@ #include +#include + namespace ipc { /// @brief Compute an inital barrier stiffness using the barrier potential gradient. @@ -104,7 +106,7 @@ double semi_implicit_stiffness( /// @param collisions Normal collisions or collision candidates. /// @param vertex_masses Lumped vertex masses. /// @param hess Hessian of the elasticity energy function. -/// @param dmin Minimum distance between elements. +/// @param dmin Deprecated minimum distance between elements. If given, it overrides each stencil's own: a collision's dmin, or a candidate's two primitives' thicknesses (CollisionMesh::set_thickness()). /// @return The semi-implicit stiffness's. template Eigen::VectorXd semi_implicit_stiffness( @@ -113,6 +115,6 @@ Eigen::VectorXd semi_implicit_stiffness( const StencilsT& collisions, Eigen::ConstRef vertex_masses, const Eigen::SparseMatrix& hess, - const double dmin); + const std::optional dmin = std::nullopt); } // namespace ipc diff --git a/src/ipc/broad_phase/aabb.cpp b/src/ipc/broad_phase/aabb.cpp index 2a5edefcd..04dbd5521 100644 --- a/src/ipc/broad_phase/aabb.cpp +++ b/src/ipc/broad_phase/aabb.cpp @@ -84,6 +84,44 @@ void build_vertex_boxes( }); } +void build_vertex_boxes( + Eigen::ConstRef vertices, + AABBs& vertex_boxes, + const double inflation_radius, + Eigen::ConstRef vertex_inflation) +{ + IPC_TOOLKIT_PROFILE_BLOCK("build_vertex_boxes"); + assert(vertex_inflation.size() == vertices.rows()); + + vertex_boxes.resize(vertices.rows()); + + tbb::parallel_for(0, index_t(vertices.rows()), [&](index_t i) { + vertex_boxes[i] = AABB::from_point( + vertices.row(i), inflation_radius + vertex_inflation[i]); + vertex_boxes[i].vertex_ids = { { i, -1, -1 } }; + }); +} + +void build_vertex_boxes( + Eigen::ConstRef vertices_t0, + Eigen::ConstRef vertices_t1, + AABBs& vertex_boxes, + const double inflation_radius, + Eigen::ConstRef vertex_inflation) +{ + IPC_TOOLKIT_PROFILE_BLOCK("build_vertex_boxes"); + assert(vertex_inflation.size() == vertices_t0.rows()); + + vertex_boxes.resize(vertices_t0.rows()); + + tbb::parallel_for(0, index_t(vertices_t0.rows()), [&](index_t i) { + vertex_boxes[i] = AABB::from_point( + vertices_t0.row(i), vertices_t1.row(i), + inflation_radius + vertex_inflation[i]); + vertex_boxes[i].vertex_ids = { { i, -1, -1 } }; + }); +} + void build_edge_boxes( const AABBs& vertex_boxes, Eigen::ConstRef edges, diff --git a/src/ipc/broad_phase/aabb.hpp b/src/ipc/broad_phase/aabb.hpp index bc38bfcd4..0a790c426 100644 --- a/src/ipc/broad_phase/aabb.hpp +++ b/src/ipc/broad_phase/aabb.hpp @@ -121,6 +121,30 @@ void build_vertex_boxes( AABBs& vertex_boxes, const double inflation_radius = 0); +/// @brief Build one AABB per vertex position, each inflated by its own radius. +/// @param[in] vertices Vertex positions (rowwise). +/// @param[out] vertex_boxes Vertex AABBs. +/// @param[in] inflation_radius Radius of a sphere around the points which the AABBs enclose. +/// @param[in] vertex_inflation Additional radius of each vertex's sphere (size vertices.rows()). +void build_vertex_boxes( + Eigen::ConstRef vertices, + AABBs& vertex_boxes, + const double inflation_radius, + Eigen::ConstRef vertex_inflation); + +/// @brief Build one AABB per moving vertex, each inflated by its own radius. +/// @param vertices_t0 Vertex positions at t=0 (rowwise). +/// @param vertices_t1 Vertex positions at t=1 (rowwise). +/// @param vertex_boxes Vertex AABBs. +/// @param inflation_radius Radius of a capsule around the temporal edges which the AABBs enclose. +/// @param vertex_inflation Additional radius of each vertex's capsule (size vertices_t0.rows()). +void build_vertex_boxes( + Eigen::ConstRef vertices_t0, + Eigen::ConstRef vertices_t1, + AABBs& vertex_boxes, + const double inflation_radius, + Eigen::ConstRef vertex_inflation); + /// @brief Build one AABB per edge. /// @param vertex_boxes Vertex AABBs. /// @param edges Edges (rowwise). diff --git a/src/ipc/candidates/CMakeLists.txt b/src/ipc/candidates/CMakeLists.txt index a7155feef..6feb252c7 100644 --- a/src/ipc/candidates/CMakeLists.txt +++ b/src/ipc/candidates/CMakeLists.txt @@ -1,6 +1,7 @@ set(SOURCES candidates.cpp candidates.hpp + min_distance.hpp collision_stencil.hpp stencil_adapter.hpp stencil_mixin.hpp diff --git a/src/ipc/candidates/candidates.cpp b/src/ipc/candidates/candidates.cpp index 5dbeb5211..fc1cdb599 100644 --- a/src/ipc/candidates/candidates.cpp +++ b/src/ipc/candidates/candidates.cpp @@ -2,8 +2,10 @@ #include #include +#include #include #include +#include #include #include #include @@ -44,6 +46,30 @@ namespace { { return E_padded.leftCols(2); } + + /// @brief How much to inflate each vertex's box for the mesh's thickness. + /// + /// The broad phase builds an edge's or face's box as the union of its + /// vertices' boxes, so a vertex's box is inflated by the largest thickness + /// of the vertex and the edges and faces it bounds. The box of every + /// primitive then covers the primitive's own thickness. + Eigen::VectorXd vertex_box_thickness(const CollisionMesh& mesh) + { + Eigen::VectorXd thickness = mesh.vertex_thicknesses(); + for (int i = 0; i < mesh.edges().rows(); i++) { + for (int j = 0; j < 2; j++) { + double& t = thickness[mesh.edges()(i, j)]; + t = std::max(t, mesh.edge_thickness(i)); + } + } + for (int i = 0; i < mesh.faces().rows(); i++) { + for (int j = 0; j < 3; j++) { + double& t = thickness[mesh.faces()(i, j)]; + t = std::max(t, mesh.face_thickness(i)); + } + } + return thickness; + } } // namespace void Candidates::build( @@ -65,8 +91,26 @@ void Candidates::build( clear(); + // Each vertex's box is inflated by its thickness on top of + // inflation_radius, so the boxes of two primitives overlap whenever they + // are within the sum of their thicknesses plus 2 * inflation_radius. + const Eigen::VectorXd box_thickness = + mesh.has_thickness() ? vertex_box_thickness(mesh) : Eigen::VectorXd(); + const auto build_broad_phase = [&](Eigen::ConstRef V, + Eigen::ConstRef E, + Eigen::ConstRef F, + Eigen::ConstRef V_t) { + if (V_t.size() == 0) { + broad_phase->build(V, E, F, inflation_radius); + return; + } + AABBs boxes; + build_vertex_boxes(V, boxes, inflation_radius, V_t); + broad_phase->build(boxes, E, F, uint8_t(dim)); + }; + broad_phase->can_vertices_collide = mesh.can_collide; - broad_phase->build(vertices, mesh.edges(), mesh.faces(), inflation_radius); + build_broad_phase(vertices, mesh.edges(), mesh.faces(), box_thickness); broad_phase->detect_collision_candidates(*this, all_types); has_all_types = all_types; @@ -80,9 +124,11 @@ void Candidates::build( // Codim. vertices to codim. vertices: if (detect_codim && mesh.num_codim_vertices()) { broad_phase->clear(); - broad_phase->build( + build_broad_phase( vertices(mesh.codim_vertices(), Eigen::all), // - Eigen::MatrixXi(), Eigen::MatrixXi(), inflation_radius); + Eigen::MatrixXi(), Eigen::MatrixXi(), + mesh.has_thickness() ? box_thickness(mesh.codim_vertices()) + : Eigen::VectorXd()); broad_phase->detect_vertex_vertex_candidates(vv_candidates); for (auto& [vi, vj] : vv_candidates) { @@ -100,12 +146,13 @@ void Candidates::build( // Extract the vertices of the codim. edges Eigen::MatrixXd CE_V; // vertices of codim. edges Eigen::MatrixXi CE; // codim. edges (indices into CEV) + Eigen::VectorXi J; // map from CE_V to vertices { - Eigen::VectorXi _I, _J; // unused mappings + Eigen::VectorXi _I; // unused mapping igl::remove_unreferenced( vertices, pad_edges(mesh.edges()(mesh.codim_edges(), Eigen::all)), CE_V, - CE, _I, _J); + CE, _I, J); CE = unpad_edges(CE); } @@ -116,12 +163,19 @@ void Candidates::build( CE.array() += nCV; // Offset indices to account for codim. vertices + Eigen::VectorXd V_thickness; + if (mesh.has_thickness()) { + V_thickness.resize(V.rows()); + V_thickness.head(nCV) = box_thickness(mesh.codim_vertices()); + V_thickness.tail(CE_V.rows()) = box_thickness(J); + } + // TODO: Can we reuse the broad phase from above? broad_phase->clear(); // Ignore c-edge to c-edge and c-vertex to c-vertex broad_phase->can_vertices_collide = make_codim_cross_filter(nCV) & mesh.can_collide; - broad_phase->build(V, CE, Eigen::MatrixXi(), inflation_radius); + build_broad_phase(V, CE, Eigen::MatrixXi(), V_thickness); broad_phase->detect_edge_vertex_candidates(ev_candidates); for (auto& [ei, vi] : ev_candidates) { @@ -134,7 +188,8 @@ void Candidates::build( // Planes to vertices: for (const auto& plane : mesh.planes) { for (index_t vi = 0; vi < mesh.num_vertices(); ++vi) { - if (plane.signedDistance(vertices.row(vi)) < inflation_radius) { + if (plane.signedDistance(vertices.row(vi)) + < inflation_radius + mesh.vertex_thickness(vi)) { pv_candidates.emplace_back(plane, vi); } } @@ -160,18 +215,37 @@ void Candidates::build( clear(); + // See the static build() for the thickness. + const Eigen::VectorXd box_thickness = + mesh.has_thickness() ? vertex_box_thickness(mesh) : Eigen::VectorXd(); + const auto build_broad_phase = [&](Eigen::ConstRef V_t0, + Eigen::ConstRef V_t1, + Eigen::ConstRef E, + Eigen::ConstRef F, + Eigen::ConstRef V_t) { + if (V_t.size() == 0) { + broad_phase->build(V_t0, V_t1, E, F, inflation_radius); + return; + } + AABBs boxes; + build_vertex_boxes(V_t0, V_t1, boxes, inflation_radius, V_t); + broad_phase->build(boxes, E, F, uint8_t(dim)); + }; + broad_phase->can_vertices_collide = mesh.can_collide; - broad_phase->build( - vertices_t0, vertices_t1, mesh.edges(), mesh.faces(), inflation_radius); + build_broad_phase( + vertices_t0, vertices_t1, mesh.edges(), mesh.faces(), box_thickness); broad_phase->detect_collision_candidates(*this); // Codim. vertices to codim. vertices: if (mesh.num_codim_vertices()) { broad_phase->clear(); - broad_phase->build( + build_broad_phase( vertices_t0(mesh.codim_vertices(), Eigen::all), vertices_t1(mesh.codim_vertices(), Eigen::all), // - Eigen::MatrixXi(), Eigen::MatrixXi(), inflation_radius); + Eigen::MatrixXi(), Eigen::MatrixXi(), + mesh.has_thickness() ? box_thickness(mesh.codim_vertices()) + : Eigen::VectorXd()); broad_phase->detect_vertex_vertex_candidates(vv_candidates); for (auto& [vi, vj] : vv_candidates) { @@ -188,8 +262,9 @@ void Candidates::build( // Extract the vertices of the codim. edges Eigen::MatrixXd CE_V_t0, CE_V_t1; // vertices of codim. edges Eigen::MatrixXi CE; // codim. edges (indices into CEV) + Eigen::VectorXi J; // map from CE_V to vertices { - Eigen::VectorXi _I, J; + Eigen::VectorXi _I; // unused mapping igl::remove_unreferenced( vertices_t0, pad_edges(mesh.edges()(mesh.codim_edges(), Eigen::all)), @@ -210,12 +285,19 @@ void Candidates::build( CE.array() += nCV; // Offset indices to account for codim. vertices + Eigen::VectorXd V_thickness; + if (mesh.has_thickness()) { + V_thickness.resize(V_t0.rows()); + V_thickness.head(nCV) = box_thickness(mesh.codim_vertices()); + V_thickness.tail(CE_V_t0.rows()) = box_thickness(J); + } + // TODO: Can we reuse the broad phase from above? broad_phase->clear(); // Ignore c-edge to c-edge and c-vertex to c-vertex broad_phase->can_vertices_collide = make_codim_cross_filter(nCV) & mesh.can_collide; - broad_phase->build(V_t0, V_t1, CE, Eigen::MatrixXi(), inflation_radius); + build_broad_phase(V_t0, V_t1, CE, Eigen::MatrixXi(), V_thickness); broad_phase->detect_edge_vertex_candidates(ev_candidates); for (auto& [ei, vi] : ev_candidates) { @@ -230,7 +312,8 @@ void Candidates::build( for (index_t vi = 0; vi < mesh.num_vertices(); ++vi) { const double d_t0 = plane.signedDistance(vertices_t0.row(vi)); const double d_t1 = plane.signedDistance(vertices_t1.row(vi)); - if (d_t0 < inflation_radius || d_t1 < inflation_radius) { + const double r = inflation_radius + mesh.vertex_thickness(vi); + if (d_t0 < r || d_t1 < r) { pv_candidates.emplace_back(plane, vi); } } @@ -241,11 +324,14 @@ bool Candidates::is_step_collision_free( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance, + const std::optional min_distance, const NarrowPhaseCCD& narrow_phase_ccd) const { assert(vertices_t0.rows() == mesh.num_vertices()); assert(vertices_t1.rows() == mesh.num_vertices()); + if (min_distance.has_value()) { + warn_min_distance_deprecated(); + } // Narrow phase for (size_t i = 0; i < size(); i++) { @@ -254,7 +340,8 @@ bool Candidates::is_step_collision_free( return candidate.ccd( candidate.dof(vertices_t0, mesh.edges(), mesh.faces()), candidate.dof(vertices_t1, mesh.edges(), mesh.faces()), // - toi, min_distance, /*tmax=*/1.0, narrow_phase_ccd); + toi, pair_min_distance(min_distance, candidate, mesh), + /*tmax=*/1.0, narrow_phase_ccd); }); if (is_collision) { @@ -269,12 +356,15 @@ double Candidates::compute_collision_free_stepsize( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance, + const std::optional min_distance, const NarrowPhaseCCD& narrow_phase_ccd) const { assert(vertices_t0.rows() == mesh.num_vertices()); assert(vertices_t1.rows() == mesh.num_vertices()); IPC_TOOLKIT_PROFILE_BLOCK("Candidates::compute_collision_free_stepsize"); + if (min_distance.has_value()) { + warn_min_distance_deprecated(); + } if (empty()) { return 1; // No possible collisions, so can take full step. @@ -290,7 +380,8 @@ double Candidates::compute_collision_free_stepsize( return candidate.ccd( candidate.dof(vertices_t0, mesh.edges(), mesh.faces()), candidate.dof(vertices_t1, mesh.edges(), mesh.faces()), // - toi, min_distance, tmax, narrow_phase_ccd); + toi, pair_min_distance(min_distance, candidate, mesh), tmax, + narrow_phase_ccd); }); if (are_colliding) { @@ -359,7 +450,7 @@ double Candidates::compute_cfl_stepsize( Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, const double dhat, - const double min_distance, + const std::optional min_distance, BroadPhase* broad_phase, const NarrowPhaseCCD& narrow_phase_ccd) const { @@ -385,22 +476,26 @@ Eigen::VectorXd Candidates::compute_per_vertex_safe_distances( const CollisionMesh& mesh, Eigen::ConstRef vertices, const double inflation_radius, - const double min_distance) const + const std::optional min_distance) const { - assert(inflation_radius >= min_distance); + if (min_distance.has_value()) { + warn_min_distance_deprecated(); + } + assert(inflation_radius >= min_distance.value_or(0.0)); // Initialize atomic distances for each vertex std::vector> min_distances(mesh.num_vertices()); for (size_t i = 0; i < mesh.num_vertices(); ++i) { min_distances[i].store( - inflation_radius - min_distance, std::memory_order_relaxed); + inflation_radius - min_distance.value_or(0.0), + std::memory_order_relaxed); } tbb::parallel_for(size_t(0), size(), [&](size_t i) { visit(i, [&](const auto& candidate) { const double d = sqrt(candidate.compute_distance( vertices, mesh.edges(), mesh.faces())) - - min_distance; + - pair_min_distance(min_distance, candidate, mesh); // Compute the distance for each vertex in the candidate for (auto vid : candidate.vertex_ids(mesh.edges(), mesh.faces())) { diff --git a/src/ipc/candidates/candidates.hpp b/src/ipc/candidates/candidates.hpp index f12854068..83c684498 100644 --- a/src/ipc/candidates/candidates.hpp +++ b/src/ipc/candidates/candidates.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -30,7 +31,7 @@ class Candidates { /// @brief Initialize the set of discrete collision detection candidates. /// @param mesh The surface of the collision mesh. /// @param vertices Surface vertex positions (rowwise). - /// @param inflation_radius Amount to inflate the bounding boxes. + /// @param inflation_radius Amount to inflate the bounding boxes. With a mesh thickness (CollisionMesh::set_thickness()), each box is inflated by its primitive's thickness on top of this. /// @param broad_phase Broad phase method to use. /// @param all_types Also detect the candidate types only ESP needs (see /// BroadPhase::detect_collision_candidates). Recorded in @@ -49,7 +50,7 @@ class Candidates { /// @param mesh The surface of the collision mesh. /// @param vertices_t0 Surface vertex starting positions (rowwise). /// @param vertices_t1 Surface vertex ending positions (rowwise). - /// @param inflation_radius Amount to inflate the bounding boxes. + /// @param inflation_radius Amount to inflate the bounding boxes. With a mesh thickness, each box is inflated by its primitive's thickness on top of this. /// @param broad_phase Broad phase method to use. void build( const CollisionMesh& mesh, @@ -160,14 +161,14 @@ class Candidates { /// @param mesh The collision mesh. /// @param vertices_t0 Surface vertex starting positions (rowwise). /// @param vertices_t1 Surface vertex ending positions (rowwise). - /// @param min_distance The minimum distance allowable between any two elements. + /// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @param narrow_phase_ccd The narrow phase CCD algorithm to use. /// @returns True if any collisions occur. bool is_step_collision_free( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; @@ -176,14 +177,14 @@ class Candidates { /// @param mesh The collision mesh. /// @param vertices_t0 Surface vertex starting positions (rowwise). Assumed to be intersection free. /// @param vertices_t1 Surface vertex ending positions (rowwise). - /// @param min_distance The minimum distance allowable between any two elements. + /// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @param narrow_phase_ccd The narrow phase CCD algorithm to use. /// @returns A step-size \f$\in [0, 1]\f$ that is collision free. A value of 1.0 if a full step and 0.0 is no step. double compute_collision_free_stepsize( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; @@ -202,7 +203,7 @@ class Candidates { /// @param vertices_t0 Surface vertex starting positions (rowwise). /// @param vertices_t1 Surface vertex ending positions (rowwise). /// @param dhat Barrier activation distance. - /// @param min_distance The minimum distance allowable between any two elements. + /// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @param broad_phase The broad phase algorithm to use. /// @param narrow_phase_ccd The narrow phase CCD algorithm to use. /// @returns A step-size \f$\in [0, 1]\f$ that is collision free. @@ -211,7 +212,7 @@ class Candidates { Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, const double dhat, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr, const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; @@ -221,13 +222,13 @@ class Candidates { /// @param mesh The collision mesh. /// @param vertices Collision mesh vertex positions (rowwise). /// @param inflation_radius The inflation radius used to build the candidates. - /// @param min_distance The minimum distance allowable between any two elements. + /// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @return A vector of minimum distances, one for each vertex. Eigen::VectorXd compute_per_vertex_safe_distances( const CollisionMesh& mesh, Eigen::ConstRef vertices, const double inflation_radius, - const double min_distance = 0.0) const; + const std::optional min_distance = std::nullopt) const; // clang-format off // /// @brief Compute the maximum distance every vertex can move (independently) without colliding with any other element. diff --git a/src/ipc/candidates/edge_edge.hpp b/src/ipc/candidates/edge_edge.hpp index 4ad12e13c..52e8dae40 100644 --- a/src/ipc/candidates/edge_edge.hpp +++ b/src/ipc/candidates/edge_edge.hpp @@ -75,6 +75,20 @@ class EdgeEdgeCandidate : public StencilMixin { const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; + /// @brief Get the sum of the two primitives' thicknesses. + /// @see CollisionMesh::set_thickness() + /// @param vertex_thickness Thickness of each vertex. + /// @param edge_thickness Thickness of each edge. + /// @param face_thickness Thickness of each face. + /// @return How much farther apart than the minimum distance the pair rests. + double thickness( + Eigen::ConstRef /*vertex_thickness*/, + Eigen::ConstRef edge_thickness, + Eigen::ConstRef /*face_thickness*/) const + { + return edge_thickness[edge0_id] + edge_thickness[edge1_id]; + } + // ------------------------------------------------------------------------ /// @brief The distance type known a priori; AUTO for a plain candidate. diff --git a/src/ipc/candidates/edge_vertex.hpp b/src/ipc/candidates/edge_vertex.hpp index 74c147c42..a91102d28 100644 --- a/src/ipc/candidates/edge_vertex.hpp +++ b/src/ipc/candidates/edge_vertex.hpp @@ -72,6 +72,20 @@ class EdgeVertexCandidate : public StencilMixin { const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; + /// @brief Get the sum of the two primitives' thicknesses. + /// @see CollisionMesh::set_thickness() + /// @param vertex_thickness Thickness of each vertex. + /// @param edge_thickness Thickness of each edge. + /// @param face_thickness Thickness of each face. + /// @return How much farther apart than the minimum distance the pair rests. + double thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef edge_thickness, + Eigen::ConstRef /*face_thickness*/) const + { + return edge_thickness[edge_id] + vertex_thickness[vertex_id]; + } + // ------------------------------------------------------------------------ constexpr static PointEdgeDistanceType known_dtype() diff --git a/src/ipc/candidates/face_vertex.hpp b/src/ipc/candidates/face_vertex.hpp index cb7b5a8d9..e9bbb5622 100644 --- a/src/ipc/candidates/face_vertex.hpp +++ b/src/ipc/candidates/face_vertex.hpp @@ -77,6 +77,20 @@ class FaceVertexCandidate : public StencilMixin { const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; + /// @brief Get the sum of the two primitives' thicknesses. + /// @see CollisionMesh::set_thickness() + /// @param vertex_thickness Thickness of each vertex. + /// @param edge_thickness Thickness of each edge. + /// @param face_thickness Thickness of each face. + /// @return How much farther apart than the minimum distance the pair rests. + double thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef /*edge_thickness*/, + Eigen::ConstRef face_thickness) const + { + return face_thickness[face_id] + vertex_thickness[vertex_id]; + } + // ------------------------------------------------------------------------ constexpr static PointTriangleDistanceType known_dtype() diff --git a/src/ipc/candidates/min_distance.hpp b/src/ipc/candidates/min_distance.hpp new file mode 100644 index 000000000..e0974df0f --- /dev/null +++ b/src/ipc/candidates/min_distance.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include +#include + +#include +#include + +namespace ipc { + +/// @brief Warn (at most once per process) that min_distance is deprecated. +inline void warn_min_distance_deprecated() +{ + static std::once_flag flag; + std::call_once(flag, [] { + logger().warn( + "The min_distance parameter is deprecated: give the collision mesh " + "a thickness with CollisionMesh::set_thickness() instead. A " + "min_distance overrides the mesh's thickness."); + }); +} + +/// @brief Warn (at most once per process) that dmin is deprecated. +inline void warn_dmin_deprecated() +{ + static std::once_flag flag; + std::call_once(flag, [] { + logger().warn( + "The dmin parameter is deprecated: give the collision mesh a " + "thickness with CollisionMesh::set_thickness() instead. A dmin " + "overrides the mesh's thickness."); + }); +} + +/// @brief The distance a pair of primitives must keep. +/// @param min_distance The deprecated uniform minimum distance, which overrides the mesh's thickness when given. +/// @param candidate The pair of primitives. +/// @param mesh The collision mesh whose primitives the candidate indexes. +/// @return min_distance if given, else the pair's thickness (zero without one). +template +double pair_min_distance( + const std::optional& min_distance, + const Candidate& candidate, + const CollisionMesh& mesh) +{ + return min_distance.has_value() ? *min_distance + : mesh.pair_thickness(candidate); +} + +} // namespace ipc diff --git a/src/ipc/candidates/plane_vertex.hpp b/src/ipc/candidates/plane_vertex.hpp index c26d915e2..58dfef592 100644 --- a/src/ipc/candidates/plane_vertex.hpp +++ b/src/ipc/candidates/plane_vertex.hpp @@ -73,6 +73,20 @@ class PlaneVertexCandidate : public StencilMixin { const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; + /// @brief Get the sum of the two primitives' thicknesses. + /// @see CollisionMesh::set_thickness() + /// @param vertex_thickness Thickness of each vertex. + /// @param edge_thickness Thickness of each edge. + /// @param face_thickness Thickness of each face. + /// @return How much farther apart than the minimum distance the pair rests. + double thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef /*edge_thickness*/, + Eigen::ConstRef /*face_thickness*/) const + { + return vertex_thickness[vertex_id]; // The plane has no thickness. + } + /// @brief The plane of the candidate. Eigen::Hyperplane plane; diff --git a/src/ipc/candidates/vertex_vertex.hpp b/src/ipc/candidates/vertex_vertex.hpp index 70ec2b53d..38e5eb16b 100644 --- a/src/ipc/candidates/vertex_vertex.hpp +++ b/src/ipc/candidates/vertex_vertex.hpp @@ -67,6 +67,20 @@ class VertexVertexCandidate : public StencilMixin { const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD) const; + /// @brief Get the sum of the two primitives' thicknesses. + /// @see CollisionMesh::set_thickness() + /// @param vertex_thickness Thickness of each vertex. + /// @param edge_thickness Thickness of each edge. + /// @param face_thickness Thickness of each face. + /// @return How much farther apart than the minimum distance the pair rests. + double thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef /*edge_thickness*/, + Eigen::ConstRef /*face_thickness*/) const + { + return vertex_thickness[vertex0_id] + vertex_thickness[vertex1_id]; + } + // ------------------------------------------------------------------------ bool operator==(const VertexVertexCandidate& other) const; diff --git a/src/ipc/collision_mesh.cpp b/src/ipc/collision_mesh.cpp index 11874f8e8..ef8be816a 100644 --- a/src/ipc/collision_mesh.cpp +++ b/src/ipc/collision_mesh.cpp @@ -392,17 +392,17 @@ void CollisionMesh::init_areas() } namespace { - void check_areas( - Eigen::ConstRef areas, + void check_non_negative( + Eigen::ConstRef values, const Eigen::Index expected_size, const char* name) { - if (areas.size() != expected_size) { + if (values.size() != expected_size) { log_and_throw_error( - "{} has size {} but expected {}", name, areas.size(), + "{} has size {} but expected {}", name, values.size(), expected_size); } - if (!areas.allFinite() || (areas.array() < 0).any()) { + if (!values.allFinite() || (values.array() < 0).any()) { log_and_throw_error( "{} must be finite and non-negative", name); } @@ -412,7 +412,7 @@ namespace { void CollisionMesh::set_vertex_areas( Eigen::ConstRef vertex_areas) { - check_areas(vertex_areas, num_vertices(), "vertex_areas"); + check_non_negative(vertex_areas, num_vertices(), "vertex_areas"); if (are_area_jacobians_initialized()) { log_and_throw_error( "Cannot set custom vertex areas after init_area_jacobians(): the area Jacobians would no longer match."); @@ -423,7 +423,7 @@ void CollisionMesh::set_vertex_areas( void CollisionMesh::set_edge_areas(Eigen::ConstRef edge_areas) { - check_areas(edge_areas, num_edges(), "edge_areas"); + check_non_negative(edge_areas, num_edges(), "edge_areas"); if (are_area_jacobians_initialized()) { log_and_throw_error( "Cannot set custom edge areas after init_area_jacobians(): the area Jacobians would no longer match."); @@ -432,6 +432,86 @@ void CollisionMesh::set_edge_areas(Eigen::ConstRef edge_areas) m_has_custom_areas = true; } +// ============================================================================ + +void CollisionMesh::set_thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef edge_thickness, + Eigen::ConstRef face_thickness) +{ + check_non_negative(vertex_thickness, num_vertices(), "vertex_thickness"); + check_non_negative(edge_thickness, num_edges(), "edge_thickness"); + check_non_negative(face_thickness, num_faces(), "face_thickness"); + m_vertex_thickness = vertex_thickness; + m_edge_thickness = edge_thickness; + m_face_thickness = face_thickness; + m_max_thickness = std::max( + { m_vertex_thickness.size() ? m_vertex_thickness.maxCoeff() : 0.0, + m_edge_thickness.size() ? m_edge_thickness.maxCoeff() : 0.0, + m_face_thickness.size() ? m_face_thickness.maxCoeff() : 0.0 }); +} + +void CollisionMesh::set_vertex_thickness( + Eigen::ConstRef vertex_thickness) +{ + check_non_negative(vertex_thickness, num_vertices(), "vertex_thickness"); + + Eigen::VectorXd edge_thickness(num_edges()); + for (int i = 0; i < m_edges.rows(); i++) { + edge_thickness[i] = std::max( + vertex_thickness[m_edges(i, 0)], vertex_thickness[m_edges(i, 1)]); + } + + Eigen::VectorXd face_thickness(num_faces()); + for (int i = 0; i < m_faces.rows(); i++) { + face_thickness[i] = std::max( + { vertex_thickness[m_faces(i, 0)], vertex_thickness[m_faces(i, 1)], + vertex_thickness[m_faces(i, 2)] }); + } + + set_thickness(vertex_thickness, edge_thickness, face_thickness); +} + +void CollisionMesh::set_face_thickness( + Eigen::ConstRef face_thickness, + Eigen::ConstRef vertex_thickness) +{ + check_non_negative(face_thickness, num_faces(), "face_thickness"); + + Eigen::VectorXd vertex = Eigen::VectorXd::Zero(num_vertices()); + if (vertex_thickness.size() > 0) { + check_non_negative( + vertex_thickness, num_vertices(), "vertex_thickness"); + vertex = vertex_thickness; + } + + // An edge no face reaches stays negative until the pass below. + Eigen::VectorXd edge = Eigen::VectorXd::Constant(num_edges(), -1); + for (int i = 0; i < m_faces.rows(); i++) { + for (int j = 0; j < 3; j++) { + vertex[m_faces(i, j)] = + std::max(vertex[m_faces(i, j)], face_thickness[i]); + edge[m_faces_to_edges(i, j)] = + std::max(edge[m_faces_to_edges(i, j)], face_thickness[i]); + } + } + for (int i = 0; i < m_edges.rows(); i++) { + if (edge[i] < 0) { + edge[i] = std::max(vertex[m_edges(i, 0)], vertex[m_edges(i, 1)]); + } + } + + set_thickness(vertex, edge, face_thickness); +} + +void CollisionMesh::clear_thickness() +{ + m_vertex_thickness.resize(0); + m_edge_thickness.resize(0); + m_face_thickness.resize(0); + m_max_thickness = 0; +} + void CollisionMesh::init_area_jacobians() { if (m_has_custom_areas) { diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp index e5559abf7..ffdbad051 100644 --- a/src/ipc/collision_mesh.hpp +++ b/src/ipc/collision_mesh.hpp @@ -399,6 +399,112 @@ class CollisionMesh { && m_edge_area_jacobian.size() == num_edges(); } + // ----------------------------------------------------------------------- + // Thickness + + /// @brief Determine if the primitives have a thickness (see set_thickness()). + bool has_thickness() const { return m_vertex_thickness.size() > 0; } + + /// @brief Get the thickness of a vertex (zero without a thickness). + /// @param vi Vertex ID. + /// @return Thickness of vertex vi. + double vertex_thickness(const index_t vi) const + { + return has_thickness() ? m_vertex_thickness[vi] : 0.0; + } + + /// @brief Get the thickness of an edge (zero without a thickness). + /// @param ei Edge ID. + /// @return Thickness of edge ei. + double edge_thickness(const index_t ei) const + { + return has_thickness() ? m_edge_thickness[ei] : 0.0; + } + + /// @brief Get the thickness of a face (zero without a thickness). + /// @param fi Face ID. + /// @return Thickness of face fi. + double face_thickness(const index_t fi) const + { + return has_thickness() ? m_face_thickness[fi] : 0.0; + } + + /// @brief Get the thickness of the vertices (empty without a thickness). + const Eigen::VectorXd& vertex_thicknesses() const + { + return m_vertex_thickness; + } + + /// @brief Get the thickness of the edges (empty without a thickness). + const Eigen::VectorXd& edge_thicknesses() const { return m_edge_thickness; } + + /// @brief Get the thickness of the faces (empty without a thickness). + const Eigen::VectorXd& face_thicknesses() const { return m_face_thickness; } + + /// @brief Get the largest thickness of any primitive (zero without a thickness). + double max_thickness() const { return m_max_thickness; } + + /// @brief Get the sum of a pair's two primitives' thicknesses (zero without a thickness). + /// @tparam Candidate A candidate or collision type with a thickness(vertex_thickness, edge_thickness, face_thickness) method. + /// @param candidate The pair of primitives, whose IDs index this mesh. + /// @return How much farther apart than the minimum distance the pair rests. + template + double pair_thickness(const Candidate& candidate) const + { + return has_thickness() + ? candidate.thickness( + m_vertex_thickness, m_edge_thickness, m_face_thickness) + : 0.0; + } + + /// @brief Give every vertex, edge, and face a thickness. + /// + /// Each primitive is offset outward by its thickness, so two primitives + /// touch when they are the sum of their thicknesses apart. Collisions + /// built from this mesh add that sum to their minimum distance + /// (NormalCollision::dmin), the broad phase inflates each primitive's box + /// by its thickness, and CCD stops each pair at that sum plus the minimum + /// distance. A sheet of total thickness h, for example, gives each of its + /// primitives h / 2. + /// + /// @note Collisions and candidates built before this call keep their old distances. + /// @note CCD stops each edge-edge and face-vertex pair at its own thickness, while the collision it reduces to (e.g., vertex-vertex at a face's corner) keeps that sub-primitive's thickness. Keep each vertex and edge no thicker than the thickest face or edge containing it, as set_vertex_thickness() and set_face_thickness() do, or a collision-free step can end inside such a collision's minimum distance. + /// @param vertex_thickness Thickness of each vertex (size num_vertices()). Must be finite and non-negative, as must the others. + /// @param edge_thickness Thickness of each edge (size num_edges()). + /// @param face_thickness Thickness of each face (size num_faces()). + void set_thickness( + Eigen::ConstRef vertex_thickness, + Eigen::ConstRef edge_thickness, + Eigen::ConstRef face_thickness); + + /// @brief Give each vertex a thickness, which its edges and faces inherit. + /// + /// An edge or face takes the largest thickness of its vertices. This suits + /// thickness given per vertex, such as rods of varying radius. + /// + /// @see set_thickness() + /// @param vertex_thickness Thickness of each vertex (size num_vertices()). + void + set_vertex_thickness(Eigen::ConstRef vertex_thickness); + + /// @brief Give each face a thickness, which its edges and vertices inherit. + /// + /// An edge or vertex takes the largest thickness of the faces it bounds, so + /// a seam between a thick and a thin panel is as thick as the thick one. An + /// edge on no face takes the larger thickness of its two vertices. This + /// suits thickness that is piecewise constant over the surface, such as + /// garment panels of different fabrics. + /// + /// @see set_thickness() + /// @param face_thickness Thickness of each face (size num_faces()). + /// @param vertex_thickness Least thickness of each vertex (size num_vertices()), or empty for zero. It is the whole thickness of a vertex on no face. + void set_face_thickness( + Eigen::ConstRef face_thickness, + Eigen::ConstRef vertex_thickness = Eigen::VectorXd()); + + /// @brief Remove the primitives' thickness. + void clear_thickness(); + // ----------------------------------------------------------------------- /// @brief Construct a vector of bools indicating whether each vertex is on the surface. @@ -536,6 +642,15 @@ class CollisionMesh { std::vector> m_vertex_area_jacobian; /// @brief The rows of the Jacobian of the edge areas vector. std::vector> m_edge_area_jacobian; + + /// @brief Thickness of each vertex (empty without a thickness). + Eigen::VectorXd m_vertex_thickness; + /// @brief Thickness of each edge (empty without a thickness). + Eigen::VectorXd m_edge_thickness; + /// @brief Thickness of each face (empty without a thickness). + Eigen::VectorXd m_face_thickness; + /// @brief Largest thickness of any primitive. + double m_max_thickness = 0; }; } // namespace ipc diff --git a/src/ipc/collisions/normal/normal_collisions.cpp b/src/ipc/collisions/normal/normal_collisions.cpp index 89c8382ac..545aad90d 100644 --- a/src/ipc/collisions/normal/normal_collisions.cpp +++ b/src/ipc/collisions/normal/normal_collisions.cpp @@ -1,5 +1,6 @@ #include "normal_collisions.hpp" +#include #include #include #include @@ -19,13 +20,14 @@ void NormalCollisions::build( const CollisionMesh& mesh, Eigen::ConstRef vertices, const double dhat, - const double dmin, + const std::optional dmin, BroadPhase* broad_phase) { assert(vertices.rows() == mesh.num_vertices()); IPC_TOOLKIT_PROFILE_BLOCK("NormalCollisions::build"); - const double inflation_radius = 0.5 * (dhat + dmin); + // Candidates::build() adds each primitive's thickness on top of this. + const double inflation_radius = 0.5 * (dhat + dmin.value_or(0.0)); Candidates candidates; candidates.build(mesh, vertices, inflation_radius, broad_phase); @@ -38,7 +40,7 @@ void NormalCollisions::build( const CollisionMesh& mesh, Eigen::ConstRef vertices, const double dhat, - const double dmin) + const std::optional dmin) { assert(vertices.rows() == mesh.num_vertices()); IPC_TOOLKIT_PROFILE_BLOCK("NormalCollisions::build(candidates)"); @@ -53,13 +55,30 @@ void NormalCollisions::build( "NormalCollisions::build: 3D candidates must be built with all_types = false"); } + if (dmin.has_value()) { + warn_dmin_deprecated(); + } + clear(); + // The minimum distance of a pair: the deprecated dmin if given, else the + // sum of the two primitives' thicknesses. + const auto pair_dmin = [&dmin](const double thickness) { + return dmin.has_value() ? *dmin : thickness; + }; + // Cull the candidates by measuring the distance and dropping those that are - // greater than dhat. - auto is_active = [offset_sqr = sqr(dmin + dhat)](double distance_sqr) { - return distance_sqr < offset_sqr; + // greater than dhat past the pair's minimum distance. + auto is_active = [&](double distance_sqr, double thickness) { + return distance_sqr < sqr(pair_dmin(thickness) + dhat); }; + // The improved max approximation's corrections are culled with the + // thickest pair's band. That keeps every correction a thinner pair's + // collision needs, and a correction outside its own band adds nothing + // because its potential is zero there. + auto is_active_correction = + [offset_sqr = sqr(pair_dmin(2 * mesh.max_thickness()) + dhat)]( + double distance_sqr) { return distance_sqr < offset_sqr; }; tbb::enumerable_thread_specific storage( use_area_weighting(), enable_shape_derivatives(), @@ -93,7 +112,7 @@ void NormalCollisions::build( if (!candidates.ev_candidates.empty()) { // Convert edge-vertex to vertex-vertex const auto vv_candidates = candidates.edge_vertex_to_vertex_vertex( - mesh, vertices, is_active); + mesh, vertices, is_active_correction); tbb::parallel_for(size_t(0), vv_candidates.size(), [&](size_t i) { storage.local() @@ -104,8 +123,8 @@ void NormalCollisions::build( if (!candidates.ee_candidates.empty()) { // Convert edge-edge to edge-vertex - const auto ev_candidates = - candidates.edge_edge_to_edge_vertex(mesh, vertices, is_active); + const auto ev_candidates = candidates.edge_edge_to_edge_vertex( + mesh, vertices, is_active_correction); tbb::parallel_for(size_t(0), ev_candidates.size(), [&](size_t i) { storage.local().add_edge_edge_negative_edge_vertex_collision( @@ -116,7 +135,7 @@ void NormalCollisions::build( if (!candidates.fv_candidates.empty()) { // Convert face-vertex to edge-vertex const auto ev_candidates = candidates.face_vertex_to_edge_vertex( - mesh, vertices, is_active); + mesh, vertices, is_active_correction); tbb::parallel_for(size_t(0), ev_candidates.size(), [&](size_t i) { storage.local().add_face_vertex_negative_edge_vertex_collision( @@ -125,7 +144,7 @@ void NormalCollisions::build( // Convert face-vertex to vertex-vertex const auto vv_candidates = candidates.face_vertex_to_vertex_vertex( - mesh, vertices, is_active); + mesh, vertices, is_active_correction); tbb::parallel_for(size_t(0), vv_candidates.size(), [&](size_t i) { storage.local() @@ -144,7 +163,7 @@ void NormalCollisions::build( continue; } const double d = plane.signedDistance(vertices.row(vertex_id)); - if (d < dhat + dmin) { + if (d < dhat + pair_dmin(mesh.vertex_thickness(vertex_id))) { const double weight = use_area_weighting() ? mesh.vertex_area(vertex_id) : 1; @@ -162,10 +181,16 @@ void NormalCollisions::build( // logger().debug(to_string(mesh, vertices)); - for (size_t ci = 0; ci < size(); ci++) { - NormalCollision& collision = (*this)[ci]; - collision.dmin = dmin; - } + const auto set_dmin = [&](auto& collisions) { + for (auto& collision : collisions) { + collision.dmin = pair_dmin(mesh.pair_thickness(collision)); + } + }; + set_dmin(vv_collisions); + set_dmin(ev_collisions); + set_dmin(ee_collisions); + set_dmin(fv_collisions); + set_dmin(pv_collisions); } void NormalCollisions::set_use_area_weighting(const bool use_area_weighting) diff --git a/src/ipc/collisions/normal/normal_collisions.hpp b/src/ipc/collisions/normal/normal_collisions.hpp index b9cbebe2e..4b3b220cb 100644 --- a/src/ipc/collisions/normal/normal_collisions.hpp +++ b/src/ipc/collisions/normal/normal_collisions.hpp @@ -11,6 +11,7 @@ #include +#include #include namespace ipc { @@ -37,16 +38,17 @@ class NormalCollisions { NormalCollisions() = default; /// @brief Initialize the set of collisions used to compute the barrier potential. + /// @note Each collision's minimum distance (NormalCollision::dmin) is the sum of its two primitives' thicknesses (CollisionMesh::set_thickness()), and it is active within dhat of that. /// @param mesh The collision mesh. /// @param vertices Vertices of the collision mesh. /// @param dhat The activation distance of the barrier. - /// @param dmin Minimum distance. + /// @param dmin Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every collision keeps this minimum distance and the mesh's thickness is ignored; otherwise each keeps its two primitives' thicknesses (zero without a thickness). /// @param broad_phase Broad-phase method to use. void build( const CollisionMesh& mesh, Eigen::ConstRef vertices, const double dhat, - const double dmin = 0, + const std::optional dmin = std::nullopt, BroadPhase* broad_phase = nullptr); /// @brief Initialize the set of collisions used to compute the barrier potential. @@ -54,7 +56,8 @@ class NormalCollisions { /// @param mesh The collision mesh. /// @param vertices Vertices of the collision mesh. /// @param dhat The activation distance of the barrier. - /// @param dmin Minimum distance. + /// @param dmin Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every collision keeps this minimum distance and the mesh's thickness is ignored; otherwise each keeps its two primitives' thicknesses (zero without a thickness). + /// @note With a mesh thickness, the candidates must be built from the same mesh so their boxes cover it (see the other overload). /// @throws std::invalid_argument In 3D, if the candidates were built with /// all_types = true (see Candidates::has_all_types). void build( @@ -62,7 +65,7 @@ class NormalCollisions { const CollisionMesh& mesh, Eigen::ConstRef vertices, const double dhat, - const double dmin = 0); + const std::optional dmin = std::nullopt); // ------------------------------------------------------------------------ diff --git a/src/ipc/collisions/normal/normal_collisions_builder.cpp b/src/ipc/collisions/normal/normal_collisions_builder.cpp index 94e98b749..feab1988c 100644 --- a/src/ipc/collisions/normal/normal_collisions_builder.cpp +++ b/src/ipc/collisions/normal/normal_collisions_builder.cpp @@ -11,6 +11,80 @@ namespace ipc { +namespace { + // Summed thickness of the primitives a collision keeps once its distance + // type reduces it to a sub-stencil (e.g., a face-vertex pair closest to an + // edge becomes an edge-vertex collision). + + double point_edge_thickness( + const CollisionMesh& mesh, + const EdgeVertexCandidate& candidate, + const PointEdgeDistanceType dtype) + { + const auto& [ei, vi] = candidate; + switch (dtype) { + case PointEdgeDistanceType::P_E0: + case PointEdgeDistanceType::P_E1: + return mesh.pair_thickness( + VertexVertexCandidate(vi, mesh.edges()(ei, int(dtype)))); + default: + return mesh.pair_thickness(candidate); + } + } + + double edge_edge_thickness( + const CollisionMesh& mesh, + const EdgeEdgeCandidate& candidate, + const EdgeEdgeDistanceType dtype) + { + const auto& [eai, ebi] = candidate; + const auto [ea0i, ea1i, eb0i, eb1i] = + candidate.vertex_ids(mesh.edges(), mesh.faces()); + switch (dtype) { + case EdgeEdgeDistanceType::EA0_EB0: + return mesh.pair_thickness(VertexVertexCandidate(ea0i, eb0i)); + case EdgeEdgeDistanceType::EA0_EB1: + return mesh.pair_thickness(VertexVertexCandidate(ea0i, eb1i)); + case EdgeEdgeDistanceType::EA1_EB0: + return mesh.pair_thickness(VertexVertexCandidate(ea1i, eb0i)); + case EdgeEdgeDistanceType::EA1_EB1: + return mesh.pair_thickness(VertexVertexCandidate(ea1i, eb1i)); + case EdgeEdgeDistanceType::EA_EB0: + return mesh.pair_thickness(EdgeVertexCandidate(eai, eb0i)); + case EdgeEdgeDistanceType::EA_EB1: + return mesh.pair_thickness(EdgeVertexCandidate(eai, eb1i)); + case EdgeEdgeDistanceType::EA0_EB: + return mesh.pair_thickness(EdgeVertexCandidate(ebi, ea0i)); + case EdgeEdgeDistanceType::EA1_EB: + return mesh.pair_thickness(EdgeVertexCandidate(ebi, ea1i)); + default: + return mesh.pair_thickness(candidate); + } + } + + double point_triangle_thickness( + const CollisionMesh& mesh, + const FaceVertexCandidate& candidate, + const PointTriangleDistanceType dtype) + { + const auto& [fi, vi] = candidate; + switch (dtype) { + case PointTriangleDistanceType::P_T0: + case PointTriangleDistanceType::P_T1: + case PointTriangleDistanceType::P_T2: + return mesh.pair_thickness( + VertexVertexCandidate(vi, mesh.faces()(fi, int(dtype)))); + case PointTriangleDistanceType::P_E0: + case PointTriangleDistanceType::P_E1: + case PointTriangleDistanceType::P_E2: + return mesh.pair_thickness(EdgeVertexCandidate( + mesh.faces_to_edges()(fi, int(dtype) - 3), vi)); + default: + return mesh.pair_thickness(candidate); + } + } +} // namespace + NormalCollisionsBuilder::NormalCollisionsBuilder( const bool _use_area_weighting, const bool _enable_shape_derivatives, @@ -29,7 +103,7 @@ void NormalCollisionsBuilder::add_vertex_vertex_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const VertexVertexCandidate& candidate, - const std::function& is_active) + const std::function& is_active) { const auto& [vi, vj] = candidate; @@ -41,7 +115,7 @@ void NormalCollisionsBuilder::add_vertex_vertex_collision( const double distance = point_point_distance(vertices.row(vi), vertices.row(vj)); - if (!is_active(distance)) { + if (!is_active(distance, mesh.pair_thickness(candidate))) { return; } @@ -76,7 +150,7 @@ void NormalCollisionsBuilder::add_edge_vertex_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const EdgeVertexCandidate& candidate, - const std::function& is_active) + const std::function& is_active) { const auto& [ei, vi] = candidate; @@ -90,7 +164,7 @@ void NormalCollisionsBuilder::add_edge_vertex_collision( const PointEdgeDistanceType dtype = point_edge_distance_type(v, e0, e1); const double distance_sqr = point_edge_distance(v, e0, e1, dtype); - if (!is_active(distance_sqr) + if (!is_active(distance_sqr, point_edge_thickness(mesh, candidate, dtype)) || (use_ogc && !ogc::is_edge_vertex_feasible( mesh, vertices, candidate, dtype))) { @@ -161,7 +235,7 @@ void NormalCollisionsBuilder::add_edge_edge_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const EdgeEdgeCandidate& candidate, - const std::function& is_active) + const std::function& is_active) { const auto& [eai, ebi] = candidate; @@ -183,7 +257,13 @@ void NormalCollisionsBuilder::add_edge_edge_collision( const double distance_sqr = edge_edge_distance(ea0, ea1, eb0, eb1, actual_dtype); - if (!is_active(distance_sqr) + // A nearly parallel pair is mollified into an edge-edge collision whatever + // its distance type, so it is active within either thickness. + const double thickness = std::max( + edge_edge_thickness(mesh, candidate, actual_dtype), + mesh.pair_thickness(candidate)); + + if (!is_active(distance_sqr, thickness) || (use_ogc && !ogc::is_edge_edge_feasible( mesh, vertices, candidate, actual_dtype))) { @@ -276,7 +356,7 @@ void NormalCollisionsBuilder::add_face_vertex_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const FaceVertexCandidate& candidate, - const std::function& is_active) + const std::function& is_active) { const auto& [fi, vi] = candidate; if (skip_obstacles && mesh.is_obstacle_vertex(vi)) { @@ -294,7 +374,8 @@ void NormalCollisionsBuilder::add_face_vertex_collision( point_triangle_distance_type(v, f0, f1, f2); const double distance_sqr = point_triangle_distance(v, f0, f1, f2, dtype); - if (!is_active(distance_sqr) + if (!is_active( + distance_sqr, point_triangle_thickness(mesh, candidate, dtype)) || (use_ogc && !ogc::is_face_vertex_feasible( mesh, vertices, candidate, dtype))) { diff --git a/src/ipc/collisions/normal/normal_collisions_builder.hpp b/src/ipc/collisions/normal/normal_collisions_builder.hpp index 56815ee13..f9748fcb6 100644 --- a/src/ipc/collisions/normal/normal_collisions_builder.hpp +++ b/src/ipc/collisions/normal/normal_collisions_builder.hpp @@ -16,6 +16,10 @@ namespace ipc { class NormalCollisionsBuilder { public: + // The add_*_collision() functions take is_active(distance_sqr, thickness), + // where thickness is the summed thickness of the two primitives that the + // collision would keep (CollisionMesh::set_thickness()). + NormalCollisionsBuilder( const bool use_area_weighting, const bool enable_shape_derivatives, @@ -26,25 +30,25 @@ class NormalCollisionsBuilder { const CollisionMesh& mesh, Eigen::ConstRef vertices, const VertexVertexCandidate& candidate, - const std::function& is_active); + const std::function& is_active); void add_edge_vertex_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const EdgeVertexCandidate& candidate, - const std::function& is_active); + const std::function& is_active); void add_edge_edge_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const EdgeEdgeCandidate& candidate, - const std::function& is_active); + const std::function& is_active); void add_face_vertex_collision( const CollisionMesh& mesh, Eigen::ConstRef vertices, const FaceVertexCandidate& candidate, - const std::function& is_active); + const std::function& is_active); // ------------------------------------------------------------------------ // Duplicate removal functions diff --git a/src/ipc/collisions/tangential/tangential_collision.cpp b/src/ipc/collisions/tangential/tangential_collision.cpp index a38c445a6..b888fb3ea 100644 --- a/src/ipc/collisions/tangential/tangential_collision.cpp +++ b/src/ipc/collisions/tangential/tangential_collision.cpp @@ -27,6 +27,7 @@ void TangentialCollision::init( closest_point = compute_closest_point(positions); tangent_basis = compute_tangent_basis(positions); normal_force_magnitude = normal_force; + dmin = collision.dmin; } } // namespace ipc diff --git a/src/ipc/collisions/tangential/tangential_collision.hpp b/src/ipc/collisions/tangential/tangential_collision.hpp index d5cf7376b..0457e48a4 100644 --- a/src/ipc/collisions/tangential/tangential_collision.hpp +++ b/src/ipc/collisions/tangential/tangential_collision.hpp @@ -93,6 +93,9 @@ class TangentialCollision : virtual public CollisionStencil { /// @brief Normal force magnitude double normal_force_magnitude = 0; + /// @brief Minimum distance of the normal collision this was built from (NormalCollision::dmin). + double dmin = 0; + /// @brief GCPCollision instance to compute normal force magnitude and its derivatives std::shared_ptr gcp_collision; diff --git a/src/ipc/ipc.cpp b/src/ipc/ipc.cpp index cf22f1685..805f1a078 100644 --- a/src/ipc/ipc.cpp +++ b/src/ipc/ipc.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -23,7 +24,7 @@ bool is_step_collision_free( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance, + const std::optional min_distance, BroadPhase* broad_phase, const NarrowPhaseCCD& narrow_phase_ccd) { @@ -35,7 +36,7 @@ bool is_step_collision_free( Candidates candidates; candidates.build( mesh, vertices_t0, vertices_t1, - /*inflation_radius=*/0.5 * min_distance, broad_phase); + /*inflation_radius=*/0.5 * min_distance.value_or(0.0), broad_phase); // Narrow phase return candidates.is_step_collision_free( @@ -48,7 +49,7 @@ double compute_collision_free_stepsize( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance, + const std::optional min_distance, BroadPhase* broad_phase, const NarrowPhaseCCD& narrow_phase_ccd) { @@ -63,7 +64,13 @@ double compute_collision_free_stepsize( } #ifdef IPC_TOOLKIT_WITH_CUDA - if (broad_phase->name() == "SweepAndTiniestQueue") { + // The GPU pipeline has no per-primitive thickness, so a mesh with one + // takes the general path below, which still uses this broad phase. + if (broad_phase->name() == "SweepAndTiniestQueue" + && (min_distance.has_value() || !mesh.has_thickness())) { + if (min_distance.has_value()) { + warn_min_distance_deprecated(); + } if (vertices_t0.cols() != 3) { log_and_throw_error( "Sweep and Tiniest Queue is only supported in 3D!"); @@ -94,8 +101,8 @@ double compute_collision_free_stepsize( // Broad phase Candidates candidates; candidates.build( - mesh, vertices_t0, vertices_t1, /*inflation_radius=*/0.5 * min_distance, - broad_phase); + mesh, vertices_t0, vertices_t1, + /*inflation_radius=*/0.5 * min_distance.value_or(0.0), broad_phase); // Narrow phase return candidates.compute_collision_free_stepsize( diff --git a/src/ipc/ipc.hpp b/src/ipc/ipc.hpp index 6b7d0ccfb..c57cebbaa 100644 --- a/src/ipc/ipc.hpp +++ b/src/ipc/ipc.hpp @@ -7,6 +7,8 @@ #include #include +#include + /// Incremental Potential Contact functions namespace ipc { @@ -18,7 +20,7 @@ namespace ipc { /// @param mesh The collision mesh. /// @param vertices_t0 Surface vertex vertices at start as rows of a matrix. /// @param vertices_t1 Surface vertex vertices at end as rows of a matrix. -/// @param min_distance The minimum distance allowable between any two elements. +/// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @param broad_phase The broad phase method to use. /// @param narrow_phase_ccd The narrow phase CCD algorithm to use. /// @returns True if any collisions occur. @@ -26,7 +28,7 @@ bool is_step_collision_free( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr, const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD); @@ -38,7 +40,7 @@ bool is_step_collision_free( /// @param mesh The collision mesh. /// @param vertices_t0 Vertex vertices at start as rows of a matrix. Assumes vertices_t0 is intersection free. /// @param vertices_t1 Surface vertex vertices at end as rows of a matrix. -/// @param min_distance The minimum distance allowable between any two elements. +/// @param min_distance Deprecated: give the mesh a thickness (CollisionMesh::set_thickness()) instead. If given, every pair keeps this distance and the mesh's thickness is ignored; otherwise each pair keeps its two primitives' thicknesses (zero without a thickness). /// @param broad_phase The broad phase method to use. /// @param narrow_phase_ccd The narrow phase CCD algorithm to use. /// @returns A step-size \f$\in [0, 1]\f$ that is collision free. A value of 1.0 if a full step and 0.0 is no step. @@ -46,7 +48,7 @@ double compute_collision_free_stepsize( const CollisionMesh& mesh, Eigen::ConstRef vertices_t0, Eigen::ConstRef vertices_t1, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr, const NarrowPhaseCCD& narrow_phase_ccd = DEFAULT_NARROW_PHASE_CCD); diff --git a/src/ipc/ogc/trust_region.cpp b/src/ipc/ogc/trust_region.cpp index 3caeee7bf..ee444e181 100644 --- a/src/ipc/ogc/trust_region.cpp +++ b/src/ipc/ogc/trust_region.cpp @@ -17,7 +17,7 @@ void TrustRegion::warm_start_time_step( Eigen::ConstRef pred_x, NormalCollisions& collisions, const double dhat, - const double min_distance, + const std::optional min_distance, BroadPhase* broad_phase) { // Assign a default broad phase if none is provided. @@ -67,7 +67,7 @@ void TrustRegion::update( const CollisionMesh& mesh, Eigen::ConstRef x, NormalCollisions& collisions, - const double min_distance, + const std::optional min_distance, BroadPhase* broad_phase) { assert(0 < relaxed_radius_scaling && relaxed_radius_scaling < 1); @@ -83,7 +83,8 @@ void TrustRegion::update( } candidates.build( - mesh, vertices, trust_region_inflation_radius + min_distance, + mesh, vertices, + trust_region_inflation_radius + min_distance.value_or(0.0), broad_phase); // Compute the trust region distances for the reduced set of collision diff --git a/src/ipc/ogc/trust_region.hpp b/src/ipc/ogc/trust_region.hpp index f676bf6c2..6b69bce64 100644 --- a/src/ipc/ogc/trust_region.hpp +++ b/src/ipc/ogc/trust_region.hpp @@ -5,6 +5,8 @@ #include #include +#include + namespace ipc::ogc { /// @brief A trust region for filtering optimization steps. @@ -51,7 +53,7 @@ struct TrustRegion { /// @param pred_x Predicted vertex positions. /// @param collisions Collisions to be initialized. /// @param dhat The offset distance for contact. - /// @param min_distance Minimum distance between elements. + /// @param min_distance Deprecated minimum distance between elements, which overrides the mesh's thickness when given (see Candidates::compute_per_vertex_safe_distances()). /// @param broad_phase Broad phase collision detection. void warm_start_time_step( const CollisionMesh& mesh, @@ -59,33 +61,33 @@ struct TrustRegion { Eigen::ConstRef pred_x, NormalCollisions& collisions, const double dhat, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr); /// @brief Update the trust region based on the current positions. /// @param[in] mesh The collision mesh. /// @param[in] x Current vertex positions. /// @param[out] collisions Collisions to be updated. - /// @param[in] min_distance Minimum distance between elements. + /// @param[in] min_distance Deprecated minimum distance between elements, which overrides the mesh's thickness when given (see Candidates::compute_per_vertex_safe_distances()). /// @param[in] broad_phase Broad phase collision detection. void update( const CollisionMesh& mesh, Eigen::ConstRef x, NormalCollisions& collisions, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr); /// @brief Update the trust region if should_update_trust_region is true. /// @param[in] mesh The collision mesh. /// @param[in] x Current vertex positions. /// @param[out] collisions Collisions to be updated. - /// @param[in] min_distance Minimum distance between elements. + /// @param[in] min_distance Deprecated minimum distance between elements, which overrides the mesh's thickness when given (see Candidates::compute_per_vertex_safe_distances()). /// @param[in] broad_phase Broad phase collision detection. void update_if_needed( const CollisionMesh& mesh, Eigen::ConstRef x, NormalCollisions& collisions, - const double min_distance = 0.0, + const std::optional min_distance = std::nullopt, BroadPhase* broad_phase = nullptr) { if (should_update_trust_region) { diff --git a/src/ipc/potentials/tangential_potential.cpp b/src/ipc/potentials/tangential_potential.cpp index 965113689..db6feb1f4 100644 --- a/src/ipc/potentials/tangential_potential.cpp +++ b/src/ipc/potentials/tangential_potential.cpp @@ -1,5 +1,6 @@ #include "friction_potential.hpp" +#include #include #include #include @@ -45,7 +46,7 @@ Eigen::VectorXd TangentialPotential::force( Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, const NormalPotential& normal_potential, - const double dmin, + const std::optional dmin, const bool no_mu) const { if (collisions.empty()) { @@ -77,7 +78,7 @@ Eigen::SparseMatrix TangentialPotential::force_jacobian( Eigen::ConstRef velocities, const NormalPotential& normal_potential, const DiffWRT wrt, - const double dmin) const + const std::optional dmin) const { if (collisions.empty()) { return Eigen::SparseMatrix( @@ -329,9 +330,12 @@ VectorMax12d TangentialPotential::force( Eigen::ConstRef lagged_displacements, // = u Eigen::ConstRef velocities, // = v const NormalPotential& normal_potential, - const double dmin, + const std::optional dmin, const bool no_mu) const { + if (dmin.has_value()) { + warn_dmin_deprecated(); + } // x is the rest position // u is the displacement at the beginning of the lagged solve // v is the current velocity @@ -359,7 +363,8 @@ VectorMax12d TangentialPotential::force( // Compute N(x + u) const double N = normal_potential.force_magnitude( - collision.compute_distance(lagged_positions), dmin); + collision.compute_distance(lagged_positions), + dmin.value_or(collision.dmin)); // Compute P const MatrixMax P = @@ -408,8 +413,11 @@ MatrixMax12d TangentialPotential::force_jacobian( Eigen::ConstRef velocities, // = v const NormalPotential& normal_potential, const DiffWRT wrt, - const double dmin) const + const std::optional dmin) const { + if (dmin.has_value()) { + warn_dmin_deprecated(); + } // x is the rest position // u is the displacement at the beginning of the lagged solve // v is the current velocity @@ -438,7 +446,8 @@ MatrixMax12d TangentialPotential::force_jacobian( // Compute N const double N = normal_potential.force_magnitude( - collision.compute_distance(lagged_positions), dmin); + collision.compute_distance(lagged_positions), + dmin.value_or(collision.dmin)); // Compute ∇N VectorMax12d grad_N; @@ -446,7 +455,8 @@ MatrixMax12d TangentialPotential::force_jacobian( // ∇ₓN = ∇ᵤN grad_N = normal_potential.force_magnitude_gradient( collision.compute_distance(lagged_positions), - collision.compute_distance_gradient(lagged_positions), dmin); + collision.compute_distance_gradient(lagged_positions), + dmin.value_or(collision.dmin)); assert(grad_N.array().isFinite().all()); } @@ -590,9 +600,12 @@ Eigen::VectorXd TangentialPotential::gcp_force( Eigen::ConstRef rest_positions, Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, - const double dmin, + const std::optional dmin, const bool no_mu) const { + if (dmin.has_value()) { + warn_dmin_deprecated(); + } if (collisions.empty()) { return Eigen::VectorXd::Zero(velocities.size()); } @@ -621,9 +634,12 @@ Eigen::SparseMatrix TangentialPotential::gcp_force_jacobian( Eigen::ConstRef velocities, const GCPParameters& params, const DiffWRT wrt, - const double dmin, + const std::optional dmin, const bool no_mu) const { + if (dmin.has_value()) { + warn_dmin_deprecated(); + } if (collisions.empty()) { return Eigen::SparseMatrix( velocities.size(), velocities.size()); diff --git a/src/ipc/potentials/tangential_potential.hpp b/src/ipc/potentials/tangential_potential.hpp index 6f917583f..c35deeaa3 100644 --- a/src/ipc/potentials/tangential_potential.hpp +++ b/src/ipc/potentials/tangential_potential.hpp @@ -3,6 +3,8 @@ #include #include +#include + namespace ipc { /// @brief A tangential dissipative potential. @@ -37,7 +39,7 @@ class TangentialPotential : public Potential { /// @param lagged_displacements Previous displacements of the vertices (rowwise). /// @param velocities Current displacements of the vertices (rowwise). /// @param normal_potential Normal potential (used for normal force magnitude). - /// @param dmin Minimum distance (used for normal force magnitude). + /// @param dmin Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin (TangentialCollision::dmin, its primitives' thicknesses). /// @param no_mu whether to not multiply by mu /// @return The friction force. Eigen::VectorXd force( @@ -47,7 +49,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, const NormalPotential& normal_potential, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const; /// @brief Compute the Jacobian of the friction force wrt the velocities. @@ -58,7 +60,7 @@ class TangentialPotential : public Potential { /// @param velocities Current displacements of the vertices (rowwise). /// @param normal_potential Normal potential (used for normal force magnitude). /// @param wrt The variable to take the derivative with respect to. - /// @param dmin Minimum distance (used for normal force magnitude). + /// @param dmin Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin (TangentialCollision::dmin, its primitives' thicknesses). /// @return The Jacobian of the friction force wrt the velocities. Eigen::SparseMatrix force_jacobian( const TangentialCollisions& collisions, @@ -68,7 +70,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef velocities, const NormalPotential& normal_potential, const DiffWRT wrt, - const double dmin = 0) const; + const std::optional dmin = std::nullopt) const; Eigen::VectorXd gcp_force( const TangentialCollisions& collisions, @@ -76,7 +78,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef rest_positions, Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const; Eigen::SparseMatrix gcp_force_jacobian( @@ -87,7 +89,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef velocities, const GCPParameters& params, const DiffWRT wrt, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const; // -- Single collision methods --------------------------------------------- @@ -125,7 +127,7 @@ class TangentialPotential : public Potential { /// @param lagged_displacements Previous displacements of the vertices (rowwise). /// @param velocities Current displacements of the vertices (rowwise). /// @param normal_potential Normal potential (used for normal force magnitude). - /// @param dmin Minimum distance (used for normal force magnitude). + /// @param dmin Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin (TangentialCollision::dmin, its primitives' thicknesses). /// @param no_mu Whether to not multiply by mu /// @return Friction force VectorMax12d force( @@ -134,7 +136,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, const NormalPotential& normal_potential, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const; //< whether to not multiply by mu /// @brief Compute the friction force Jacobian. @@ -144,7 +146,7 @@ class TangentialPotential : public Potential { /// @param velocities Current displacements of the vertices (rowwise). /// @param normal_potential Normal potential (used for normal force magnitude). /// @param wrt Variable to differentiate the friction force with respect to. - /// @param dmin Minimum distance (used for normal force magnitude). + /// @param dmin Deprecated minimum distance (used for normal force magnitude). If given, it overrides each collision's own dmin (TangentialCollision::dmin, its primitives' thicknesses). /// @return Friction force Jacobian MatrixMax12d force_jacobian( const TangentialCollision& collision, @@ -153,7 +155,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef velocities, const NormalPotential& normal_potential, const DiffWRT wrt, - const double dmin = 0) const; + const std::optional dmin = std::nullopt) const; VectorMaxNd gcp_force( const TangentialCollision& collision, @@ -207,7 +209,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef rest_positions, Eigen::ConstRef lagged_displacements, Eigen::ConstRef velocities, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const { return gcp_force( @@ -226,7 +228,7 @@ class TangentialPotential : public Potential { Eigen::ConstRef velocities, const GCPParameters& params, const DiffWRT wrt, - const double dmin = 0, + const std::optional dmin = std::nullopt, const bool no_mu = false) const { return gcp_force_jacobian( diff --git a/tests/src/tests/collisions/CMakeLists.txt b/tests/src/tests/collisions/CMakeLists.txt index bf61942ef..7a5afeb30 100644 --- a/tests/src/tests/collisions/CMakeLists.txt +++ b/tests/src/tests/collisions/CMakeLists.txt @@ -3,6 +3,7 @@ set(SOURCES test_normal_collisions.cpp test_obstacles.cpp test_plane_vertex_collisions.cpp + test_thickness.cpp # Benchmarks diff --git a/tests/src/tests/collisions/test_thickness.cpp b/tests/src/tests/collisions/test_thickness.cpp new file mode 100644 index 000000000..c9b92d251 --- /dev/null +++ b/tests/src/tests/collisions/test_thickness.cpp @@ -0,0 +1,388 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +using namespace ipc; + +namespace { + +/// Two triangles facing each other, the upper one gap above the lower and +/// over its interior. +void two_triangles( + const double gap, + Eigen::MatrixXd& V, + Eigen::MatrixXi& E, + Eigen::MatrixXi& F) +{ + V.resize(6, 3); + V << 0, 0, 0, // + 2, 0, 0, // + 0, 2, 0, // + 0.2, 0.2, gap, 1.2, 0.2, gap, 0.2, 1.2, gap; + F.resize(2, 3); + F << 0, 1, 2, 3, 4, 5; + igl::edges(F, E); +} + +/// Every collision keeps its own two primitives' thicknesses apart. +void check_dmin(const NormalCollisions& collisions, const CollisionMesh& mesh) +{ + for (const auto& c : collisions.vv_collisions) { + CHECK(c.dmin == Catch::Approx(mesh.pair_thickness(c))); + } + for (const auto& c : collisions.ev_collisions) { + CHECK(c.dmin == Catch::Approx(mesh.pair_thickness(c))); + } + for (const auto& c : collisions.ee_collisions) { + CHECK(c.dmin == Catch::Approx(mesh.pair_thickness(c))); + } + for (const auto& c : collisions.fv_collisions) { + CHECK(c.dmin == Catch::Approx(mesh.pair_thickness(c))); + } + for (const auto& c : collisions.pv_collisions) { + CHECK(c.dmin == Catch::Approx(mesh.pair_thickness(c))); + } +} + +} // namespace + +TEST_CASE( + "A uniform thickness matches a uniform dmin", "[collisions][thickness]") +{ + const auto broad_phase = GENERATE(tests::BroadPhaseGenerator::create()); +#ifndef NDEBUG + if (broad_phase->name() == "BruteForce") { + SKIP("Brute force is too slow in debug mode"); + } +#endif + const NormalCollisions::CollisionSetType collision_set_type = GENERATE( + NormalCollisions::CollisionSetType::IPC, + NormalCollisions::CollisionSetType::IMPROVED_MAX_APPROX, + NormalCollisions::CollisionSetType::OGC); + CAPTURE(broad_phase->name(), collision_set_type); + + Eigen::MatrixXd V; + Eigen::MatrixXi E, F; + REQUIRE(tests::load_mesh("two-cubes-close.ply", V, E, F)); + + const double dhat = 1e-1; + const double thickness = 1e-2; + + // Every pair of primitives is 2 * thickness thick, so contacts rest + // exactly where a dmin of 2 * thickness puts them. + const CollisionMesh plain_mesh(V, E, F); + NormalCollisions expected; + expected.set_collision_set_type(collision_set_type); + expected.build(plain_mesh, V, dhat, 2 * thickness, broad_phase.get()); + REQUIRE(!expected.empty()); + + CollisionMesh mesh(V, E, F); + mesh.set_vertex_thickness( + Eigen::VectorXd::Constant(mesh.num_vertices(), thickness)); + NormalCollisions collisions; + collisions.set_collision_set_type(collision_set_type); + collisions.build(mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + + CHECK(collisions.vv_collisions.size() == expected.vv_collisions.size()); + CHECK(collisions.ev_collisions.size() == expected.ev_collisions.size()); + CHECK(collisions.ee_collisions.size() == expected.ee_collisions.size()); + CHECK(collisions.fv_collisions.size() == expected.fv_collisions.size()); + check_dmin(collisions, mesh); + + const BarrierPotential barrier(dhat, /*stiffness=*/1.0); + // The improved max approximation's sum of positive and negative terms + // already differs between two builds of the same set (~1e-5 in energy and + // ~1e-3 in gradient, as its merge order is not deterministic), so it is + // compared more loosely. + const double tol = collision_set_type + == NormalCollisions::CollisionSetType::IMPROVED_MAX_APPROX + ? 1e-2 + : 1e-8; + CHECK( + barrier(collisions, mesh, V) + == Catch::Approx(barrier(expected, plain_mesh, V)).epsilon(tol)); + CHECK(barrier.gradient(collisions, mesh, V) + .isApprox(barrier.gradient(expected, plain_mesh, V), tol)); +} + +TEST_CASE( + "Contacts rest at the sum of their thicknesses", "[collisions][thickness]") +{ + const auto broad_phase = GENERATE(tests::BroadPhaseGenerator::create()); + CAPTURE(broad_phase->name()); + + const double dhat = 0.25; + const double lower = 0.2, upper = 0.1; // face thicknesses + + Eigen::MatrixXd V; + Eigen::MatrixXi E, F; + + SECTION("Within dhat of the thicknesses") + { + two_triangles(/*gap=*/0.5, V, E, F); + + const CollisionMesh plain_mesh(V, E, F); + NormalCollisions plain; + plain.build( + plain_mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + CHECK(plain.empty()); // 0.5 > dhat + + CollisionMesh mesh(V, E, F); + mesh.set_face_thickness(Eigen::Vector2d(lower, upper)); + NormalCollisions collisions; + collisions.build( + mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + REQUIRE(!collisions.fv_collisions.empty()); + check_dmin(collisions, mesh); + for (const auto& fv : collisions.fv_collisions) { + CHECK(fv.dmin == Catch::Approx(lower + upper)); + } + + const BarrierPotential barrier(dhat, /*stiffness=*/1.0); + CHECK(barrier(collisions, mesh, V) > 0); + + // A dmin (deprecated) overrides the thickness. + collisions.build(mesh, V, dhat, /*dmin=*/0.0, broad_phase.get()); + CHECK(collisions.empty()); + collisions.build(mesh, V, dhat, /*dmin=*/0.4, broad_phase.get()); + REQUIRE(!collisions.empty()); + for (size_t i = 0; i < collisions.size(); i++) { + CHECK(collisions[i].dmin == 0.4); + } + } + + SECTION("Beyond dhat of the thicknesses") + { + two_triangles(/*gap=*/lower + upper + dhat + 1e-3, V, E, F); + + CollisionMesh mesh(V, E, F); + mesh.set_face_thickness(Eigen::Vector2d(lower, upper)); + NormalCollisions collisions; + collisions.build( + mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + CHECK(collisions.empty()); + } + + SECTION("Each primitive type keeps its own thickness") + { + two_triangles(/*gap=*/0.5, V, E, F); + + CollisionMesh mesh(V, E, F); + Eigen::VectorXd vertex(6), edge(E.rows()), face(2); + vertex << 0.12, 0.12, 0.12, 0.1, 0.1, 0.1; + for (int i = 0; i < E.rows(); i++) { + edge[i] = E(i, 0) < 3 ? 0.15 : 0.1; + } + face << 0.2, 0.1; + mesh.set_thickness(vertex, edge, face); + + NormalCollisions collisions; + collisions.build( + mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + REQUIRE(!collisions.fv_collisions.empty()); + check_dmin(collisions, mesh); + } +} + +TEST_CASE("Codim. vertices with a thickness", "[collisions][thickness][codim]") +{ + const auto broad_phase = GENERATE(tests::BroadPhaseGenerator::create()); + CAPTURE(broad_phase->name()); + + Eigen::MatrixXd V(2, 3); + V << 0, 0, 0, 1.2, 0, 0; + CollisionMesh mesh(V); + REQUIRE(mesh.num_codim_vertices() == 2); + + const double dhat = 0.3; + NormalCollisions collisions; + collisions.build(mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + CHECK(collisions.empty()); + + mesh.set_vertex_thickness(Eigen::Vector2d(0.5, 0.5)); + collisions.build(mesh, V, dhat, /*dmin=*/std::nullopt, broad_phase.get()); + REQUIRE(collisions.vv_collisions.size() == 1); + CHECK(collisions.vv_collisions[0].dmin == Catch::Approx(1.0)); +} + +TEST_CASE("Plane-vertex with a thickness", "[collisions][thickness][plane]") +{ + Eigen::MatrixXd V(1, 3); + V << 0, 0.25, 0; + CollisionMesh mesh(V); + mesh.planes.emplace_back(Eigen::Vector3d::UnitY(), 0); + + const double dhat = 0.15; + NormalCollisions collisions; + collisions.build(mesh, V, dhat); + CHECK(collisions.empty()); + + mesh.set_vertex_thickness(Eigen::VectorXd::Constant(1, 0.2)); + collisions.build(mesh, V, dhat); + REQUIRE(collisions.pv_collisions.size() == 1); + CHECK(collisions.pv_collisions[0].dmin == Catch::Approx(0.2)); + + V(0, 1) = 0.4; + collisions.build(mesh, V, dhat); + CHECK(collisions.empty()); +} + +TEST_CASE("CCD stops at the thickness", "[ccd][thickness]") +{ +#if (defined(WIN32) || defined(_WIN32) || defined(__WIN32)) && !defined(NDEBUG) + SKIP("'CCD stops at the thickness' test is skipped in debug mode"); +#endif + + const auto broad_phase = GENERATE(tests::BroadPhaseGenerator::create()); + CAPTURE(broad_phase->name()); + + Eigen::MatrixXd V0, V1; + Eigen::MatrixXi E, F; + two_triangles(/*gap=*/1, V0, E, F); + V1 = V0; + V1.bottomRows(3).col(2).setConstant(-1); // the upper one falls through + + CollisionMesh mesh(V0, E, F); + const double thickness = 0.1; // per face, so a pair rests 0.2 apart + + const double plain_step = compute_collision_free_stepsize( + mesh, V0, V1, /*min_distance=*/std::nullopt, broad_phase.get()); + CHECK(plain_step <= 0.5); + CHECK(plain_step > 0.3); + + mesh.set_face_thickness(Eigen::Vector2d::Constant(thickness)); + const double step = compute_collision_free_stepsize( + mesh, V0, V1, /*min_distance=*/std::nullopt, broad_phase.get()); + // Contact at a gap of 0.2, so z = 0.2 at t = 0.4. + CHECK(step <= 0.4); + CHECK(step > 0.3); + CHECK(step < plain_step); + + // A min_distance (deprecated) overrides the thickness. + CHECK( + compute_collision_free_stepsize( + mesh, V0, V1, /*min_distance=*/0.0, broad_phase.get()) + == Catch::Approx(plain_step)); + + Eigen::MatrixXd V_close = V0; + V_close.bottomRows(3).col(2).setConstant(0.15); + CHECK(!is_step_collision_free( + mesh, V0, V_close, /*min_distance=*/std::nullopt, broad_phase.get())); + CHECK(is_step_collision_free( + mesh, V0, V_close, /*min_distance=*/0.0, broad_phase.get())); + mesh.clear_thickness(); + CHECK(is_step_collision_free( + mesh, V0, V_close, /*min_distance=*/std::nullopt, broad_phase.get())); +} + +TEST_CASE( + "Per-vertex safe distances with a thickness", "[candidates][thickness]") +{ + Eigen::MatrixXd V; + Eigen::MatrixXi E, F; + two_triangles(/*gap=*/0.5, V, E, F); + + CollisionMesh mesh(V, E, F); + mesh.set_face_thickness(Eigen::Vector2d(0.2, 0.1)); + + const double inflation_radius = 0.5; + Candidates candidates; + candidates.build(mesh, V, inflation_radius); + REQUIRE(!candidates.empty()); + + const Eigen::VectorXd safe = + candidates.compute_per_vertex_safe_distances(mesh, V, inflation_radius); + // Vertex 3 is 0.5 above the lower face's interior and the pair is 0.3 + // thick. + CHECK(safe[3] == Catch::Approx(0.2)); + CHECK(safe.minCoeff() >= 0); +} + +TEST_CASE("Friction keeps each collision's thickness", "[friction][thickness]") +{ + Eigen::MatrixXd V; + Eigen::MatrixXi E, F; + two_triangles(/*gap=*/0.5, V, E, F); + + CollisionMesh mesh(V, E, F); + mesh.set_face_thickness(Eigen::Vector2d(0.2, 0.1)); + + const double dhat = 0.25; + const BarrierPotential barrier(dhat, /*stiffness=*/1.0); + NormalCollisions collisions; + collisions.build(mesh, V, dhat); + REQUIRE(!collisions.empty()); + + TangentialCollisions tangential_collisions; + tangential_collisions.build(mesh, V, collisions, barrier, /*mu=*/0.5); + // Mollified edge-edge pairs have no friction, so there may be fewer. + // Every pair here is 0.2 + 0.1 thick. + REQUIRE(!tangential_collisions.empty()); + for (size_t i = 0; i < tangential_collisions.size(); i++) { + CHECK(tangential_collisions[i].dmin == Catch::Approx(0.3)); + CHECK(tangential_collisions[i].normal_force_magnitude > 0); + } + + // Without a dmin, the friction force uses each collision's own, so it + // matches passing the pairs' (uniform) thickness explicitly. + const FrictionPotential friction(/*eps_v=*/1e-3); + const Eigen::MatrixXd U = Eigen::MatrixXd::Zero(V.rows(), V.cols()); + Eigen::MatrixXd velocities = U; + velocities.col(0).setConstant(1e-4); + const Eigen::VectorXd force = + friction.force(tangential_collisions, mesh, V, U, velocities, barrier); + CHECK(force.squaredNorm() > 0); + CHECK(force.isApprox(friction.force( + tangential_collisions, mesh, V, U, velocities, barrier, + /*dmin=*/0.3))); + CHECK(!force.isApprox(friction.force( + tangential_collisions, mesh, V, U, velocities, barrier, + /*dmin=*/0.0))); +} + +TEST_CASE( + "Semi-implicit stiffness keeps each stencil's thickness", + "[barrier][thickness]") +{ + Eigen::MatrixXd V; + Eigen::MatrixXi E, F; + two_triangles(/*gap=*/0.5, V, E, F); + + CollisionMesh mesh(V, E, F); + mesh.set_face_thickness(Eigen::Vector2d(0.2, 0.1)); + + NormalCollisions collisions; + collisions.build(mesh, V, /*dhat=*/0.25); + REQUIRE(!collisions.empty()); + Candidates candidates; + candidates.build(mesh, V, /*inflation_radius=*/0.125); + REQUIRE(!candidates.empty()); + + const Eigen::VectorXd masses = Eigen::VectorXd::Ones(V.rows()); + Eigen::SparseMatrix hess(V.size(), V.size()); + hess.setIdentity(); + + // Every pair here is 0.2 + 0.1 thick, so the default matches passing it. + CHECK(semi_implicit_stiffness(mesh, V, collisions, masses, hess) + .isApprox(semi_implicit_stiffness( + mesh, V, collisions, masses, hess, /*dmin=*/0.3))); + CHECK(semi_implicit_stiffness(mesh, V, candidates, masses, hess) + .isApprox(semi_implicit_stiffness( + mesh, V, candidates, masses, hess, /*dmin=*/0.3))); + CHECK(!semi_implicit_stiffness(mesh, V, candidates, masses, hess) + .isApprox(semi_implicit_stiffness( + mesh, V, candidates, masses, hess, /*dmin=*/0.0))); +} diff --git a/tests/src/tests/potential/test_esp_potential.cpp b/tests/src/tests/potential/test_esp_potential.cpp index 46124e9ac..1c3102e67 100644 --- a/tests/src/tests/potential/test_esp_potential.cpp +++ b/tests/src/tests/potential/test_esp_potential.cpp @@ -1950,6 +1950,10 @@ TEST_CASE( "ESP potential 2D finite differences", "[esp_potential], [esp_potential_2d]") { +#if (defined(WIN32) || defined(_WIN32) || defined(__WIN32)) && !defined(NDEBUG) + SKIP("'ESP potential 2D finite differences' test is skipped in debug mode"); +#endif + const int quadrature_order = GENERATE(1, 14); Eigen::MatrixXd V; diff --git a/tests/src/tests/test_collision_mesh.cpp b/tests/src/tests/test_collision_mesh.cpp index eeff97f66..1086d9fa2 100644 --- a/tests/src/tests/test_collision_mesh.cpp +++ b/tests/src/tests/test_collision_mesh.cpp @@ -3,6 +3,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -334,3 +338,99 @@ TEST_CASE("User-provided collision mesh areas", "[collision_mesh][areas]") } } } + +TEST_CASE("Collision mesh thickness", "[collision_mesh][thickness]") +{ + // Two triangles sharing the edge 1-2, as two panels sewn at a seam, plus + // a lone edge 4-5 on no face. + Eigen::MatrixXd V(6, 3); + V << 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 3, 0, 0, 4, 0, 0; + Eigen::MatrixXi F(2, 3); + F << 0, 1, 2, 1, 3, 2; + Eigen::MatrixXi E(6, 2); + E << 0, 1, 1, 2, 2, 0, 1, 3, 3, 2, 4, 5; + + CollisionMesh mesh(V, E, F); + CHECK(!mesh.has_thickness()); + CHECK(mesh.vertex_thickness(0) == 0); + CHECK(mesh.max_thickness() == 0); + + const auto edge_id = [&](const index_t a, const index_t b) { + for (int i = 0; i < mesh.num_edges(); i++) { + const auto& e = mesh.edges().row(i); + if ((e(0) == a && e(1) == b) || (e(0) == b && e(1) == a)) { + return index_t(i); + } + } + return index_t(-1); + }; + + SECTION("From faces") + { + Eigen::VectorXd vertex = Eigen::VectorXd::Zero(6); + vertex[5] = 0.7; + mesh.set_face_thickness(Eigen::Vector2d(1.5, 2.5), vertex); + REQUIRE(mesh.has_thickness()); + + CHECK(mesh.face_thickness(0) == 1.5); + CHECK(mesh.face_thickness(1) == 2.5); + // Vertex 0 is on the thin face only, vertex 3 on the thick one, and + // the seam vertices on both. + CHECK(mesh.vertex_thickness(0) == 1.5); + CHECK(mesh.vertex_thickness(3) == 2.5); + CHECK(mesh.vertex_thickness(1) == 2.5); + CHECK(mesh.vertex_thickness(2) == 2.5); + CHECK(mesh.edge_thickness(edge_id(0, 1)) == 1.5); + CHECK(mesh.edge_thickness(edge_id(1, 3)) == 2.5); + CHECK(mesh.edge_thickness(edge_id(1, 2)) == 2.5); // the seam + // On no face: the vertices keep their own, the edge the larger. + CHECK(mesh.vertex_thickness(4) == 0); + CHECK(mesh.vertex_thickness(5) == 0.7); + CHECK(mesh.edge_thickness(edge_id(4, 5)) == 0.7); + CHECK(mesh.max_thickness() == 2.5); + + // A pair rests at the sum of its two sides. + CHECK(mesh.pair_thickness(FaceVertexCandidate(1, 0)) == 2.5 + 1.5); + CHECK( + mesh.pair_thickness(EdgeVertexCandidate(edge_id(4, 5), 0)) + == 0.7 + 1.5); + CHECK( + mesh.pair_thickness(EdgeEdgeCandidate(edge_id(0, 1), edge_id(4, 5))) + == 1.5 + 0.7); + CHECK(mesh.pair_thickness(VertexVertexCandidate(0, 4)) == 1.5); + CHECK( + FaceVertexCandidate(1, 0).thickness( + mesh.vertex_thicknesses(), mesh.edge_thicknesses(), + mesh.face_thicknesses()) + == 2.5 + 1.5); + + mesh.clear_thickness(); + CHECK(!mesh.has_thickness()); + CHECK(mesh.pair_thickness(FaceVertexCandidate(1, 0)) == 0); + } + + SECTION("From vertices") + { + Eigen::VectorXd vertex(6); + vertex << 1, 2, 3, 4, 5, 6; + mesh.set_vertex_thickness(vertex); + CHECK(mesh.vertex_thicknesses() == vertex); + CHECK(mesh.edge_thickness(edge_id(0, 1)) == 2); + CHECK(mesh.edge_thickness(edge_id(4, 5)) == 6); + CHECK(mesh.face_thickness(0) == 3); + CHECK(mesh.face_thickness(1) == 4); + CHECK(mesh.max_thickness() == 6); + } + + SECTION("Invalid thicknesses are rejected") + { + CHECK_THROWS(mesh.set_vertex_thickness(Eigen::VectorXd::Ones(3))); + CHECK_THROWS( + mesh.set_vertex_thickness(Eigen::VectorXd::Constant(6, -1))); + CHECK_THROWS(mesh.set_face_thickness(Eigen::VectorXd::Ones(3))); + CHECK_THROWS(mesh.set_thickness( + Eigen::VectorXd::Ones(6), Eigen::VectorXd::Ones(6), + Eigen::VectorXd::Constant(2, std::nan("")))); + CHECK(!mesh.has_thickness()); + } +}