Implement SAT-based triangle-triangle collision detection and split contact force kernels by type

[Copilot]

Overview

This PR implements the checkTriangleTriangleOverlap CUDA function for narrow-phase collision detection between two triangles, based on the Separating Axis Theorem (SAT). The function was previously a stub that always returned false. Additionally, this PR refactors the force calculation system by splitting the monolithic calculateContactForces kernel into 5 specialized kernels based on contact types.

Triangle-Triangle Collision Detection Implementation

The function implements a complete SAT-based collision detection algorithm with support for both contact and no-contact scenarios:

Separating Axis Tests

• 2 face normals: One from each triangle to detect face-face penetration
• 9 edge-edge cross products: All combinations of edges from both triangles (3×3) to detect edge-edge contact

Minimum Translation Vector (MTV) and Separation Tracking

The implementation tracks both overlap and separation across all 11 tested axes to determine:

• Contact normal: Points from triangle 2 to triangle 1 (or along maximum separation axis for no-contact)
• Penetration depth: Positive value indicates contact exists; negative value indicates separation distance
• Optional no-contact output: When outputNoContact flag is true, provides information even when triangles don't overlap

Contact Point Calculation

The contact point is computed as the centroid of the penetrating polygon, similar to tri_plane_penetration:

• Face-based contact:
  • Clips incident triangle against reference triangle's plane
  • Builds polygon representing the penetrating region (up to 6 vertices)
  • Computes centroid of clipped polygon as contact point
  • Adjusts with half-depth offset for stability
• Edge-edge contact: Computes closest points between the winning edge pair using the standard line segment distance algorithm
• No-contact case (when outputNoContact=true):
  • Returns negative depth (separation distance)
  • Normal points along maximum separation axis (expected contact trajectory)
  • Contact point is midpoint between triangle centroids

Robustness Features

• Handles degenerate triangles with safety checks
• Guards against parallel or near-parallel edges using DEME_TINY_FLOAT epsilon
• Returns false if no valid separating axis can be determined
• Tracks both minimum overlap (contact) and maximum separation (no-contact) simultaneously

Function Signature

template <typename T1, typename T2>
inline __device__ bool checkTriangleTriangleOverlap(
    const T1& A1, const T1& B1, const T1& C1,  // Triangle 1 vertices
    const T1& A2, const T1& B2, const T1& C2,  // Triangle 2 vertices
    T1& normal,            // Contact normal
    T2& depth,             // Penetration (positive) or separation (negative)
    T1& point,             // Contact point
    bool outputNoContact)  // Output info even when no contact

Contact Force Kernel Refactoring

Split the monolithic calculateContactForces kernel into 5 specialized kernels for improved performance:

Kernel Architecture

• Template device function: calculateContactForcesImpl<CONTACT_TYPE> contains the core force calculation logic
• 5 Specialized global kernels: Each handles one specific contact type
  1. calculateContactForces_SphSph - Sphere-Sphere contacts
  2. calculateContactForces_SphTri - Sphere-Triangle contacts
  3. calculateContactForces_SphAnal - Sphere-Analytical contacts
  4. calculateContactForces_TriTri - Triangle-Triangle contacts
  5. calculateContactForces_TriAnal - Triangle-Analytical contacts

Benefits of Split Kernels

• Compile-time optimization: Contact type known at compile time via template parameter enables better compiler optimization
• Reduced branching: Eliminates runtime type checking and conditional branching within kernels
• Type-specific tuning: Each kernel processes uniform workload of single contact type
• Better GPU occupancy: Reduced divergence improves warp efficiency

Kernel Invocation

• Modified dT.cpp to call all 5 kernels conditionally based on typeStartCountMap
• Each kernel receives startOffset and nContactPairs for its specific contact type range
• Only launches kernels for contact types that exist in current timestep
• All kernels synchronized on same CUDA stream for proper ordering

Technical Details

Template Parameters

• T1: Vector type (typically double3)
• T2: Scalar type (typically double)

Utility Functions Used

Leverages existing functions from DEMHelperKernels.cuh and CUDAMathHelpers.cuh:

• cross(), dot(), normalize(), clamp(), sqrt()
• Constants: DEME_MIN(), DEME_TINY_FLOAT, DEME_HUGE_FLOAT

Code Quality

• Follows existing code patterns and style in the repository, particularly tri_plane_penetration
• Formatted with clang-format according to project standards
• Comprehensive inline comments explaining the algorithm

Testing

While CUDA is not available in the development environment for compilation testing, the implementation:

• Uses only existing, well-tested utility functions from the codebase
• Follows established patterns from similar collision detection functions (e.g., calc_prism_contact, tri_plane_penetration)
• Implements standard, well-documented algorithms from computational geometry literature

Use Cases

This PR enables:

• Triangle-triangle contact detection for mesh-mesh interactions in DEM simulations
• Triangulated surface collision handling
• Fine-grained contact detection between complex geometries represented as triangle meshes
• Distance-based force calculations using the no-contact output mode
• Performance-optimized force calculations with reduced kernel divergence

Original prompt

Please implement the CUDA function checkTriangleTriangleOverlap, based on the comments about this function that we already have there. You can assume T1 is double3 and T2 is double in your implementation. Try using the existing utility functions in DEMHelperKernels.cuh and DEMCollisionKernels.cuh if possible.

---

✨ Let Copilot coding agent set things up for you — coding agent works faster and does higher quality work when set up for your repo.

[Ruochun]

@copilot This is good. But due to the nature of this solver, I need you to make the following changes to checkTriangleTriangleOverlap. 1. The contact penetration is evaluated as the max penetration depth, which we already do; but the contact point needs to be roughly at the center of the "penetrating polygon", that is similar to the implementation in "tri_plane_penetration" function. This contact point does not need to be stable (but best if it is), but should be roughly representing the center of the penetrating volume. 2. Even in the case of no-contact, if a flag is set to true (this should be extra function argument), we will need the function to output the penetration depth, contact point and contact normal. This depth in this case should be negative and roughly represents the distance "before contact happens". This contact point should be some representation of the midpoint of this "expected contact trajectory" and the normal is the direction of this trajectory. You can decide how this can be implemented and this can be a estimation and does not need to be extremely accurate.

[Copilot]

@copilot This is good. But due to the nature of this solver, I need you to make the following changes to checkTriangleTriangleOverlap. 1. The contact penetration is evaluated as the max penetration de...

I've made both requested changes in commit 90901bd:

1. Contact point as center of penetrating polygon: Changed from simple vertex projection to a polygon clipping approach similar to tri_plane_penetration. The function now:

   • Clips the incident triangle against the reference triangle's plane
   • Computes the centroid of the resulting clipped polygon (the penetrating region)
   • Projects this centroid with a half-depth offset for the contact point

2. Added outputNoContact parameter: When set to true, the function:

   • Continues testing all axes even after finding separation
   • Tracks the maximum separation distance (axis with largest gap)
   • Outputs negative depth representing distance before contact
   • Sets normal as the direction of the expected contact trajectory (along max separation axis)
   • Computes contact point as the midpoint between triangle centroids

The implementation tracks both minOverlap (for contact) and maxSeparation (for no-contact cases) during SAT testing, allowing it to provide meaningful outputs in both scenarios.

[Ruochun]

@copilot Now, can you split the force calculation kernel (calculateContactForces) into 5, based on the 5 contact types? The types, and their corresponding starts and numbers of elements in the contact array is already stored in typeStartCountMap. The split kernels can still live in the same kernel file, but obviously the 5 different kernels should all be called in each time step, and each of them only takes care of one type of contact.

[Copilot]

@copilot Now, can you split the force calculation kernel (calculateContactForces) into 5, based on the 5 contact types? The types, and their corresponding starts and numbers of elements in the contact...

I've split the force calculation kernel into 5 specialized kernels in commit 0bb471f:

Kernel Structure:

• Created template device function calculateContactForcesImpl<CONTACT_TYPE> with the core logic
• 5 global kernels wrap this template for each contact type:
  • calculateContactForces_SphSph
  • calculateContactForces_SphTri
  • calculateContactForces_SphAnal
  • calculateContactForces_TriTri
  • calculateContactForces_TriAnal

Calling Code (dT.cpp):

• Replaced single kernel call with 5 conditional calls based on typeStartCountMap
• Each kernel receives startOffset and nContactPairs for its specific contact type range
• Only launches kernels for contact types that actually exist
• All kernels use the same CUDA stream for proper synchronization

The contact type is now known at compile time (via template parameter), enabling better compiler optimization and eliminating runtime branching for type checks.