diff --git a/src/DEM/APIPrivate.cpp b/src/DEM/APIPrivate.cpp index 4a39993f..0f93147f 100644 --- a/src/DEM/APIPrivate.cpp +++ b/src/DEM/APIPrivate.cpp @@ -2142,8 +2142,13 @@ inline void DEMSolver::equipKernelIncludes(std::unordered_map (start offset, count), contact type -> list of [(program bundle name, +// kernel name)] +inline void DEMDynamicThread::dispatchCalcForceKernels( + const std::unordered_map>& typeStartCountMap, + const std::unordered_map, std::string>>>& + typeKernelMap) { + // For each contact type that exists, call its corresponding kernel(s) + for (size_t i = 0; i < m_numExistingTypes; i++) { + contact_t contact_type = existingContactTypes[i]; + const auto& start_count = typeStartCountMap.at(contact_type); + // Offset and count being contactPairs_t is very important, as CUDA kernel arguments cannot safely implicitly + // convert type (from size_t to unsigned int, for example) + contactPairs_t startOffset = start_count.first; + contactPairs_t count = start_count.second; + + // For this contact type, get its list of (program bundle name, kernel name) + if (typeKernelMap.count(contact_type) == 0) { + DEME_ERROR("Contact type %d has no associated force kernel in to execute!", contact_type); + } + const auto& kernelList = typeKernelMap.at(contact_type); + for (const auto& [progName, kernelName] : kernelList) { + size_t blocks = (count + DT_FORCE_CALC_NTHREADS_PER_BLOCK - 1) / DT_FORCE_CALC_NTHREADS_PER_BLOCK; + if (blocks > 0) { + progName->kernel(kernelName) + .instantiate() + .configure(dim3(blocks), dim3(DT_FORCE_CALC_NTHREADS_PER_BLOCK), 0, streamInfo.stream) + .launch(&simParams, &granData, startOffset, count); + } + } + } + DEME_GPU_CALL(cudaStreamSynchronize(streamInfo.stream)); +} + +void DEMDynamicThread::calculateForces() { // Reset force (acceleration) arrays for this time step size_t nContactPairs = *solverScratchSpace.numContacts; @@ -2171,18 +2204,14 @@ inline void DEMDynamicThread::calculateForces() { } timers.GetTimer("Clear force array").stop(); - size_t blocks_needed_for_contacts = - (nContactPairs + DT_FORCE_CALC_NTHREADS_PER_BLOCK - 1) / DT_FORCE_CALC_NTHREADS_PER_BLOCK; // If no contact then we don't have to calculate forces. Note there might still be forces, coming from prescription // or other sources. - if (blocks_needed_for_contacts > 0) { + if (nContactPairs > 0) { timers.GetTimer("Calculate contact forces").start(); - // a custom kernel to compute forces - cal_force_kernels->kernel("calculateContactForces") - .instantiate() - .configure(dim3(blocks_needed_for_contacts), dim3(DT_FORCE_CALC_NTHREADS_PER_BLOCK), 0, streamInfo.stream) - .launch(&simParams, &granData, nContactPairs); - DEME_GPU_CALL(cudaStreamSynchronize(streamInfo.stream)); + + // Call specialized kernels for each contact type that exists + dispatchCalcForceKernels(typeStartCountMap, contactTypeKernelMap); + // displayDeviceFloat3(granData->contactForces, nContactPairs); // displayDeviceArray(granData->contactType, nContactPairs); // std::cout << "===========================" << std::endl; @@ -2195,7 +2224,7 @@ inline void DEMDynamicThread::calculateForces() { collectContactForcesThruCub(collect_force_kernels, granData, nContactPairs, simParams->nOwnerBodies, contactPairArr_isFresh, streamInfo.stream, solverScratchSpace, timers); } else { - blocks_needed_for_contacts = + size_t blocks_needed_for_contacts = (nContactPairs + DEME_MAX_THREADS_PER_BLOCK - 1) / DEME_MAX_THREADS_PER_BLOCK; // This does both acc and ang acc collect_force_kernels->kernel("forceToAcc") @@ -2283,10 +2312,10 @@ inline void DEMDynamicThread::unpack_impl() { typeStartOffsets.toHost(); for (size_t i = 0; i < m_numExistingTypes; i++) { DEME_DEBUG_PRINTF("Contact type %d starts at offset %u", existingContactTypes[i], typeStartOffsets[i]); - typeStartCountMap[existingContactTypes[i]] = - std::make_pair(typeStartOffsets[i], - (i + 1 < m_numExistingTypes ? typeStartOffsets[i + 1] : *solverScratchSpace.numContacts) - - typeStartOffsets[i]); + typeStartCountMap[existingContactTypes[i]] = std::make_pair( + typeStartOffsets[i], + (i + 1 < m_numExistingTypes ? typeStartOffsets[i + 1] : (contactPairs_t)*solverScratchSpace.numContacts) - + typeStartOffsets[i]); } // Debug output of the map // for (const auto& entry : typeStartCountMap) { @@ -2579,6 +2608,18 @@ void DEMDynamicThread::jitifyKernels(const std::unordered_map(std::move(JitHelper::buildProgram( "DEMMiscKernels", JitHelper::KERNEL_DIR / "DEMMiscKernels.cu", Subs, JitifyOptions))); } + + // For now, the contact type to kernel map is known and hard-coded after jitification + contactTypeKernelMap = {// Sphere-Sphere contact + {SPHERE_SPHERE_CONTACT, {{cal_force_kernels, "calculateContactForces_SphSph"}}}, + // Sphere-Triangle contact + {SPHERE_TRIANGLE_CONTACT, {{cal_force_kernels, "calculateContactForces_SphTri"}}}, + // Sphere-Analytical contact + {SPHERE_ANALYTICAL_CONTACT, {{cal_force_kernels, "calculateContactForces_SphAnal"}}}, + // Triangle-Triangle contact + {TRIANGLE_TRIANGLE_CONTACT, {{cal_force_kernels, "calculateContactForces_TriTri"}}}, + // Triangle-Analytical contact + {TRIANGLE_ANALYTICAL_CONTACT, {{cal_force_kernels, "calculateContactForces_TriAnal"}}}}; } float* DEMDynamicThread::inspectCall(const std::shared_ptr& inspection_kernel, diff --git a/src/DEM/dT.h b/src/DEM/dT.h index 49efeaea..30dd363a 100644 --- a/src/DEM/dT.h +++ b/src/DEM/dT.h @@ -296,6 +296,9 @@ class DEMDynamicThread { size_t m_numExistingTypes = 0; // A map that records the contact for each contact type currently existing std::unordered_map> typeStartCountMap; + // A map that records the corresponding jitify program bundle and kernel name for each contact type + std::unordered_map, std::string>>> + contactTypeKernelMap; // dT's timers std::vector timer_names = {"Clear force array", "Calculate contact forces", "Optional force reduction", @@ -673,8 +676,13 @@ class DEMDynamicThread { // Migrate contact history to fit the structure of the newly received contact array inline void migrateEnduringContacts(); + // Impl of calculateForces + inline void dispatchCalcForceKernels( + const std::unordered_map>& typeStartCountMap, + const std::unordered_map, std::string>>>& + typeKernelMap); // Update clump-based acceleration array based on sphere-based force array - inline void calculateForces(); + void calculateForces(); // Update clump pos/oriQ and vel/omega based on acceleration inline void integrateOwnerMotions(); diff --git a/src/demo/DEMdemo_SingleSphereCollide.cpp b/src/demo/DEMdemo_SingleSphereCollide.cpp index 8afc9c54..d53b045a 100644 --- a/src/demo/DEMdemo_SingleSphereCollide.cpp +++ b/src/demo/DEMdemo_SingleSphereCollide.cpp @@ -115,8 +115,9 @@ int main() { // Testing modifying jitify options and force model prerequisites auto jitify_options = DEMSim.GetJitifyOptions(); - jitify_options.pop_back(); // Remove a warning suppression option - // DEMSim.SetJitifyOptions(jitify_options); + jitify_options.pop_back(); // Remove C++ std17 option + jitify_options.push_back("-std=c++20"); + DEMSim.SetJitifyOptions(jitify_options); // Then set it my_force_model->DefineCustomModelPrerequisites( "float3 __device__ GetContactForce(float3 AOwner, float3 BOwner, float3 ALinVel, float3 BLinVel, " "float3 ARotVel, float3 BRotVel, float delta_time, float delta_tan_x, float delta_tan_y, " diff --git a/src/kernel/DEMCalcForceKernels.cu b/src/kernel/DEMCalcForceKernels.cu index e2882baf..6e44509e 100644 --- a/src/kernel/DEMCalcForceKernels.cu +++ b/src/kernel/DEMCalcForceKernels.cu @@ -41,301 +41,359 @@ inline __device__ void equipOwnerPosRot(deme::DEMSimParams* simParams, bodyPos.z = ownerPos.z + (double)relPos.z; } -__global__ void calculateContactForces(deme::DEMSimParams* simParams, deme::DEMDataDT* granData, size_t nContactPairs) { - deme::contactPairs_t myContactID = blockIdx.x * blockDim.x + threadIdx.x; - if (myContactID < nContactPairs) { - // Identify contact type first - deme::contact_t ContactType = granData->contactType[myContactID]; - // The following quantities are always calculated, regardless of force model - double3 contactPnt; - float3 B2A; // Unit vector pointing from body B to body A (contact normal) - double overlapDepth; - double3 AOwnerPos, bodyAPos, BOwnerPos, bodyBPos; - float AOwnerMass, ARadius, BOwnerMass, BRadius; - float4 AOriQ, BOriQ; - deme::materialsOffset_t bodyAMatType, bodyBMatType; - // The user-specified extra margin size (how much we should be lenient in determining `in-contact') - float extraMarginSize = 0.; - // Triangle A's three points are defined outside, as may be reused in B's acquisition and penetration calc. - double3 triANode1, triANode2, triANode3; - // Then allocate the optional quantities that will be needed in the force model (note: this one can't be in a - // curly bracket, obviously...) - _forceModelIngredientDefinition_; - // Take care of 2 bodies in order, bodyA first, grab location and velocity to local cache - // Decompose ContactType to get the types of A and B - deme::geoType_t AType = deme::decodeTypeA(ContactType); - deme::geoType_t BType = deme::decodeTypeB(ContactType); - - // ---------------------------------------------------------------- - // Based on A's type, equip info - // ---------------------------------------------------------------- - if (AType == deme::GEO_T_SPHERE) { - deme::bodyID_t sphereID = granData->idGeometryA[myContactID]; - deme::bodyID_t myOwner = granData->ownerClumpBody[sphereID]; - - float3 myRelPos; - float myRadius; - // Get my component offset info from either jitified arrays or global memory - // Outputs myRelPos, myRadius - // Use an input named exactly `sphereID' which is the id of this sphere component - { _componentAcqStrat_; } - - // Get my mass info from either jitified arrays or global memory - // Outputs myMass - // Use an input named exactly `myOwner' which is the id of this owner - { - float myMass; - _massAcqStrat_; - AOwnerMass = myMass; - } +// Template device function for contact force calculation - will be called by 5 specialized kernels +template +__device__ __forceinline__ void calculateContactForcesImpl(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t myContactID) { + // Contact type is known at compile time + deme::contact_t ContactType = CONTACT_TYPE; + // The following quantities are always calculated, regardless of force model + double3 contactPnt; + float3 B2A; // Unit vector pointing from body B to body A (contact normal) + double overlapDepth; + double3 AOwnerPos, bodyAPos, BOwnerPos, bodyBPos; + float AOwnerMass, ARadius, BOwnerMass, BRadius; + float4 AOriQ, BOriQ; + deme::materialsOffset_t bodyAMatType, bodyBMatType; + // The user-specified extra margin size (how much we should be lenient in determining `in-contact') + float extraMarginSize = 0.; + // Triangle A's three points are defined outside, as may be reused in B's acquisition and penetration calc. + double3 triANode1, triANode2, triANode3; + // Then allocate the optional quantities that will be needed in the force model (note: this one can't be in a + // curly bracket, obviously...) + _forceModelIngredientDefinition_; + // Take care of 2 bodies in order, bodyA first, grab location and velocity to local cache + // Decompose ContactType to get the types of A and B (known at compile time) + constexpr deme::geoType_t AType = (CONTACT_TYPE >> 4); + constexpr deme::geoType_t BType = (CONTACT_TYPE & 0xF); + + // ---------------------------------------------------------------- + // Based on A's type, equip info + // ---------------------------------------------------------------- + if constexpr (AType == deme::GEO_T_SPHERE) { + deme::bodyID_t sphereID = granData->idGeometryA[myContactID]; + deme::bodyID_t myOwner = granData->ownerClumpBody[sphereID]; + + float3 myRelPos; + float myRadius; + // Get my component offset info from either jitified arrays or global memory + // Outputs myRelPos, myRadius + // Use an input named exactly `sphereID' which is the id of this sphere component + { _componentAcqStrat_; } + + // Get my mass info from either jitified arrays or global memory + // Outputs myMass + // Use an input named exactly `myOwner' which is the id of this owner + { + float myMass; + _massAcqStrat_; + AOwnerMass = myMass; + } - // Optional force model ingredients are loaded here... - _forceModelIngredientAcqForA_; - _forceModelGeoWildcardAcqForASph_; - - equipOwnerPosRot(simParams, granData, myOwner, myRelPos, AOwnerPos, bodyAPos, AOriQ); - - ARadius = myRadius; - bodyAMatType = granData->sphereMaterialOffset[sphereID]; - extraMarginSize = granData->familyExtraMarginSize[AOwnerFamily]; - } else if (AType == deme::GEO_T_TRIANGLE) { - // Geometry ID here is called sphereID, although it is not a sphere, it's more like triID. But naming it - // sphereID makes the acquisition process cleaner. - deme::bodyID_t sphereID = granData->idGeometryA[myContactID]; - deme::bodyID_t myOwner = granData->ownerMesh[sphereID]; - //// TODO: Is this OK? - ARadius = DEME_HUGE_FLOAT; - bodyAMatType = granData->triMaterialOffset[sphereID]; - - // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick - // the larger of the 2 familyExtraMarginSize. - extraMarginSize = granData->familyExtraMarginSize[AOwnerFamily]; - - triANode1 = to_double3(granData->relPosNode1[sphereID]); - triANode2 = to_double3(granData->relPosNode2[sphereID]); - triANode3 = to_double3(granData->relPosNode3[sphereID]); - - // Get my mass info from either jitified arrays or global memory - // Outputs myMass - // Use an input named exactly `myOwner' which is the id of this owner - { - float myMass; - _massAcqStrat_; - AOwnerMass = myMass; - } - _forceModelIngredientAcqForA_; - _forceModelGeoWildcardAcqForATri_; - - // bodyAPos is for a place holder for the outcome triANode1 position - equipOwnerPosRot(simParams, granData, myOwner, triANode1, AOwnerPos, bodyAPos, AOriQ); - triANode1 = bodyAPos; - // Do this to node 2 and 3 as well - applyOriQToVector3(triANode2.x, triANode2.y, triANode2.z, AOriQ.w, AOriQ.x, AOriQ.y, AOriQ.z); - triANode2 += AOwnerPos; - applyOriQToVector3(triANode3.x, triANode3.y, triANode3.z, AOriQ.w, AOriQ.x, AOriQ.y, AOriQ.z); - triANode3 += AOwnerPos; - // Assign the correct bodyAPos - bodyAPos = triangleCentroid(triANode1, triANode2, triANode3); - } else { - // Currently, we only support sphere and mesh for body A + // Optional force model ingredients are loaded here... + _forceModelIngredientAcqForA_; + _forceModelGeoWildcardAcqForASph_; + + equipOwnerPosRot(simParams, granData, myOwner, myRelPos, AOwnerPos, bodyAPos, AOriQ); + + ARadius = myRadius; + bodyAMatType = granData->sphereMaterialOffset[sphereID]; + extraMarginSize = granData->familyExtraMarginSize[AOwnerFamily]; + } else if constexpr (AType == deme::GEO_T_TRIANGLE) { + // Geometry ID here is called sphereID, although it is not a sphere, it's more like triID. But naming it + // sphereID makes the acquisition process cleaner. + deme::bodyID_t sphereID = granData->idGeometryA[myContactID]; + deme::bodyID_t myOwner = granData->ownerMesh[sphereID]; + //// TODO: Is this OK? + ARadius = DEME_HUGE_FLOAT; + bodyAMatType = granData->triMaterialOffset[sphereID]; + + // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick + // the larger of the 2 familyExtraMarginSize. + extraMarginSize = granData->familyExtraMarginSize[AOwnerFamily]; + + triANode1 = to_double3(granData->relPosNode1[sphereID]); + triANode2 = to_double3(granData->relPosNode2[sphereID]); + triANode3 = to_double3(granData->relPosNode3[sphereID]); + + // Get my mass info from either jitified arrays or global memory + // Outputs myMass + // Use an input named exactly `myOwner' which is the id of this owner + { + float myMass; + _massAcqStrat_; + AOwnerMass = myMass; + } + _forceModelIngredientAcqForA_; + _forceModelGeoWildcardAcqForATri_; + + // bodyAPos is for a place holder for the outcome triANode1 position + equipOwnerPosRot(simParams, granData, myOwner, triANode1, AOwnerPos, bodyAPos, AOriQ); + triANode1 = bodyAPos; + // Do this to node 2 and 3 as well + applyOriQToVector3(triANode2.x, triANode2.y, triANode2.z, AOriQ.w, AOriQ.x, AOriQ.y, AOriQ.z); + triANode2 += AOwnerPos; + applyOriQToVector3(triANode3.x, triANode3.y, triANode3.z, AOriQ.w, AOriQ.x, AOriQ.y, AOriQ.z); + triANode3 += AOwnerPos; + // Assign the correct bodyAPos + bodyAPos = triangleCentroid(triANode1, triANode2, triANode3); + } else { + // Currently, we only support sphere and mesh for body A + ContactType = deme::NOT_A_CONTACT; + } + + // ---------------------------------------------------------------- + // Then B, location and velocity, depending on type + // ---------------------------------------------------------------- + if constexpr (BType == deme::GEO_T_SPHERE) { + deme::bodyID_t sphereID = granData->idGeometryB[myContactID]; + deme::bodyID_t myOwner = granData->ownerClumpBody[sphereID]; + + float3 myRelPos; + float myRadius; + // Get my component offset info from either jitified arrays or global memory + // Outputs myRelPos, myRadius + // Use an input named exactly `sphereID' which is the id of this sphere component + { _componentAcqStrat_; } + + // Get my mass info from either jitified arrays or global memory + // Outputs myMass + // Use an input named exactly `myOwner' which is the id of this owner + { + float myMass; + _massAcqStrat_; + BOwnerMass = myMass; + } + _forceModelIngredientAcqForB_; + _forceModelGeoWildcardAcqForBSph_; + + equipOwnerPosRot(simParams, granData, myOwner, myRelPos, BOwnerPos, bodyBPos, BOriQ); + + BRadius = myRadius; + bodyBMatType = granData->sphereMaterialOffset[sphereID]; + + // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick + // the larger of the 2 familyExtraMarginSize. + extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) + ? extraMarginSize + : granData->familyExtraMarginSize[BOwnerFamily]; + + // If B is a sphere, then A can only be a sphere + checkSpheresOverlap(bodyAPos.x, bodyAPos.y, bodyAPos.z, ARadius, bodyBPos.x, bodyBPos.y, + bodyBPos.z, BRadius, contactPnt.x, contactPnt.y, contactPnt.z, B2A.x, B2A.y, + B2A.z, overlapDepth); + // If overlapDepth is negative then it might still be considered in contact, if the extra margins of A and B + // combined is larger than abs(overlapDepth) + if (overlapDepth < -extraMarginSize) { ContactType = deme::NOT_A_CONTACT; } - // ---------------------------------------------------------------- - // Then B, location and velocity, depending on type - // ---------------------------------------------------------------- - if (BType == deme::GEO_T_SPHERE) { - deme::bodyID_t sphereID = granData->idGeometryB[myContactID]; - deme::bodyID_t myOwner = granData->ownerClumpBody[sphereID]; - - float3 myRelPos; - float myRadius; - // Get my component offset info from either jitified arrays or global memory - // Outputs myRelPos, myRadius - // Use an input named exactly `sphereID' which is the id of this sphere component - { _componentAcqStrat_; } - - // Get my mass info from either jitified arrays or global memory - // Outputs myMass - // Use an input named exactly `myOwner' which is the id of this owner - { - float myMass; - _massAcqStrat_; - BOwnerMass = myMass; + } else if constexpr (BType == deme::GEO_T_TRIANGLE) { + // Geometry ID here is called sphereID, although it is not a sphere, it's more like triID. But naming it + // sphereID makes the acquisition process cleaner. + deme::bodyID_t sphereID = granData->idGeometryB[myContactID]; + deme::bodyID_t myOwner = granData->ownerMesh[sphereID]; + //// TODO: Is this OK? + BRadius = DEME_HUGE_FLOAT; + bodyBMatType = granData->triMaterialOffset[sphereID]; + + // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick + // the larger of the 2 familyExtraMarginSize. + extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) + ? extraMarginSize + : granData->familyExtraMarginSize[BOwnerFamily]; + // extraMarginSize here is purely family-based extra margin, so it can be used to determine if the user + // potentially needs remote (non-contact) force calculation. + bool needsNonContactPenetrationCalc = (extraMarginSize > 0.); + + double3 triBNode1 = to_double3(granData->relPosNode1[sphereID]); + double3 triBNode2 = to_double3(granData->relPosNode2[sphereID]); + double3 triBNode3 = to_double3(granData->relPosNode3[sphereID]); + + // Get my mass info from either jitified arrays or global memory + // Outputs myMass + // Use an input named exactly `myOwner' which is the id of this owner + { + float myMass; + _massAcqStrat_; + BOwnerMass = myMass; + } + _forceModelIngredientAcqForB_; + _forceModelGeoWildcardAcqForBTri_; + + // bodyBPos is for a place holder for the outcome triBNode1 position + equipOwnerPosRot(simParams, granData, myOwner, triBNode1, BOwnerPos, bodyBPos, BOriQ); + triBNode1 = bodyBPos; + // Do this to node 2 and 3 as well + applyOriQToVector3(triBNode2.x, triBNode2.y, triBNode2.z, BOriQ.w, BOriQ.x, BOriQ.y, BOriQ.z); + triBNode2 += BOwnerPos; + applyOriQToVector3(triBNode3.x, triBNode3.y, triBNode3.z, BOriQ.w, BOriQ.x, BOriQ.y, BOriQ.z); + triBNode3 += BOwnerPos; + // Assign the correct bodyBPos + bodyBPos = triangleCentroid(triBNode1, triBNode2, triBNode3); + + // If B is a triangle, then A can be a sphere or a triangle. + if constexpr (AType == deme::GEO_T_SPHERE) { + double3 contact_normal; + bool in_contact = checkTriSphereOverlap(triBNode1, triBNode2, triBNode3, bodyAPos, ARadius, + contact_normal, overlapDepth, contactPnt); + B2A = to_float3(contact_normal); + + // Sphere--triangle is a bit tricky. Extra margin should only take effect when it comes from the + // positive direction of the mesh facet. If not, sphere-setting-on-needle case will give huge + // penetration since in that case, overlapDepth is very negative and this will be considered in-contact. + // So the cases we exclude are: too far away while at the positive direction; not in contact while at + // the negative side. + // Also checkTriSphereOverlap gives positive number for overlapping cases + if ((overlapDepth < -extraMarginSize) || (!in_contact && overlapDepth > 0.)) { + ContactType = deme::NOT_A_CONTACT; } - _forceModelIngredientAcqForB_; - _forceModelGeoWildcardAcqForBSph_; - - equipOwnerPosRot(simParams, granData, myOwner, myRelPos, BOwnerPos, bodyBPos, BOriQ); - - BRadius = myRadius; - bodyBMatType = granData->sphereMaterialOffset[sphereID]; - - // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick - // the larger of the 2 familyExtraMarginSize. - extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) - ? extraMarginSize - : granData->familyExtraMarginSize[BOwnerFamily]; - - // If B is a sphere, then A can only be a sphere - checkSpheresOverlap(bodyAPos.x, bodyAPos.y, bodyAPos.z, ARadius, bodyBPos.x, bodyBPos.y, - bodyBPos.z, BRadius, contactPnt.x, contactPnt.y, contactPnt.z, B2A.x, - B2A.y, B2A.z, overlapDepth); - // If overlapDepth is negative then it might still be considered in contact, if the extra margins of A and B - // combined is larger than abs(overlapDepth) + } else if constexpr (AType == deme::GEO_T_TRIANGLE) { + // Triangle--triangle contact, a bit more complex... + double3 contact_normal; + checkTriangleTriangleOverlap(triANode1, triANode2, triANode3, triBNode1, triBNode2, + triBNode3, contact_normal, overlapDepth, contactPnt, + needsNonContactPenetrationCalc); + B2A = to_float3(contact_normal); + // Fix ContactType if needed if (overlapDepth < -extraMarginSize) { ContactType = deme::NOT_A_CONTACT; } + } - } else if (BType == deme::GEO_T_TRIANGLE) { - // Geometry ID here is called sphereID, although it is not a sphere, it's more like triID. But naming it - // sphereID makes the acquisition process cleaner. - deme::bodyID_t sphereID = granData->idGeometryB[myContactID]; - deme::bodyID_t myOwner = granData->ownerMesh[sphereID]; - //// TODO: Is this OK? - BRadius = DEME_HUGE_FLOAT; - bodyBMatType = granData->triMaterialOffset[sphereID]; - - // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick - // the larger of the 2 familyExtraMarginSize. - extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) - ? extraMarginSize - : granData->familyExtraMarginSize[BOwnerFamily]; - - double3 triBNode1 = to_double3(granData->relPosNode1[sphereID]); - double3 triBNode2 = to_double3(granData->relPosNode2[sphereID]); - double3 triBNode3 = to_double3(granData->relPosNode3[sphereID]); - - // Get my mass info from either jitified arrays or global memory - // Outputs myMass - // Use an input named exactly `myOwner' which is the id of this owner - { - float myMass; - _massAcqStrat_; - BOwnerMass = myMass; - } - _forceModelIngredientAcqForB_; - _forceModelGeoWildcardAcqForBTri_; - - // bodyBPos is for a place holder for the outcome triBNode1 position - equipOwnerPosRot(simParams, granData, myOwner, triBNode1, BOwnerPos, bodyBPos, BOriQ); - triBNode1 = bodyBPos; - // Do this to node 2 and 3 as well - applyOriQToVector3(triBNode2.x, triBNode2.y, triBNode2.z, BOriQ.w, BOriQ.x, BOriQ.y, BOriQ.z); - triBNode2 += BOwnerPos; - applyOriQToVector3(triBNode3.x, triBNode3.y, triBNode3.z, BOriQ.w, BOriQ.x, BOriQ.y, BOriQ.z); - triBNode3 += BOwnerPos; - // Assign the correct bodyBPos - bodyBPos = triangleCentroid(triBNode1, triBNode2, triBNode3); - - // If B is a triangle, then A can be a sphere or a triangle. But this branching is not too bad, as most - // threads in this block will have the same ContactType. - if (AType == deme::GEO_T_SPHERE) { - double3 contact_normal; - bool in_contact = checkTriSphereOverlap( - triBNode1, triBNode2, triBNode3, bodyAPos, ARadius, contact_normal, overlapDepth, contactPnt); - B2A = to_float3(contact_normal); - - // Sphere--triangle is a bit tricky. Extra margin should only take effect when it comes from the - // positive direction of the mesh facet. If not, sphere-setting-on-needle case will give huge - // penetration since in that case, overlapDepth is very negative and this will be considered in-contact. - // So the cases we exclude are: too far away while at the positive direction; not in contact while at - // the negative side. - // Also checkTriSphereOverlap gives positive number for overlapping cases - if ((overlapDepth < -extraMarginSize) || (!in_contact && overlapDepth > 0.)) { - ContactType = deme::NOT_A_CONTACT; - } - } else if (AType == deme::GEO_T_TRIANGLE) { - // Triangle--triangle contact is not supported yet, so we just skip it. + } else if constexpr (BType == deme::GEO_T_ANALYTICAL) { + // Geometry ID here is called sphereID, although it is not a sphere, it's more like analyticalID. But naming + // it sphereID makes the acquisition process cleaner. + deme::objID_t sphereID = granData->idGeometryB[myContactID]; + deme::bodyID_t myOwner = objOwner[sphereID]; + // If B is analytical entity, its owner, relative location, material info is jitified. + bodyBMatType = objMaterial[sphereID]; + BOwnerMass = objMass[sphereID]; + //// TODO: Is this OK? + BRadius = DEME_HUGE_FLOAT; + float3 myRelPos; + float3 bodyBRot; + myRelPos.x = objRelPosX[sphereID]; + myRelPos.y = objRelPosY[sphereID]; + myRelPos.z = objRelPosZ[sphereID]; + _forceModelIngredientAcqForB_; + _forceModelGeoWildcardAcqForBAnal_; + + equipOwnerPosRot(simParams, granData, myOwner, myRelPos, BOwnerPos, bodyBPos, BOriQ); + + // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick + // the larger of the 2 familyExtraMarginSize. + extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) + ? extraMarginSize + : granData->familyExtraMarginSize[BOwnerFamily]; + + // B's orientation (such as plane normal) is rotated with its owner too + bodyBRot.x = objRotX[sphereID]; + bodyBRot.y = objRotY[sphereID]; + bodyBRot.z = objRotZ[sphereID]; + applyOriQToVector3(bodyBRot.x, bodyBRot.y, bodyBRot.z, BOriQ.w, BOriQ.x, BOriQ.y, BOriQ.z); + + // If B is an analytical entity, then A can be a sphere or a triangle. + if constexpr (AType == deme::GEO_T_SPHERE) { + // Note for this test on dT side we don't enlarge entities + checkSphereEntityOverlap(bodyAPos, ARadius, objType[sphereID], bodyBPos, bodyBRot, + objSize1[sphereID], objSize2[sphereID], objSize3[sphereID], + objNormal[sphereID], 0.0, contactPnt, B2A, overlapDepth); + // Fix ContactType if needed + if (overlapDepth < -extraMarginSize) { ContactType = deme::NOT_A_CONTACT; } - - } else if (BType == deme::GEO_T_ANALYTICAL) { - // Geometry ID here is called sphereID, although it is not a sphere, it's more like analyticalID. But naming - // it sphereID makes the acquisition process cleaner. - deme::objID_t sphereID = granData->idGeometryB[myContactID]; - deme::bodyID_t myOwner = objOwner[sphereID]; - // If B is analytical entity, its owner, relative location, material info is jitified. - bodyBMatType = objMaterial[sphereID]; - BOwnerMass = objMass[sphereID]; - //// TODO: Is this OK? - BRadius = DEME_HUGE_FLOAT; - float3 myRelPos; - float3 bodyBRot; - myRelPos.x = objRelPosX[sphereID]; - myRelPos.y = objRelPosY[sphereID]; - myRelPos.z = objRelPosZ[sphereID]; - _forceModelIngredientAcqForB_; - _forceModelGeoWildcardAcqForBAnal_; - - equipOwnerPosRot(simParams, granData, myOwner, myRelPos, BOwnerPos, bodyBPos, BOriQ); - - // As the grace margin, the distance (negative overlap) just needs to be within the grace margin. So we pick - // the larger of the 2 familyExtraMarginSize. - extraMarginSize = (extraMarginSize > granData->familyExtraMarginSize[BOwnerFamily]) - ? extraMarginSize - : granData->familyExtraMarginSize[BOwnerFamily]; - - // B's orientation (such as plane normal) is rotated with its owner too - bodyBRot.x = objRotX[sphereID]; - bodyBRot.y = objRotY[sphereID]; - bodyBRot.z = objRotZ[sphereID]; - applyOriQToVector3(bodyBRot.x, bodyBRot.y, bodyBRot.z, BOriQ.w, BOriQ.x, BOriQ.y, - BOriQ.z); - - // If B is an analytical entity, then A can be a sphere or a triangle. But this branching is not too bad, as - // most threads in this block will have the same ContactType. - if (AType == deme::GEO_T_SPHERE) { - // Note for this test on dT side we don't enlarge entities - checkSphereEntityOverlap( - bodyAPos, ARadius, objType[sphereID], bodyBPos, bodyBRot, objSize1[sphereID], objSize2[sphereID], - objSize3[sphereID], objNormal[sphereID], 0.0, contactPnt, B2A, overlapDepth); - // Fix ContactType if needed - if (overlapDepth < -extraMarginSize) { - ContactType = deme::NOT_A_CONTACT; - } - } else if (AType == deme::GEO_T_TRIANGLE) { - calcTriEntityOverlap( - triANode1, triANode2, triANode3, objType[sphereID], bodyBPos, bodyBRot, objSize1[sphereID], - objSize2[sphereID], objSize3[sphereID], objNormal[sphereID], contactPnt, B2A, overlapDepth); - // Fix ContactType if needed - if (overlapDepth < -extraMarginSize) { - ContactType = deme::NOT_A_CONTACT; - } + } else if constexpr (AType == deme::GEO_T_TRIANGLE) { + calcTriEntityOverlap(triANode1, triANode2, triANode3, objType[sphereID], bodyBPos, + bodyBRot, objSize1[sphereID], objSize2[sphereID], objSize3[sphereID], + objNormal[sphereID], contactPnt, B2A, overlapDepth); + // Fix ContactType if needed + if (overlapDepth < -extraMarginSize) { + ContactType = deme::NOT_A_CONTACT; } } + } - _forceModelContactWildcardAcq_; - if (ContactType != deme::NOT_A_CONTACT) { - float3 force = make_float3(0, 0, 0); - float3 torque_only_force = make_float3(0, 0, 0); - // Local position of the contact point is always a piece of info we require... regardless of force model - float3 locCPA = to_float3(contactPnt - AOwnerPos); - float3 locCPB = to_float3(contactPnt - BOwnerPos); - // Now map this contact point location to bodies' local ref - applyOriQToVector3(locCPA.x, locCPA.y, locCPA.z, AOriQ.w, -AOriQ.x, -AOriQ.y, - -AOriQ.z); - applyOriQToVector3(locCPB.x, locCPB.y, locCPB.z, BOriQ.w, -BOriQ.x, -BOriQ.y, - -BOriQ.z); - // The following part, the force model, is user-specifiable - // NOTE!! "force" and all wildcards must be properly set by this piece of code - { _DEMForceModel_; } - - // Write contact location values back to global memory - _contactInfoWrite_; - - // If force model modifies owner wildcards, write them back here - _forceModelOwnerWildcardWrite_; - - // Optionally, the forces can be reduced to acc right here (may be faster) - _forceCollectInPlaceStrat_; - } else { - // The contact is no longer active, so we need to destroy its contact history recording - _forceModelContactWildcardDestroy_; - } + _forceModelContactWildcardAcq_; + if (ContactType != deme::NOT_A_CONTACT) { + float3 force = make_float3(0, 0, 0); + float3 torque_only_force = make_float3(0, 0, 0); + // Local position of the contact point is always a piece of info we require... regardless of force model + float3 locCPA = to_float3(contactPnt - AOwnerPos); + float3 locCPB = to_float3(contactPnt - BOwnerPos); + // Now map this contact point location to bodies' local ref + applyOriQToVector3(locCPA.x, locCPA.y, locCPA.z, AOriQ.w, -AOriQ.x, -AOriQ.y, -AOriQ.z); + applyOriQToVector3(locCPB.x, locCPB.y, locCPB.z, BOriQ.w, -BOriQ.x, -BOriQ.y, -BOriQ.z); + // The following part, the force model, is user-specifiable + // NOTE!! "force" and all wildcards must be properly set by this piece of code + { _DEMForceModel_; } + + // Write contact location values back to global memory + _contactInfoWrite_; + + // If force model modifies owner wildcards, write them back here + _forceModelOwnerWildcardWrite_; + + // Optionally, the forces can be reduced to acc right here (may be faster) + _forceCollectInPlaceStrat_; + } else { + // The contact is no longer active, so we need to destroy its contact history recording + _forceModelContactWildcardDestroy_; + } + + // Updated contact wildcards need to be write back to global mem. It is here because contact wildcard may need + // to be destroyed for non-contact, so it has to go last. + _forceModelContactWildcardWrite_; +} + +// 5 specialized kernels for different contact types +__global__ void calculateContactForces_SphSph(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t startOffset, + deme::contactPairs_t nContactPairs) { + deme::contactPairs_t myContactID = startOffset + blockIdx.x * blockDim.x + threadIdx.x; + if (myContactID < startOffset + nContactPairs) { + calculateContactForcesImpl(simParams, granData, myContactID); + } +} + +__global__ void calculateContactForces_SphTri(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t startOffset, + deme::contactPairs_t nContactPairs) { + deme::contactPairs_t myContactID = startOffset + blockIdx.x * blockDim.x + threadIdx.x; + if (myContactID < startOffset + nContactPairs) { + calculateContactForcesImpl(simParams, granData, myContactID); + } +} + +__global__ void calculateContactForces_SphAnal(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t startOffset, + deme::contactPairs_t nContactPairs) { + deme::contactPairs_t myContactID = startOffset + blockIdx.x * blockDim.x + threadIdx.x; + if (myContactID < startOffset + nContactPairs) { + calculateContactForcesImpl(simParams, granData, myContactID); + } +} + +__global__ void calculateContactForces_TriTri(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t startOffset, + deme::contactPairs_t nContactPairs) { + deme::contactPairs_t myContactID = startOffset + blockIdx.x * blockDim.x + threadIdx.x; + if (myContactID < startOffset + nContactPairs) { + calculateContactForcesImpl(simParams, granData, myContactID); + } +} - // Updated contact wildcards need to be write back to global mem. It is here because contact wildcard may need - // to be destroyed for non-contact, so it has to go last. - _forceModelContactWildcardWrite_; +__global__ void calculateContactForces_TriAnal(deme::DEMSimParams* simParams, + deme::DEMDataDT* granData, + deme::contactPairs_t startOffset, + deme::contactPairs_t nContactPairs) { + deme::contactPairs_t myContactID = startOffset + blockIdx.x * blockDim.x + threadIdx.x; + if (myContactID < startOffset + nContactPairs) { + calculateContactForcesImpl(simParams, granData, myContactID); } } diff --git a/src/kernel/DEMCollisionKernels.cuh b/src/kernel/DEMCollisionKernels.cuh index 63b49e0f..0cb98740 100644 --- a/src/kernel/DEMCollisionKernels.cuh +++ b/src/kernel/DEMCollisionKernels.cuh @@ -485,10 +485,347 @@ inline __device__ bool checkTriangleTriangleOverlap(const T1& A1, const T1& A2, const T1& B2, const T1& C2, - T1& normal, ///< contact normal - T2& depth, ///< penetration (positive if in contact) - T1& point) { ///< contact point - return false; + T1& normal, ///< contact normal + T2& depth, ///< penetration (positive if in contact) + T1& point, ///< contact point + bool outputNoContact) { ///< output info even when no contact + // Triangle 1 vertices + const T1 tri1[3] = {A1, B1, C1}; + // Triangle 2 vertices + const T1 tri2[3] = {A2, B2, C2}; + + // Compute face normals + T1 edge1_0 = B1 - A1; + T1 edge1_1 = C1 - B1; + T1 edge1_2 = A1 - C1; + T1 normal1 = cross(edge1_0, C1 - A1); + + T1 edge2_0 = B2 - A2; + T1 edge2_1 = C2 - B2; + T1 edge2_2 = A2 - C2; + T1 normal2 = cross(edge2_0, C2 - A2); + + // Track the MTV (minimum translation vector) or maximum separation + T2 minOverlap = DEME_HUGE_FLOAT; + T1 mtv_axis; + int axisType = -1; // 0: face1, 1: face2, 2+: edge-edge + int edgeIndexA = -1, edgeIndexB = -1; + bool foundSeparation = false; + T2 maxSeparation = -DEME_HUGE_FLOAT; // Track largest separation gap + + // Store edges for later use + T1 edges1[3] = {edge1_0, edge1_1, edge1_2}; + T1 edges2[3] = {edge2_0, edge2_1, edge2_2}; + + // Test face normal of triangle 1 + { + T2 len2 = dot(normal1, normal1); + if (len2 > DEME_TINY_FLOAT) { + T1 axis = normal1 * rsqrt(len2); + + // Project triangle 1 vertices + T2 min1 = dot(tri1[0], axis); + T2 max1 = min1; + for (int i = 1; i < 3; ++i) { + T2 proj = dot(tri1[i], axis); + if (proj < min1) + min1 = proj; + if (proj > max1) + max1 = proj; + } + + // Project triangle 2 vertices + T2 min2 = dot(tri2[0], axis); + T2 max2 = min2; + for (int i = 1; i < 3; ++i) { + T2 proj = dot(tri2[i], axis); + if (proj < min2) + min2 = proj; + if (proj > max2) + max2 = proj; + } + + // Check overlap or separation + if (max1 < min2 || max2 < min1) { + // Separating axis found + T2 separation = (max1 < min2) ? (min2 - max1) : (min1 - max2); + if (separation > maxSeparation) { + maxSeparation = separation; + mtv_axis = (max1 < min2) ? axis : (axis * T2(-1.0)); + axisType = 0; + foundSeparation = true; + } + if (!outputNoContact) + return false; + } else { + // Calculate overlap + T2 overlap = DEME_MIN(max1 - min2, max2 - min1); + if (overlap < minOverlap) { + minOverlap = overlap; + mtv_axis = axis; + axisType = 0; + // Ensure MTV points from tri2 to tri1 + if (max1 - min2 < max2 - min1) { + mtv_axis = axis * T2(-1.0); + } + } + } + } + } + + // Test face normal of triangle 2 + { + T2 len2 = dot(normal2, normal2); + if (len2 > DEME_TINY_FLOAT) { + T1 axis = normal2 * rsqrt(len2); + + // Project triangle 1 vertices + T2 min1 = dot(tri1[0], axis); + T2 max1 = min1; + for (int i = 1; i < 3; ++i) { + T2 proj = dot(tri1[i], axis); + if (proj < min1) + min1 = proj; + if (proj > max1) + max1 = proj; + } + + // Project triangle 2 vertices + T2 min2 = dot(tri2[0], axis); + T2 max2 = min2; + for (int i = 1; i < 3; ++i) { + T2 proj = dot(tri2[i], axis); + if (proj < min2) + min2 = proj; + if (proj > max2) + max2 = proj; + } + + // Check overlap or separation + if (max1 < min2 || max2 < min1) { + // Separating axis found + T2 separation = (max1 < min2) ? (min2 - max1) : (min1 - max2); + if (separation > maxSeparation) { + maxSeparation = separation; + mtv_axis = (max1 < min2) ? axis : (axis * T2(-1.0)); + axisType = 1; + foundSeparation = true; + } + if (!outputNoContact) + return false; + } else { + // Calculate overlap + T2 overlap = DEME_MIN(max1 - min2, max2 - min1); + if (overlap < minOverlap) { + minOverlap = overlap; + mtv_axis = axis; + axisType = 1; + // Ensure MTV points from tri2 to tri1 + if (max1 - min2 < max2 - min1) { + mtv_axis = axis * T2(-1.0); + } + } + } + } + } + + // Test 9 edge-edge cross products + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + T1 axis = cross(edges1[i], edges2[j]); + T2 len2 = dot(axis, axis); + + if (len2 > DEME_TINY_FLOAT) { + axis = axis * rsqrt(len2); + + // Project triangle 1 vertices + T2 min1 = dot(tri1[0], axis); + T2 max1 = min1; + for (int k = 1; k < 3; ++k) { + T2 proj = dot(tri1[k], axis); + if (proj < min1) + min1 = proj; + if (proj > max1) + max1 = proj; + } + + // Project triangle 2 vertices + T2 min2 = dot(tri2[0], axis); + T2 max2 = min2; + for (int k = 1; k < 3; ++k) { + T2 proj = dot(tri2[k], axis); + if (proj < min2) + min2 = proj; + if (proj > max2) + max2 = proj; + } + + // Check overlap or separation + if (max1 < min2 || max2 < min1) { + // Separating axis found + T2 separation = (max1 < min2) ? (min2 - max1) : (min1 - max2); + if (separation > maxSeparation) { + maxSeparation = separation; + mtv_axis = (max1 < min2) ? axis : (axis * T2(-1.0)); + axisType = 2; + edgeIndexA = i; + edgeIndexB = j; + foundSeparation = true; + } + if (!outputNoContact) + return false; + } else { + // Calculate overlap + T2 overlap = DEME_MIN(max1 - min2, max2 - min1); + if (overlap < minOverlap) { + minOverlap = overlap; + mtv_axis = axis; + axisType = 2; + edgeIndexA = i; + edgeIndexB = j; + // Ensure MTV points from tri2 to tri1 + if (max1 - min2 < max2 - min1) { + mtv_axis = axis * T2(-1.0); + } + } + } + } + } + } + + // Safety check: if no valid axis was found (degenerate triangles), return false + if (axisType == -1) { + return false; + } + + // Determine if contact exists + bool inContact = !foundSeparation; + + // Set output depth and normal + if (inContact) { + depth = minOverlap; + } else { + depth = -maxSeparation; // Negative for no contact + } + normal = mtv_axis; + + // Calculate contact point as centroid of penetrating region + // Use a polygon clipping approach similar to tri_plane_penetration + if (axisType == 0 || axisType == 1) { + // Face-based contact: clip one triangle against the other's plane + const T1* refTri = (axisType == 0) ? tri1 : tri2; + const T1* incTri = (axisType == 0) ? tri2 : tri1; + T1 planeNormal = (axisType == 0) ? normalize(normal1) : normalize(normal2); + + // Compute signed distances + T2 d[3]; + for (int i = 0; i < 3; ++i) { + d[i] = dot(incTri[i] - refTri[0], planeNormal); + } + + // Build clipped polygon + T1 poly[6]; // Maximum 6 vertices for clipped triangle + int nNode = 0; + + for (int i = 0; i < 3; ++i) { + int j = (i + 1) % 3; + bool in_i = (d[i] < 0.0); + bool in_j = (d[j] < 0.0); + + // Add vertex i if it's inside + if (in_i) { + poly[nNode++] = incTri[i]; + } + + // Add intersection point if edge crosses plane + if (in_i != in_j) { + T2 t = d[i] / (d[i] - d[j]); + T1 inter = incTri[i] + (incTri[j] - incTri[i]) * t; + poly[nNode++] = inter; + } + } + + // Compute centroid of clipped polygon + T1 centroid; + centroid.x = 0.0; + centroid.y = 0.0; + centroid.z = 0.0; + + if (nNode > 0) { + for (int i = 0; i < nNode; i++) { + centroid = centroid + poly[i]; + } + centroid = centroid / T2(nNode); + } else { + // Fallback: use centroid of incident triangle + centroid = (incTri[0] + incTri[1] + incTri[2]) / T2(3.0); + } + + // Project centroid onto plane if in contact, or compute midpoint trajectory if not + if (inContact) { + T2 centroidDist = dot(centroid - refTri[0], planeNormal); + point = centroid - planeNormal * (centroidDist + depth * T2(0.5)); + } else { + // For no contact: midpoint between centroids along separation axis + T1 cent1 = (tri1[0] + tri1[1] + tri1[2]) / T2(3.0); + T1 cent2 = (tri2[0] + tri2[1] + tri2[2]) / T2(3.0); + point = (cent1 + cent2) * T2(0.5); + } + } else { + // Edge-edge: compute closest points between the winning edge pair + int nextA = (edgeIndexA + 1) % 3; + T1 edgeA_start = tri1[edgeIndexA]; + T1 edgeA = edges1[edgeIndexA]; + + int nextB = (edgeIndexB + 1) % 3; + T1 edgeB_start = tri2[edgeIndexB]; + T1 edgeB = edges2[edgeIndexB]; + + // Compute closest points on two line segments + T1 r = edgeB_start - edgeA_start; + T2 a = dot(edgeA, edgeA); + T2 e = dot(edgeB, edgeB); + T2 f = dot(edgeB, r); + + T2 s, t; + if (a <= DEME_TINY_FLOAT && e <= DEME_TINY_FLOAT) { + s = t = 0.0; + } else if (a <= DEME_TINY_FLOAT) { + s = 0.0; + t = clampBetween(f / e, T2(0.0), T2(1.0)); + } else { + T2 c = dot(edgeA, r); + if (e <= DEME_TINY_FLOAT) { + t = 0.0; + s = clampBetween(-c / a, T2(0.0), T2(1.0)); + } else { + T2 b = dot(edgeA, edgeB); + T2 denom = a * e - b * b; + + if (denom != T2(0.0)) { + s = clampBetween((b * f - c * e) / denom, T2(0.0), T2(1.0)); + } else { + s = 0.0; + } + + t = (b * s + f) / e; + + if (t < T2(0.0)) { + t = 0.0; + s = clampBetween(-c / a, T2(0.0), T2(1.0)); + } else if (t > T2(1.0)) { + t = 1.0; + s = clampBetween((b - c) / a, T2(0.0), T2(1.0)); + } + } + } + + T1 closestA = edgeA_start + edgeA * s; + T1 closestB = edgeB_start + edgeB * t; + point = (closestA + closestB) * T2(0.5); + } + + return inContact; } #endif diff --git a/src/kernel/DEMMiscKernels.cu b/src/kernel/DEMMiscKernels.cu index 918dc103..b96dadcf 100644 --- a/src/kernel/DEMMiscKernels.cu +++ b/src/kernel/DEMMiscKernels.cu @@ -45,8 +45,10 @@ __global__ void computeMarginFromAbsv(deme::DEMSimParams* simParams, unsigned int my_family = granData->familyID[ownerID]; if (!isfinite(absv)) { // May produce messy error messages, but it's still good to know what entities went wrong - DEME_ABORT_KERNEL("Absolute velocity for ownerID %llu is not finite. This happened at time %.9g.\n", - static_cast(ownerID), simParams->timeElapsed); + DEME_ABORT_KERNEL( + "Absolute velocity for ownerID %llu is infinite (and it's a worse version of " + "max-velocity-exceeded-allowance).\n", + static_cast(ownerID)); } if (absv > simParams->approxMaxVel) { absv = simParams->approxMaxVel;