Fix parallel edge-edge dtype

src/ipc/candidates/candidates.cpp

@@ -6,7 +6,6 @@
 
 #include <tbb/parallel_for.h>
 #include <tbb/blocked_range.h>
-#include <tbb/enumerable_thread_specific.h>
 #include <shared_mutex>
 
 #include <fstream>

src/ipc/collisions/collision_constraints.cpp

@@ -153,11 +153,7 @@ double CollisionConstraints::compute_potential(
  }
  });
 
- double potential = 0;
- for (const auto& local_potential : storage) {
- potential += local_potential;
- }
- return potential;
+ return storage.combine([](double a, double b) { return a + b; });
 }
 
 Eigen::VectorXd CollisionConstraints::compute_potential_gradient(
@@ -191,11 +187,8 @@ Eigen::VectorXd CollisionConstraints::compute_potential_gradient(
  }
  });
 
- Eigen::VectorXd grad = Eigen::VectorXd::Zero(vertices.size());
- for (const auto& local_grad : storage) {
- grad += local_grad;
- }
- return grad;
+ return storage.combine([](const Eigen::VectorXd& a,
+ const Eigen::VectorXd& b) { return a + b; });
 }
 
 Eigen::SparseMatrix<double> CollisionConstraints::compute_potential_hessian(
@@ -347,11 +340,7 @@ double CollisionConstraints::compute_minimum_distance(
  }
  });
 
- double min_dist = std::numeric_limits<double>::infinity();
- for (const auto& local_min_dist : storage) {
- min_dist = std::min(min_dist, local_min_dist);
- }
- return min_dist;
+ return storage.combine([](double a, double b) { return std::min(a, b); });
 }
 
 // ============================================================================

src/ipc/collisions/collision_constraints.hpp

@@ -12,8 +12,6 @@
 
 #include <Eigen/Core>
 
-#include <tbb/enumerable_thread_specific.h>
-
 #include <vector>
 
 namespace ipc {

src/ipc/distance/distance_type.cpp

@@ -91,46 +91,61 @@ EdgeEdgeDistanceType edge_edge_distance_type(
  const Eigen::Vector3d v = eb1 - eb0;
  const Eigen::Vector3d w = ea0 - eb0;
 
- const double a = u.squaredNorm(); // always >= 0
+ const double a = u.squaredNorm(); // always ≥ 0
  const double b = u.dot(v);
- const double c = v.squaredNorm(); // always >= 0
+ const double c = v.squaredNorm(); // always ≥ 0
  const double d = u.dot(w);
  const double e = v.dot(w);
- const double D = a * c - b * b; // always >= 0
+ const double D = a * c - b * b; // always ≥ 0
+
+ // Degenerate cases should not happen in practice, but we handle them
+ if (a == 0.0 && c == 0.0) {
+ return EdgeEdgeDistanceType::EA0_EB0;
+ } else if (a == 0.0) {
+ return EdgeEdgeDistanceType::EA0_EB;
+ } else if (c == 0.0) {
+ return EdgeEdgeDistanceType::EA_EB0;
+ }
+
+ // Special handling for parallel edges
+ const double parallel_tolerance = PARALLEL_THRESHOLD * std::max(1.0, a * c);
+ if (u.cross(v).squaredNorm() < parallel_tolerance) {
+ return edge_edge_parallel_distance_type(ea0, ea1, eb0, eb1);
+ }
 
- EdgeEdgeDistanceType defaultCase = EdgeEdgeDistanceType::EA_EB;
+ EdgeEdgeDistanceType default_case = EdgeEdgeDistanceType::EA_EB;
 
  // compute the line parameters of the two closest points
  const double sN = (b * e - c * d);
  double tN, tD; // tc = tN / tD
- if (sN <= 0.0) { // sc < 0 => the s=0 edge is visible
+ if (sN <= 0.0) { // sc < 0 ⟹ the s=0 edge is visible
  tN = e;
  tD = c;
- defaultCase = EdgeEdgeDistanceType::EA0_EB;
- } else if (sN >= D) { // sc > 1  => the s=1 edge is visible
+ default_case = EdgeEdgeDistanceType::EA0_EB;
+ } else if (sN >= D) { // sc > 1 ⟹ the s=1 edge is visible
  tN = e + b;
  tD = c;
- defaultCase = EdgeEdgeDistanceType::EA1_EB;
+ default_case = EdgeEdgeDistanceType::EA1_EB;
  } else {
  tN = (a * e - b * d);
- tD = D; // default tD = D >= 0
+ tD = D; // default tD = D ≥ 0
  if (tN > 0.0 && tN < tD
- && u.cross(v).squaredNorm() < PARALLEL_THRESHOLD * a * c) {
- // avoid nearly parallel EE
+ && u.cross(v).squaredNorm() < parallel_tolerance) {
+ // avoid coplanar or nearly parallel EE
  if (sN < D / 2) {
  tN = e;
  tD = c;
- defaultCase = EdgeEdgeDistanceType::EA0_EB;
+ default_case = EdgeEdgeDistanceType::EA0_EB;
  } else {
  tN = e + b;
  tD = c;
- defaultCase = EdgeEdgeDistanceType::EA1_EB;
+ default_case = EdgeEdgeDistanceType::EA1_EB;
  }
  }
- // else defaultCase stays EdgeEdgeDistanceType::EA_EB
+ // else default_case stays EdgeEdgeDistanceType::EA_EB
  }
 
- if (tN <= 0.0) { // tc < 0 => the t=0 edge is visible
+ if (tN <= 0.0) { // tc < 0 ⟹ the t=0 edge is visible
  // recompute sc for this edge
  if (-d <= 0.0) {
  return EdgeEdgeDistanceType::EA0_EB0;
@@ -139,7 +154,7 @@ EdgeEdgeDistanceType edge_edge_distance_type(
  } else {
  return EdgeEdgeDistanceType::EA_EB0;
  }
- } else if (tN >= tD) { // tc > 1  => the t=1 edge is visible
+ } else if (tN >= tD) { // tc > 1 ⟹ the t=1 edge is visible
  // recompute sc for this edge
  if ((-d + b) <= 0.0) {
  return EdgeEdgeDistanceType::EA0_EB1;
@@ -150,7 +165,44 @@ EdgeEdgeDistanceType edge_edge_distance_type(
  }
  }
 
- return defaultCase;
+ return default_case;
+}
+
+EdgeEdgeDistanceType edge_edge_parallel_distance_type(
+ const Eigen::Ref<const Eigen::Vector3d>& ea0,
+ const Eigen::Ref<const Eigen::Vector3d>& ea1,
+ const Eigen::Ref<const Eigen::Vector3d>& eb0,
+ const Eigen::Ref<const Eigen::Vector3d>& eb1)
+{
+ const Eigen::Vector3d ea = ea1 - ea0;
+ const double alpha = (eb0 - ea0).dot(ea) / ea.squaredNorm();
+ const double beta = (eb1 - eb0).dot(ea) / ea.squaredNorm();
+
+ uint8_t eac; // 0: EA0, 1: EA1, 2: EA
+ uint8_t ebc; // 0: EB0, 1: EB1, 2: EB
+ if (alpha < 0) {
+ eac = (0 <= beta && beta <= 1) ? 2 : 0;
+ ebc = (beta <= alpha) ? 0 : (beta <= 1 ? 1 : 2);
+ } else if (alpha > 1) {
+ eac = (0 <= beta && beta <= 1) ? 2 : 1;
+ ebc = (beta >= alpha) ? 0 : (0 <= beta ? 1 : 2);
+ } else {
+ eac = 2;
+ ebc = 0;
+ }
+
+ // f(0, 0) = 0000 = 0 -> EA0_EB0
+ // f(0, 1) = 0001 = 1 -> EA0_EB1
+ // f(1, 0) = 0010 = 2 -> EA1_EB0
+ // f(1, 1) = 0011 = 3 -> EA1_EB1
+ // f(2, 0) = 0100 = 4 -> EA_EB0
+ // f(2, 1) = 0101 = 5 -> EA_EB1
+ // f(0, 2) = 0110 = 6 -> EA0_EB
+ // f(1, 2) = 0111 = 7 -> EA1_EB
+ // f(2, 2) = 1000 = 8 -> EA_EB
+
+ assert(eac != 2 || ebc != 2); // This case results in a degenerate line-line
+ return EdgeEdgeDistanceType(ebc < 2 ? (eac << 1 | ebc) : (6 + eac));
 }
 
 } // namespace ipc

src/ipc/distance/distance_type.hpp

@@ -76,4 +76,16 @@ EdgeEdgeDistanceType edge_edge_distance_type(
  const Eigen::Ref<const Eigen::Vector3d>& eb0,
  const Eigen::Ref<const Eigen::Vector3d>& eb1);
 
+/// @brief Determine the closest pair between two parallel edges.
+/// @param ea0 The first vertex of the first edge.
+/// @param ea1 The second vertex of the first edge.
+/// @param eb0 The first vertex of the second edge.
+/// @param eb1 The second vertex of the second edge.
+/// @return The distance type of the edge-edge pair.
+EdgeEdgeDistanceType edge_edge_parallel_distance_type(
+ const Eigen::Ref<const Eigen::Vector3d>& ea0,
+ const Eigen::Ref<const Eigen::Vector3d>& ea1,
+ const Eigen::Ref<const Eigen::Vector3d>& eb0,
+ const Eigen::Ref<const Eigen::Vector3d>& eb1);
+
 } // namespace ipc

src/ipc/friction/friction_constraints.cpp

@@ -114,11 +114,7 @@ double FrictionConstraints::compute_potential(
  }
  });
 
- double potential = 0;
- for (const auto& local_potential : storage) {
- potential += local_potential;
- }
- return potential;
+ return storage.combine([](double a, double b) { return a + b; });
 }
 
 Eigen::VectorXd FrictionConstraints::compute_potential_gradient(
@@ -157,11 +153,8 @@ Eigen::VectorXd FrictionConstraints::compute_potential_gradient(
  }
  });
 
- Eigen::VectorXd grad = Eigen::VectorXd::Zero(ndof);
- for (const auto& local_grad : storage) {
- grad += local_grad;
- }
- return grad;
+ return storage.combine([](const Eigen::VectorXd& a,
+ const Eigen::VectorXd& b) { return a + b; });
 }
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -255,11 +248,8 @@ Eigen::VectorXd FrictionConstraints::compute_force(
  }
  });
 
- Eigen::VectorXd force = Eigen::VectorXd::Zero(velocities.size());
- for (const auto& local_force : storage) {
- force += local_force;
- }
- return force;
+ return storage.combine([](const Eigen::VectorXd& a,
+ const Eigen::VectorXd& b) { return a + b; });
 }
 
 ///////////////////////////////////////////////////////////////////////////////

src/ipc/implicits/plane.cpp

@@ -128,10 +128,8 @@ double compute_point_plane_collision_free_stepsize(
  }
  });
 
- double earliest_toi = 1;
- for (const auto& local_earliest_toi : storage) {
- earliest_toi = std::min(earliest_toi, local_earliest_toi);
- }
+ const double earliest_toi =
+ storage.combine([](double a, double b) { return std::min(a, b); });
  assert(earliest_toi >= 0 && earliest_toi <= 1.0);
  return earliest_toi;
 }

tests/distance/test_distance_type.cpp

@@ -386,30 +386,32 @@ TEST_CASE(
  "Edge-edge parallel distance type",
  "[distance][distance-type][edge-edge][parallel]")
 {
- double alpha = GENERATE(take(10, random(0.01, 0.99)));
- double s = GENERATE(take(10, random(-5.0, 5.0)));
+ const double alpha = GENERATE(take(10, random(0.01, 0.99)));
+ const double beta = GENERATE(take(10, random(1.01, 1.99)));
+ const double s = GENERATE(take(10, random(-5.0, 5.0)));
  const int n_random_edges = 20;
 
  for (int i = 0; i < n_random_edges; i++) {
  const Eigen::Vector3d ea0 = Eigen::Vector3d::Random();
  const Eigen::Vector3d ea1 = Eigen::Vector3d::Random();
- const double edge_len = (ea1 - ea0).norm();
+ const Eigen::Vector3d ea = ea1 - ea0;
 
  Eigen::Vector3d n = (ea1 - ea0).cross(Eigen::Vector3d::UnitX());
  REQUIRE(n.norm() != 0);
  n.normalize();
 
- const Eigen::Vector3d eb0 = (ea1 - ea0) * alpha + ea0 + s * n;
- const Eigen::Vector3d eb1 = (ea1 - ea0) + eb0;
- REQUIRE(
- (ea1 - ea0).dot(eb1 - eb0) == Catch::Approx(edge_len * edge_len));
+ const Eigen::Vector3d eb0 = ea0 + alpha * ea + s * n;
+ const Eigen::Vector3d eb1 = ea0 + beta * ea + s * n;
+ const Eigen::Vector3d eb = eb1 - eb0;
+
  REQUIRE(
- (ea1 - ea0).cross(eb1 - eb0).norm()
-   == Catch::Approx(0).margin(1e-14));
+ std::abs(ea.normalized().dot(eb.normalized())) == Catch::Approx(1));
+ REQUIRE(ea.cross(eb).norm() == Catch::Approx(0).margin(1e-14));
 
  const EdgeEdgeDistanceType dtype =
  edge_edge_distance_type(ea0, ea1, eb0, eb1);
 
+ CAPTURE(alpha, beta, s, dtype);
  CHECK(
  (dtype == EdgeEdgeDistanceType::EA_EB0
  || dtype == EdgeEdgeDistanceType::EA1_EB));

tests/distance/test_edge_edge.cpp

@@ -160,6 +160,13 @@ TEST_CASE("Edge-edge distance parallel", "[distance][edge-edge][parallel]")
 
  const double distance = edge_edge_distance(ea0, ea1, eb0, eb1);
  CHECK(distance == Catch::Approx(s * s).margin(1e-15));
+
+ for (int dtype = 0; dtype < int(EdgeEdgeDistanceType::EA_EB); dtype++) {
+ const double distance2 = edge_edge_distance(
+ ea0, ea1, eb0, eb1, EdgeEdgeDistanceType(dtype));
+ CAPTURE(dtype);
+ CHECK(distance <= Catch::Approx(distance2));
+ }
  }
 }