diff --git a/VERSION b/VERSION index 29cec99e..25875f01 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.18 +3.0.19 diff --git a/src/DEM/API.h b/src/DEM/API.h index db4f8a13..658cf33f 100644 --- a/src/DEM/API.h +++ b/src/DEM/API.h @@ -290,6 +290,29 @@ class DEMSolver { /// @details A negative value disables this guard when the corresponding contact kernel support is enabled. void SetTriTriContactRejectionRatio(float ratio); + /// @brief Set the overall multiplier for automatic multi-patch center offsets (default 1). + /// @details Before Initialize() only. Automatic centers use c - d*n, where c is the unweighted mean of + /// triangle centroids, n is the area-weighted unit normal, and d = alpha*max(gamma*h, beta*L). + /// h = abs(dot(c,n)) is the normal distance to the implicit local CoM (0,0,0); L is the current mesh-local + /// bounding-box diagonal, including unused vertices. Bounds and normals refresh with deformation when requested. + /// These advanced heuristic settings apply to mesh initialization and automatic deformation + /// refreshes. Single-patch origins, explicit user centers and mass properties are unaffected. Cancelling or + /// degenerate normals retain the unshifted average. Larger offsets make the contact rejection guard more + /// permissive; this is not a guarantee of valid filtering for arbitrary concave geometry or inconsistent winding. + /// @param alpha Finite, non-negative scale. Zero restores the historical unshifted multi-patch averages. + void SetMeshPatchCenterOffsetAlpha(float alpha); + /// @brief Set the minimum automatic center offset relative to mesh size (default 0.01). + /// @details Before Initialize() only. In d = alpha*max(gamma*h, beta*L), beta supplies a fallback when the + /// CoM lies near the patch plane. See SetMeshPatchCenterOffsetAlpha for geometry, scope and limitations. + /// @param beta Finite, non-negative fraction of the current local mesh bounding-box diagonal. Zero disables + /// the mesh-size fallback, so coplanar patches with h=0 remain unshifted. + void SetMeshPatchCenterOffsetBeta(float beta); + /// @brief Set the multiplier on the CoM normal distance for automatic centers (default 0.8). + /// @details Before Initialize() only. In d = alpha*max(gamma*h, beta*L), gamma scales only h, not the + /// mesh-size fallback. See SetMeshPatchCenterOffsetAlpha for geometry, scope and limitations. + /// @param gamma Finite, non-negative scale. Zero leaves only the beta*L contribution before alpha scaling. + void SetMeshPatchCenterOffsetGamma(float gamma); + /// @brief Used to force the solver to error out when there are too many spheres in a bin. A huge number can be used /// to discourage this error type. /// @param max_sph Max number of spheres in a bin. @@ -871,8 +894,9 @@ class DEMSolver { /// Device-input counterparts of the local node setters. Supply one float3 per mesh node, in original order. /// Calls are synchronous; finish producing source on its CUDA stream before calling. Topology is fixed. /// validate=false skips pointer/finite-value checks. update_patch_centers refreshes automatic centers only: - /// one patch uses the local origin, multiple patches average triangle centroids. Explicit centers are preserved. - /// Device updates lazily synchronize CPU geometry through GetCachedMesh, tracker getters, or solver output. + /// one patch uses the local origin, multiple patches use inward-offset centroid averages. Explicit centers are + /// preserved. Device updates lazily synchronize CPU geometry through GetCachedMesh, tracker getters, or solver + /// output. void SetTriNodeRelPosFromDevice(bodyID_t owner, const float3* source, int source_device, diff --git a/src/DEM/APIPublic.cpp b/src/DEM/APIPublic.cpp index f903fe9e..98a8344f 100644 --- a/src/DEM/APIPublic.cpp +++ b/src/DEM/APIPublic.cpp @@ -1646,6 +1646,28 @@ void DEMSolver::SetMaxTriTriPenetration(double max_margin) { } } +// Setup-only settings remain on dT for mesh initialization and device deformation refreshes. +void DEMSolver::SetMeshPatchCenterOffsetAlpha(float alpha) { + assertSysNotInit("SetMeshPatchCenterOffsetAlpha"); + if (!std::isfinite(alpha) || alpha < 0) + DEME_ERROR("Patch center offset alpha must be finite and non-negative."); + dT->meshPatchCenterOptions.alpha = alpha; +} + +void DEMSolver::SetMeshPatchCenterOffsetBeta(float beta) { + assertSysNotInit("SetMeshPatchCenterOffsetBeta"); + if (!std::isfinite(beta) || beta < 0) + DEME_ERROR("Patch center offset beta must be finite and non-negative."); + dT->meshPatchCenterOptions.beta = beta; +} + +void DEMSolver::SetMeshPatchCenterOffsetGamma(float gamma) { + assertSysNotInit("SetMeshPatchCenterOffsetGamma"); + if (!std::isfinite(gamma) || gamma < 0) + DEME_ERROR("Patch center offset gamma must be finite and non-negative."); + dT->meshPatchCenterOptions.gamma = gamma; +} + void DEMSolver::SetTriTriContactRejectionRatio(float ratio) { m_triTriContactRejectionRatio = ratio; if (!sys_initialized) { diff --git a/src/DEM/AuxClasses.h b/src/DEM/AuxClasses.h index 10ada198..06128667 100644 --- a/src/DEM/AuxClasses.h +++ b/src/DEM/AuxClasses.h @@ -419,7 +419,7 @@ class DEMTracker { /// @param update_patch_centers Recompute automatic centers from the full current geometry. False skips only /// this call; a later true call includes accumulated deformation. This option has NO effect for user-supplied /// centers, which are never overwritten. Automatic single-patch centers remain at the local origin; multiple - /// patches use the unweighted mean of triangle centroids. This does not change the contact-rejection heuristic. + /// patches use inward-offset centroid averages with the solver offset settings. The rejection rule is unchanged. void UpdateMeshFromDevice(const float3* source, int source_device, bool validate = true, diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index 1c44b308..ff62d2f0 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -764,10 +764,17 @@ class DEMMesh : public DEMInitializer { bool ArePatchLocationsExplicitlySet() const { return patch_locations_explicitly_set; } /// @brief Compute patch locations relative to the implicit CoM of the mesh. - /// @details For single patch: returns (0,0,0). For multiple patches: returns average of triangle centroids per - /// patch. + /// @details Single patches return (0,0,0). Multiple patches start at the unweighted mean triangle centroid c, + /// then move opposite the area-weighted normal n by alpha*max(gamma*abs(dot(c,n)), beta*L), where L is the + /// current local mesh bounding-box diagonal and the implicit CoM is the local origin. Cancelling or degenerate + /// normals leave c unchanged. This computes automatic geometry only; it does not overwrite explicit centers. + /// @param alpha Overall offset multiplier (default 1); zero restores unshifted averages. + /// @param beta Minimum offset as a fraction of L before alpha scaling (default 0.01). + /// @param gamma Multiplier on the CoM normal distance (default 0.8). + /// All parameters must be finite and non-negative. Pass the solver settings explicitly when previewing custom + /// settings; this standalone helper does not inherit settings from a solver. /// @return Vector of locations (one per patch). - std::vector ComputePatchLocations() const; + std::vector ComputePatchLocations(float alpha = 1.f, float beta = 0.01f, float gamma = 0.8f) const; enum class PatchQualityLevel : uint8_t { SAFE = 0, WARN = 1, CRITICAL = 2 }; diff --git a/src/DEM/dT.cpp b/src/DEM/dT.cpp index 5df0fe9f..40c9b7b1 100644 --- a/src/DEM/dT.cpp +++ b/src/DEM/dT.cpp @@ -1108,7 +1108,8 @@ void DEMDynamicThread::populateEntityArrays(const std::vectorpatch_locations_explicitly_set) { this_mesh_patch_locations = input_mesh_objs.at(i)->m_patch_locations; } else { - this_mesh_patch_locations = input_mesh_objs.at(i)->ComputePatchLocations(); + this_mesh_patch_locations = input_mesh_objs.at(i)->ComputePatchLocations( + meshPatchCenterOptions.alpha, meshPatchCenterOptions.beta, meshPatchCenterOptions.gamma); } // mesh_patch_owner run length is the num of patches in this mesh entity @@ -4663,6 +4664,7 @@ void DEMDynamicThread::refreshMeshPatchCenters(bodyID_t owner) { RefreshMeshPatchCenters(state.patch_offsets.data(), state.patch_triangles.data(), state.patch_offsets.size() - 1, relPosNode1.data() + state.tri_start, relPosNode2.data() + state.tri_start, relPosNode3.data() + state.tri_start, relPosPatch.data() + state.patch_start, + state.vertices.data(), state.vertices.size(), state.diagonal.data(), meshPatchCenterOptions, streamInfo.stream); state.centers_dirty = true; } diff --git a/src/DEM/dT.h b/src/DEM/dT.h index 6d7420af..810d3f0f 100644 --- a/src/DEM/dT.h +++ b/src/DEM/dT.h @@ -382,6 +382,8 @@ class DEMDynamicThread { SolverTimers timers = SolverTimers(timer_names); public: + // Setup-only policy reused for mesh initialization and every automatic deformation refresh. + MeshPatchCenterOptions meshPatchCenterOptions; // Invalidates renderer geometry after initialization, insertion, or mesh deformation (not rigid motion). std::uint64_t visualizationRevision = 0; friend class DEMSolver; diff --git a/src/DEM/python/bindings.cpp b/src/DEM/python/bindings.cpp index 5b2fe887..d5f59fbe 100644 --- a/src/DEM/python/bindings.cpp +++ b/src/DEM/python/bindings.cpp @@ -1403,6 +1403,27 @@ enabled. A zero count is a no-op.)doc", .def("SetMaxTriTriPenetration", &deme::DEMSolver::SetMaxTriTriPenetration, "Set the maximum allowed triangle-triangle penetration used as the margin added in kT contact detection. " "This value caps the penetration margin added to prevent excessively large values.") + .def("SetMeshPatchCenterOffsetAlpha", &deme::DEMSolver::SetMeshPatchCenterOffsetAlpha, + "Before Initialize only. Default 1. Overall scale; zero restores unshifted averages. " + "Automatic multi-patch centers use c-d*n, d=alpha*max(gamma*abs(dot(c,n)), beta*L), " + "with local CoM at the origin and L the current local mesh bounding-box diagonal. " + "Finite non-negative values only. Applies to mesh initialization and deformation refreshes; " + "single-patch origins and explicit centers are preserved. Larger offsets admit deeper contacts.", + py::arg("alpha")) + .def("SetMeshPatchCenterOffsetBeta", &deme::DEMSolver::SetMeshPatchCenterOffsetBeta, + "Before Initialize only. Default 0.01. Mesh-size fallback fraction; zero disables the fallback. " + "Automatic multi-patch centers use c-d*n, d=alpha*max(gamma*abs(dot(c,n)), beta*L), " + "with local CoM at the origin and L the current local mesh bounding-box diagonal. " + "Finite non-negative values only. Applies to mesh initialization and deformation refreshes; " + "single-patch origins and explicit centers are preserved. Larger offsets admit deeper contacts.", + py::arg("beta")) + .def("SetMeshPatchCenterOffsetGamma", &deme::DEMSolver::SetMeshPatchCenterOffsetGamma, + "Before Initialize only. Default 0.8. CoM normal-distance scale; zero uses only the fallback. " + "Automatic multi-patch centers use c-d*n, d=alpha*max(gamma*abs(dot(c,n)), beta*L), " + "with local CoM at the origin and L the current local mesh bounding-box diagonal. " + "Finite non-negative values only. Applies to mesh initialization and deformation refreshes; " + "single-patch origins and explicit centers are preserved. Larger offsets admit deeper contacts.", + py::arg("gamma")) .def("SetTriTriContactRejectionRatio", &deme::DEMSolver::SetTriTriContactRejectionRatio, "Set the ratio threshold for rejecting suspicious tri-tri contacts. A contact is rejected when the " "penetration depth exceeds this fraction of the contact-point-to-mesh-center distance for either mesh " @@ -2417,6 +2438,10 @@ enabled. A zero count is a no-op.)doc", "Assign every triangle to its own patch using consecutive patch IDs.") .def("GetPatchIDs", &deme::DEMMeshConnected::GetPatchIDs, "Get one patch ID per triangle.") .def("GetNumPatches", &deme::DEMMeshConnected::GetNumPatches, "Get the number of mesh patches.") + .def("ComputePatchLocations", &deme::DEMMeshConnected::ComputePatchLocations, + "Compute automatic local patch centers using explicit offset settings (solver settings are not " + "inherited). Single patches return the origin; explicit centers are not modified.", + py::arg("alpha") = 1.f, py::arg("beta") = 0.01f, py::arg("gamma") = 0.8f) .def("SetPatchLocations", &deme::DEMMeshConnected::SetPatchLocations, "Set local patch centers before initialization. Use tracker UpdateMeshPatchLocations at runtime.", py::arg("centers")) diff --git a/src/DEM/utils/MeshDeformation.hpp b/src/DEM/utils/MeshDeformation.hpp index 623143a0..5a444acf 100644 --- a/src/DEM/utils/MeshDeformation.hpp +++ b/src/DEM/utils/MeshDeformation.hpp @@ -4,6 +4,7 @@ #define DEME_MESH_DEFORMATION_HPP #include +#include "algorithms/DEMMeshPatchCenter.h" #include "core/utils/DataMigrationHelper.hpp" #include "DEM/BdrsAndObjs.h" @@ -18,7 +19,8 @@ struct MeshDeformationState { faces(host_bytes, device_bytes), patch_offsets(host_bytes, device_bytes), patch_triangles(host_bytes, device_bytes), - invalid(host_bytes, device_bytes) {} + invalid(host_bytes, device_bytes), + diagonal(host_bytes, device_bytes) {} std::shared_ptr mesh; size_t tri_start = 0; size_t patch_start = 0; @@ -32,6 +34,7 @@ struct MeshDeformationState { patch_offsets.toDevice(); patch_triangles.toDevice(); invalid.toDevice(); + diagonal.resizeDevice(1); input.resizeDevice(input.size()); device_ready = true; } @@ -44,6 +47,8 @@ struct MeshDeformationState { DualArray patch_offsets; DualArray patch_triangles; DualArray invalid; + // One scalar computed on the same stream before patch reduction; never downloaded during deformation. + DualArray diagonal; }; } // namespace deme #endif diff --git a/src/DEM/utils/MeshUtils.cpp b/src/DEM/utils/MeshUtils.cpp index 21bb05c5..d47e4d25 100644 --- a/src/DEM/utils/MeshUtils.cpp +++ b/src/DEM/utils/MeshUtils.cpp @@ -46,6 +46,7 @@ #include "DEM/BdrsAndObjs.h" #include "DEM/utils/HostSideHelpers.hpp" #include "DEM/utils/MeshUtils.hpp" +#include "algorithms/DEMMeshPatchCenter.h" #include "core/utils/Logger.hpp" #include "core/utils/WavefrontMeshLoader.hpp" #include "kernel/DEMHelperKernels.cuh" @@ -1079,49 +1080,43 @@ void DEMMesh::SetEachTriangleAsPatch() { SetPatchIDs(patch_ids); } -// Compute patch locations (relative to CoM, which is implicitly at 0,0,0) -// If not explicitly set, calculates as: -// - Single patch: (0,0,0) -// - Multiple patches: average of triangle centroids per patch -std::vector DEMMesh::ComputePatchLocations() const { - std::vector patch_locations(nPatches, make_float3(0, 0, 0)); - - if (nPatches == 1) { - // Single patch: location is just CoM (0,0,0) - return patch_locations; - } - - // Multiple patches: compute average of triangle centroids per patch - std::vector patch_triangle_counts(nPatches, 0); - - for (size_t i = 0; i < nTri; ++i) { - const int3& face = m_face_v_indices[i]; - const float3& v0 = m_vertices[face.x]; - const float3& v1 = m_vertices[face.y]; - const float3& v2 = m_vertices[face.z]; - - float3 centroid = triangleCentroid(v0, v1, v2); - patchID_t patch_id = (i < m_patch_ids.size()) ? m_patch_ids[i] : 0; - - // Validate patch_id is within bounds - if (patch_id >= 0 && patch_id < static_cast(nPatches)) { - patch_locations[patch_id].x += centroid.x; - patch_locations[patch_id].y += centroid.y; - patch_locations[patch_id].z += centroid.z; - patch_triangle_counts[patch_id]++; +// Match deformation's corrected winding and shared host/device offset policy. Bounds use all current vertices, +// including unused ones, so setup and device refresh have the same definition of mesh size. +std::vector DEMMesh::ComputePatchLocations(float alpha, float beta, float gamma) const { + if (!std::isfinite(alpha) || !std::isfinite(beta) || !std::isfinite(gamma) || alpha < 0 || beta < 0 || gamma < 0) + DEME_ERROR("Patch center offset parameters must be finite and non-negative."); + std::vector locations(nPatches, make_float3(0, 0, 0)); + if (nPatches <= 1) + return locations; + double lo[3] = {INFINITY, INFINITY, INFINITY}, hi[3] = {-INFINITY, -INFINITY, -INFINITY}; + for (const auto& v : m_vertices) { + const double xyz[3] = {v.x, v.y, v.z}; + for (int k = 0; k < 3; ++k) { + lo[k] = std::min(lo[k], xyz[k]); + hi[k] = std::max(hi[k], xyz[k]); } } - - // Average the accumulated centroids - for (unsigned int p = 0; p < nPatches; ++p) { - if (patch_triangle_counts[p] > 0) { - patch_locations[p].x /= patch_triangle_counts[p]; - patch_locations[p].y /= patch_triangle_counts[p]; - patch_locations[p].z /= patch_triangle_counts[p]; + double diagonal2 = 0; + if (!m_vertices.empty()) + for (int k = 0; k < 3; ++k) + diagonal2 += (hi[k] - lo[k]) * (hi[k] - lo[k]); + std::vector> sums(nPatches, std::array{}); + std::vector counts(nPatches, 0); + for (size_t i = 0; i < nTri; ++i) { + const int3 face = m_face_v_indices[i]; + float3 a = m_vertices[face.x], b = m_vertices[face.y], c = m_vertices[face.z]; + // Supplied normals correct winding at initialization. The device subsequently keeps this topology. + if (use_mesh_normals && dot(cross(b - a, c - a), m_normals.at(m_face_n_indices.at(i).x)) < 0) + std::swap(b, c); + const patchID_t p = i < m_patch_ids.size() ? m_patch_ids[i] : 0; + if (p >= 0 && p < static_cast(nPatches)) { + AccumulatePatchTriangle(sums[p].data(), a, b, c); + ++counts[p]; } } - - return patch_locations; + for (unsigned int p = 0; p < nPatches; ++p) + locations[p] = FinishPatchCenter(sums[p].data(), counts[p], sqrt(diagonal2), {alpha, beta, gamma}); + return locations; } // Compute volume, centroid and MOI in CoM frame (unit density). diff --git a/src/algorithms/CMakeLists.txt b/src/algorithms/CMakeLists.txt index 4c4ced9f..f5531f4e 100644 --- a/src/algorithms/CMakeLists.txt +++ b/src/algorithms/CMakeLists.txt @@ -25,6 +25,7 @@ set(algorithms_interface ${CMAKE_CURRENT_SOURCE_DIR}/DEMStaticDeviceSubroutines.h ${CMAKE_CURRENT_SOURCE_DIR}/DEMOwnerDataRetrieval.h ${CMAKE_CURRENT_SOURCE_DIR}/DEMMeshDeformation.h + ${CMAKE_CURRENT_SOURCE_DIR}/DEMMeshPatchCenter.h ) ### INTERNAL HEADERS ONLY (.h, .hpp, or .cuh) ### diff --git a/src/algorithms/DEMMeshDeformation.cu b/src/algorithms/DEMMeshDeformation.cu index 635fdbb6..c8f9c034 100644 --- a/src/algorithms/DEMMeshDeformation.cu +++ b/src/algorithms/DEMMeshDeformation.cu @@ -52,15 +52,54 @@ __global__ void scatterVertices(const float3* vertices, node3[i] = vertices[f.z]; } -// One block reduces each patch's fixed triangle list. Strided loads keep large patches parallel without atomics; -// the fixed reduction tree makes repeated updates deterministic. +// A single strided reduction computes whole-mesh bounds once, rather than rescanning the mesh per patch. +// The persistent output scalar and same-stream ordering avoid allocations and host synchronization per refresh. +__global__ void meshDiagonal(const float3* vertices, size_t count, double* diagonal) { + __shared__ double bounds[6][MESH_DEFORMATION_BLOCK]; + const unsigned int lane = threadIdx.x; + double lo[3] = {INFINITY, INFINITY, INFINITY}, hi[3] = {-INFINITY, -INFINITY, -INFINITY}; + for (size_t i = lane; i < count; i += blockDim.x) { + const double xyz[3] = {vertices[i].x, vertices[i].y, vertices[i].z}; + for (int k = 0; k < 3; ++k) { + lo[k] = fmin(lo[k], xyz[k]); + hi[k] = fmax(hi[k], xyz[k]); + } + } + for (int k = 0; k < 3; ++k) { + bounds[k][lane] = lo[k]; + bounds[k + 3][lane] = hi[k]; + } + __syncthreads(); + for (unsigned int stride = MESH_DEFORMATION_BLOCK / 2; stride; stride /= 2) { + if (lane < stride) + for (int k = 0; k < 3; ++k) { + bounds[k][lane] = fmin(bounds[k][lane], bounds[k][lane + stride]); + bounds[k + 3][lane] = fmax(bounds[k + 3][lane], bounds[k + 3][lane + stride]); + } + __syncthreads(); + } + if (lane == 0) { + double squared = 0; + if (count) + for (int k = 0; k < 3; ++k) { + const double extent = bounds[k + 3][0] - bounds[k][0]; + squared += extent * extent; + } + *diagonal = sqrt(squared); + } +} + +// One block reduces each patch's fixed triangle list. The fixed tree is deterministic and avoids atomics. +// The shared helper documents the geometry assumptions and is also used for CPU initialization. __global__ void patchCenters(const size_t* offsets, const size_t* triangles, size_t patches, const float3* node1, const float3* node2, const float3* node3, - float3* centers) { + float3* centers, + const double* diagonal, + MeshPatchCenterOptions options) { const size_t p = blockIdx.x; const unsigned int lane = threadIdx.x; const size_t count = offsets[p + 1] - offsets[p]; @@ -69,28 +108,25 @@ __global__ void patchCenters(const size_t* offsets, centers[p] = make_float3(0, 0, 0); return; } - __shared__ double sums[3][MESH_DEFORMATION_BLOCK]; - double x = 0, y = 0, z = 0; + __shared__ double sums[7][MESH_DEFORMATION_BLOCK]; + double local[7] = {}; for (size_t j = offsets[p] + lane; j < offsets[p + 1]; j += blockDim.x) { const size_t t = triangles[j]; - x += ((double)node1[t].x + node2[t].x + node3[t].x) / 3.; - y += ((double)node1[t].y + node2[t].y + node3[t].y) / 3.; - z += ((double)node1[t].z + node2[t].z + node3[t].z) / 3.; + AccumulatePatchTriangle(local, node1[t], node2[t], node3[t]); } - sums[0][lane] = x; - sums[1][lane] = y; - sums[2][lane] = z; + for (int k = 0; k < 7; ++k) + sums[k][lane] = local[k]; __syncthreads(); - for (unsigned int stride = MESH_DEFORMATION_BLOCK / 2; stride > 0; stride /= 2) { - if (lane < stride) { - sums[0][lane] += sums[0][lane + stride]; - sums[1][lane] += sums[1][lane + stride]; - sums[2][lane] += sums[2][lane + stride]; - } + for (unsigned int stride = MESH_DEFORMATION_BLOCK / 2; stride; stride /= 2) { + if (lane < stride) + for (int k = 0; k < 7; ++k) + sums[k][lane] += sums[k][lane + stride]; __syncthreads(); } if (lane == 0) { - centers[p] = make_float3((float)(sums[0][0] / count), (float)(sums[1][0] / count), (float)(sums[2][0] / count)); + for (int k = 0; k < 7; ++k) + local[k] = sums[k][0]; + centers[p] = FinishPatchCenter(local, count, *diagonal, options); } } } // namespace @@ -136,10 +172,18 @@ void RefreshMeshPatchCenters(const size_t* offsets, const float3* node2, const float3* node3, float3* centers, + const float3* vertices, + size_t vertex_count, + double* diagonal, + MeshPatchCenterOptions options, cudaStream_t stream) { + if (patches > 1) { + meshDiagonal<<<1, MESH_DEFORMATION_BLOCK, 0, stream>>>(vertices, vertex_count, diagonal); + DEME_GPU_CALL(cudaGetLastError()); + } if (patches) { patchCenters<<>>(offsets, triangles, patches, node1, node2, node3, - centers); + centers, diagonal, options); DEME_GPU_CALL(cudaGetLastError()); } } diff --git a/src/algorithms/DEMMeshDeformation.h b/src/algorithms/DEMMeshDeformation.h index d80a66ad..05f8fcd5 100644 --- a/src/algorithms/DEMMeshDeformation.h +++ b/src/algorithms/DEMMeshDeformation.h @@ -4,6 +4,7 @@ #define DEME_MESH_DEFORMATION_H #include #include +#include "algorithms/DEMMeshPatchCenter.h" namespace deme { // Validate input (and the result of an incremental update) before mutating persistent geometry. void ValidateMeshVectors(const float3* input, @@ -22,7 +23,7 @@ void DeformMeshVertices(const float3* input, float3* node2, float3* node3, cudaStream_t stream); -// Preserve the historical local-origin convention for a single patch; otherwise average triangle centroids. +// Preserve the single-patch origin; otherwise apply the shared inward-offset policy using current device geometry. void RefreshMeshPatchCenters(const size_t* offsets, const size_t* triangles, size_t patches, @@ -30,6 +31,10 @@ void RefreshMeshPatchCenters(const size_t* offsets, const float3* node2, const float3* node3, float3* centers, + const float3* vertices, + size_t vertex_count, + double* diagonal, + MeshPatchCenterOptions options, cudaStream_t stream); } // namespace deme #endif diff --git a/src/algorithms/DEMMeshPatchCenter.h b/src/algorithms/DEMMeshPatchCenter.h new file mode 100644 index 00000000..4193d8eb --- /dev/null +++ b/src/algorithms/DEMMeshPatchCenter.h @@ -0,0 +1,64 @@ +// Copyright (c) 2026, SBEL GPU Development Team +// SPDX-License-Identifier: BSD-3-Clause +#ifndef DEME_MESH_PATCH_CENTER_H +#define DEME_MESH_PATCH_CENTER_H +#include +#include +#include + +namespace deme { +// Shared host/device policy. These settings belong to the solver, not mass properties or the contact kernel. +struct MeshPatchCenterOptions { + float alpha = 1.f; + float beta = 0.01f; + float gamma = 0.8f; +}; + +// Accumulate centroid coordinates, oriented twice-area vector, and total twice-area in double precision. +// Cross products are evaluated after promotion so small triangles and large coordinates do not overflow float +// intermediates. Keeping the scalar area separately lets us detect cancellation independently of mesh scale. +__host__ __device__ inline void AccumulatePatchTriangle(double* sum, float3 a, float3 b, float3 c) { + const double ux = double(b.x) - a.x, uy = double(b.y) - a.y, uz = double(b.z) - a.z; + const double vx = double(c.x) - a.x, vy = double(c.y) - a.y, vz = double(c.z) - a.z; + const double nx = uy * vz - uz * vy, ny = uz * vx - ux * vz, nz = ux * vy - uy * vx; + sum[0] += (double(a.x) + b.x + c.x) / 3.; + sum[1] += (double(a.y) + b.y + c.y) / 3.; + sum[2] += (double(a.z) + b.z + c.z) / 3.; + sum[3] += nx; + sum[4] += ny; + sum[5] += nz; + sum[6] += sqrt(nx * nx + ny * ny + nz * nz); +} + +// A surface centroid is a poor reference for the depth/center-distance rejection guard: legitimate contacts +// close to that centroid have almost zero permitted penetration. Move automatic multi-patch centers behind +// their oriented surface instead. Both an exterior face and an inward-facing container wall move along -n; +// testing which side contains the CoM is unnecessary. Local coordinates have their implicit CoM at the origin. +// h = |(CoM-c).n| uses only normal distance, avoiding offsets caused by tangential displacement from the CoM. +// d = alpha * max(gamma*h, beta*L), where L is the CURRENT whole-mesh local bounding-box diagonal (including +// unused vertices). The fallback handles planar meshes whose CoM lies in the plane, and does not shrink with +// triangle subdivision. This is a heuristic, not a solid-interior or visibility test: inconsistent winding, +// folded/concave patches and thin walls can still make a center inappropriate. Larger offsets also admit more +// deep contacts. Explicit centers bypass this policy, and alpha=0 restores the historical averaged centers. +// A resultant normal below 1e-6 of the sum of triangle areas is unreliable; leave such a patch unshifted. +// Single-patch meshes are handled by callers and retain their local-origin center. No mass properties change. +__host__ __device__ inline float3 FinishPatchCenter(const double* sum, + size_t count, + double diagonal, + MeshPatchCenterOptions options) { + if (!count) + return make_float3(0, 0, 0); + double x = sum[0] / count, y = sum[1] / count, z = sum[2] / count; + const double norm = sqrt(sum[3] * sum[3] + sum[4] * sum[4] + sum[5] * sum[5]); + if (options.alpha > 0 && norm > 1.e-6 * sum[6]) { + const double nx = sum[3] / norm, ny = sum[4] / norm, nz = sum[5] / norm; + const double h = fabs(x * nx + y * ny + z * nz); + const double d = double(options.alpha) * fmax(double(options.gamma) * h, double(options.beta) * diagonal); + x -= d * nx; + y -= d * ny; + z -= d * nz; + } + return make_float3(float(x), float(y), float(z)); +} +} // namespace deme +#endif diff --git a/src/demo/DEMdemo_MixerCubes.cpp b/src/demo/DEMdemo_MixerCubes.cpp index 0f1e4001..c2ac9dfd 100644 --- a/src/demo/DEMdemo_MixerCubes.cpp +++ b/src/demo/DEMdemo_MixerCubes.cpp @@ -138,7 +138,7 @@ int main() { out_dir /= "DemoOutput_MixerCubes"; create_directory(out_dir); - float sim_end = 10.0; + float sim_end = 5.0; unsigned int fps = 20; float frame_time = 1.0 / fps; diff --git a/src/demo/ModularTests/DEMTest_MeshDeformation.cpp b/src/demo/ModularTests/DEMTest_MeshDeformation.cpp index 557a03d3..328c608b 100644 --- a/src/demo/ModularTests/DEMTest_MeshDeformation.cpp +++ b/src/demo/ModularTests/DEMTest_MeshDeformation.cpp @@ -9,6 +9,7 @@ #include #include "DEM/API.h" +#include "algorithms/DEMMeshDeformation.h" #include "core/utils/CudaDebugSync.hpp" using namespace deme; @@ -48,8 +49,50 @@ struct DeviceInput { void checkAutomatic(const std::shared_ptr& tracker, const std::vector& expected) { auto mesh = tracker->GetMesh(); compare(mesh->GetCoordsVertices(), expected, "device/host vertex mismatch"); - compare(mesh->GetPatchLocations(), mesh->ComputePatchLocations(), "automatic GPU patch centers are stale"); + compare(mesh->GetPatchLocations(), mesh->ComputePatchLocations(0.7f, 0.03f, 1.f), + "automatic GPU patch centers are stale"); } +// Compare the CUDA reduction to CPU geometry for a patch spanning multiple block strides, opposing normals, +// and a zero-area patch. Include an unused vertex so both paths must use whole-mesh bounds, not patch bounds. +void checkCenterReduction(int device) { + ScopedCudaDevice scope(device); + DEMMesh mesh; + mesh.m_vertices = {make_float3(0, 0, 2), make_float3(1, 0, 2), make_float3(0, 1, 2), make_float3(5, 0, 0)}; + mesh.m_face_v_indices.assign(257, make_int3(0, 1, 2)); + mesh.m_patch_ids.assign(257, 0); + mesh.m_face_v_indices.insert(mesh.m_face_v_indices.end(), + {make_int3(0, 1, 2), make_int3(0, 2, 1), make_int3(0, 0, 0)}); + mesh.m_patch_ids.insert(mesh.m_patch_ids.end(), {1, 1, 2}); + mesh.nTri = mesh.m_face_v_indices.size(); + mesh.nPatches = 3; + DualArray vertices(mesh.m_vertices), a(mesh.nTri), b(mesh.nTri), c(mesh.nTri), centers(size_t(3)); + DualArray offsets(std::vector{0, 257, 259, 260}), triangles(mesh.nTri); + DualArray diagonal(size_t(1)); + for (size_t i = 0; i < mesh.nTri; ++i) { + const auto face = mesh.m_face_v_indices[i]; + a[i] = mesh.m_vertices[face.x]; + b[i] = mesh.m_vertices[face.y]; + c[i] = mesh.m_vertices[face.z]; + triangles[i] = i; + } + vertices.toDevice(); + a.toDevice(); + b.toDevice(); + c.toDevice(); + offsets.toDevice(); + triangles.toDevice(); + for (MeshPatchCenterOptions options : + {MeshPatchCenterOptions{}, {0.f, 0.01f, 0.8f}, {2.f, 0.2f, 0.f}, {0.7f, 0.03f, 1.f}}) { + RefreshMeshPatchCenters(offsets.data(), triangles.data(), size_t(3), a.data(), b.data(), c.data(), + centers.data(), vertices.data(), vertices.size(), diagonal.data(), options, nullptr); + DEME_GPU_CALL(cudaDeviceSynchronize()); + centers.toHost(); + compare(std::vector(centers.host(), centers.host() + 3), + mesh.ComputePatchLocations(options.alpha, options.beta, options.gamma), + "CPU/GPU center reduction disagreement"); + } +} + // A cube below an upward-facing open plane produces deep back-face triangle pairings. Fresh centers reject these; // centers left at the cube's previous local position admit them. The open plane cannot mask that mistake with its // own center guard. Compare deformation against loading the target geometry directly, without changing the rule. @@ -113,12 +156,77 @@ float contactForceAfterDeformation(bool deform, bool refresh) { } } // namespace +// A cube resting shallowly on a plane puts each bottom-triangle contact at its old surface centroid. +// The legacy center guard rejects those valid contacts. Offset centers must keep them while tangential motion +// across the plane leaves the summed response unchanged. Per-triangle patches exercise the limiting patch size. +float shallowContactForce(float shift, bool legacy) { + DEMSolver solver(1); + solver.SetVerbosity("ERROR"); + if (legacy) + solver.SetMeshPatchCenterOffsetAlpha(0); + solver.InstructBoxDomainDimension(8, 8, 8); + solver.SetGravitationalAcceleration(make_float3(0)); + solver.SetMeshUniversalContact(true); + solver.SetCDUpdateFreq(1); + solver.SetInitBinNumTarget(8); + solver.SetErrorOutAvgContacts(10000); + solver.SetTimeStepSize(1.e-5); + solver.DefineContactForceModel("force = make_float3(0.f, 0.f, static_cast(overlapDepth));"); + auto material = solver.LoadMaterial({{"E", 1.e7}, {"nu", 0.3}, {"CoR", 0.5}, {"mu", 0.4}}); + auto cube = solver.AddWavefrontMeshObject((GET_DATA_PATH() / "mesh/cube.obj").string(), material); + cube->Scale(0.2f); + cube->SetEachTriangleAsPatch(); + cube->SetMaterial(material); + cube->SetInitPos(make_float3(shift, 0, 0)); + cube->SetFamily(100); + solver.SetFamilyFixed(100); + auto tracker = solver.Track(cube); + auto plane = solver.AddWavefrontMeshObject((GET_DATA_PATH() / "mesh/plane_20by20.obj").string(), material); + // One large triangle contains both cube faces' projections completely. A triangulated plane seam would + // split the contact polygons and move their centroids away from the cube's patch centers, hiding the bug. + plane->m_vertices = {make_float3(-1, -1, 0), make_float3(1, -1, 0), make_float3(0, 1, 0)}; + plane->m_face_v_indices = {make_int3(0, 1, 2)}; + plane->nTri = 1; + plane->SetPatchIDs({0}); + plane->SetInitPos(make_float3(0, 0, -0.096f)); + plane->SetFamily(101); + solver.SetFamilyFixed(101); + solver.Initialize(); + solver.SetTriTriPenetration(0.4f); + solver.DoDynamicsThenSync(1.e-4); + std::vector points, forces; + tracker->GetContactForces(points, forces); + float total = 0; + for (const auto& force : forces) + total += std::abs(force.z); + return total; +} + // Cover shared-node increments, deferred centers, manual policy, lazy readers, and host/device interleaving. int main() try { // Opt-in diagnostics preserve ordinary asynchronous execution for reproducing timing-sensitive failures. deme::SetCudaDebugSyncEnabled(std::getenv("DEME_MESH_DEBUG_SYNC") != nullptr); DEMSolver solver(1); solver.SetVerbosity("ERROR"); + // Exercise non-default settings through initialization, every deformation flavor, and later clump insertion. + solver.SetMeshPatchCenterOffsetAlpha(0.7f); + solver.SetMeshPatchCenterOffsetBeta(0.03f); + solver.SetMeshPatchCenterOffsetGamma(1.f); + using Setter = void (DEMSolver::*)(float); + const Setter setters[] = {&DEMSolver::SetMeshPatchCenterOffsetAlpha, &DEMSolver::SetMeshPatchCenterOffsetBeta, + &DEMSolver::SetMeshPatchCenterOffsetGamma}; + std::cout << "Negative parameter and coordinate checks intentionally emit ERROR diagnostics.\n"; + for (auto setter : setters) { + for (float invalid : {-1.f, std::numeric_limits::infinity(), std::numeric_limits::quiet_NaN()}) { + bool rejected = false; + try { + (solver.*setter)(invalid); + } catch (const std::exception&) { + rejected = true; + } + require(rejected, "invalid solver offset parameter accepted"); + } + } solver.InstructBoxDomainDimension(40, 40, 40); solver.SetGravitationalAcceleration(make_float3(0)); solver.SetTimeStepSize(1.e-5); @@ -127,6 +235,10 @@ int main() try { solver.AddClumps(sphere, std::vector{make_float3(0, 0, 10)}); auto mesh = solver.AddWavefrontMeshObject((GET_DATA_PATH() / "mesh/plane_20by20.obj").string(), material); mesh->Scale(0.1f); + // Multi-patch automatic centers must use initialization's supplied-normal winding correction too. + mesh->UseNormals(true); + for (auto& face : mesh->m_face_v_indices) + std::swap(face.y, face.z); mesh->SetEachTriangleAsPatch(); mesh->SetMaterial(material); mesh->SetFamily(1); @@ -143,7 +255,18 @@ int main() try { single->SetPatchLocations({make_float3(1, 2, 3)}); auto single_tracker = solver.Track(single); solver.Initialize(); + for (auto setter : setters) { + bool rejected = false; + try { + (solver.*setter)(1.f); + } catch (const std::exception&) { + rejected = true; + } + require(rejected, "runtime offset parameter mutation accepted"); + } + checkAutomatic(tracker, mesh->GetCoordsVertices()); const int device = solver.GetGPUDeviceIDs().front(); + checkCenterReduction(device); const size_t count = mesh->GetNumNodes(); require(mesh->GetNumPatches() > 1, "test requires multiple triangle patches"); DeviceInput input(std::max(count, (size_t)mesh->GetNumPatches()), device); @@ -331,6 +454,13 @@ int main() try { require(std::abs(reference_force - refreshed_force) < 1.e-5f, "deformation changed filtered patch forces"); require(std::abs(stale_force - reference_force) > 1.e-3f, "contact regression did not exercise stale-center rejection"); + const float shallow = shallowContactForce(0.f, false); + const float shifted = shallowContactForce(0.2f, false); + const float legacy = shallowContactForce(0.f, true); + std::cout << "Shallow forces: offset=" << shallow << ", shifted=" << shifted << ", legacy=" << legacy << '\n'; + require(shallow > 1.e-3f, "offset centers rejected valid shallow contacts"); + require(std::abs(shallow - shifted) < 1.e-5f, "sliding across the plane changed shallow force"); + require(legacy < shallow * 0.5f, "shallow regression did not distinguish old center policy"); std::cout << "PASS: mesh device deformation, center policies, mixed updates, output, and insertion.\n"; return 0; } catch (const std::exception& e) { diff --git a/src/demo/ModularTests/DEMTest_PatchLocations.cpp b/src/demo/ModularTests/DEMTest_PatchLocations.cpp index f181b2ba..b6bdc43a 100644 --- a/src/demo/ModularTests/DEMTest_PatchLocations.cpp +++ b/src/demo/ModularTests/DEMTest_PatchLocations.cpp @@ -17,6 +17,8 @@ #include #include #include +#include +#include using namespace deme; @@ -25,7 +27,78 @@ bool approxEqual(const float3& a, const float3& b, float tolerance = 1e-5f) { return std::abs(a.x - b.x) < tolerance && std::abs(a.y - b.y) < tolerance && std::abs(a.z - b.z) < tolerance; } +// Analytic geometry checks exercise each term independently, including normals that cannot define an offset. +// A large number of triangles is not needed to distinguish area-weighting from averaging unit normals. +void checkOffsetGeometry() { + auto require = [](bool ok, const char* message) { + if (!ok) + throw std::runtime_error(message); + }; + DEMMesh mesh; + mesh.m_vertices = {make_float3(0, 0, 2), make_float3(1, 0, 2), make_float3(0, 1, 2)}; + mesh.m_face_v_indices = {make_int3(0, 1, 2), make_int3(0, 2, 1)}; + mesh.nTri = 2; + mesh.nPatches = 2; + mesh.m_patch_ids = {0, 1}; + auto centers = mesh.ComputePatchLocations(); + require(approxEqual(centers[0], make_float3(1.f / 3, 1.f / 3, 0.4f)), "exterior offset incorrect"); + require(approxEqual(centers[1], make_float3(1.f / 3, 1.f / 3, 3.6f)), "interior offset incorrect"); + centers = mesh.ComputePatchLocations(0.5f, 0, 1); + require(std::abs(centers[0].z - 1) < 1.e-5f, "custom alpha/gamma ignored"); + centers = mesh.ComputePatchLocations(2, 0.25f, 0); + require(std::abs(centers[0].z - (2 - 0.5f * std::sqrt(2.f))) < 1.e-5f, "custom beta ignored"); + centers = mesh.ComputePatchLocations(0); + require(std::abs(centers[0].z - 2) < 1.e-5f, "alpha=0 did not restore average"); + for (float bad : {-1.f, std::numeric_limits::infinity(), std::numeric_limits::quiet_NaN()}) { + for (int k = 0; k < 3; ++k) { + float args[3] = {1, 0.01f, 0.8f}; + args[k] = bad; + bool rejected = false; + try { + mesh.ComputePatchLocations(args[0], args[1], args[2]); + } catch (const std::exception&) { + rejected = true; + } + require(rejected, "non-finite/negative geometry setting accepted"); + } + } + // Opposing equal-area faces cancel exactly; an empty patch must remain at zero as well. + mesh.m_patch_ids = {0, 0}; + centers = mesh.ComputePatchLocations(); + require(approxEqual(centers[0], make_float3(1.f / 3, 1.f / 3, 2)), "cancelling normals shifted center"); + require(approxEqual(centers[1], make_float3(0, 0, 0)), "empty patch is not zero"); + mesh.m_face_v_indices = {make_int3(0, 0, 0), make_int3(1, 1, 1)}; + mesh.m_patch_ids = {0, 1}; + centers = mesh.ComputePatchLocations(); + require(approxEqual(centers[0], mesh.m_vertices[0]), "degenerate patch shifted center"); + // Equal area opposite normals with a tiny residual should also remain unshifted, independent of units. + mesh.m_vertices.push_back(make_float3(0, 1.0000001f, 2)); + mesh.m_face_v_indices = {make_int3(0, 1, 2), make_int3(0, 3, 1)}; + mesh.m_patch_ids = {0, 0}; + centers = mesh.ComputePatchLocations(); + require(std::abs(centers[0].z - 2) < 1.e-5f, "nearly cancelling normals shifted center"); + // Triangles with +Z area 2 and +X area 1: resultant is proportional to (1,0,2), not (1,0,1). + mesh.m_vertices = {make_float3(0, 0, 0), make_float3(2, 0, 0), make_float3(0, 2, 0), make_float3(0, 0, 1)}; + mesh.m_face_v_indices = {make_int3(0, 1, 2), make_int3(0, 2, 3)}; + centers = mesh.ComputePatchLocations(1, 0, 1); + require(approxEqual(centers[0], make_float3(0.2f, 2.f / 3, -0.1f)), "normal is not area weighted"); + // Normal correction must agree with initialization when users supply opposite vertex winding. + mesh.m_normals = {make_float3(0, 0, -1), make_float3(-1, 0, 0)}; + mesh.m_face_n_indices = {make_int3(0, 0, 0), make_int3(1, 1, 1)}; + mesh.UseNormals(true); + centers = mesh.ComputePatchLocations(1, 0, 1); + require(approxEqual(centers[0], make_float3(7.f / 15, 2.f / 3, 13.f / 30)), "supplied winding ignored"); + std::cout << "PASS: offset geometry, parameters, winding, degenerate and cancelling normals.\n"; +} + int main() { + std::cout << "Negative parameter checks intentionally emit ERROR diagnostics.\n"; + try { + checkOffsetGeometry(); + } catch (const std::exception& e) { + std::cerr << "FAIL: " << e.what() << '\n'; + return 1; + } std::cout << "========================================" << std::endl; std::cout << "DEM Mesh Patch Location Test" << std::endl; std::cout << "========================================" << std::endl; @@ -55,7 +128,7 @@ int main() { test_failures++; } - // Test 2: Multiple patches - should compute centroids + // Test 2: Multiple planar patches use a mesh-size fallback behind their centroid plane. std::cout << "\n--- Test 2: Multi-Patch Mesh Computed Locations ---" << std::endl; auto mesh2 = std::make_shared(); @@ -75,11 +148,11 @@ int main() { auto locations2 = mesh2->ComputePatchLocations(); std::cout << "Number of patches: " << mesh2->GetNumPatches() << std::endl; - // Expected centroids: + // Expected centers: these centroids shifted by -0.01 * sqrt(10) along Z. // Patch 0 triangle centroid: ((0+1+0)/3, (0+0+1)/3, 0) = (0.333, 0.333, 0) // Patch 1 triangle centroid: ((2+3+2)/3, (0+0+1)/3, 0) = (2.333, 0.333, 0) - float3 expected_loc0 = make_float3(1.0f / 3.0f, 1.0f / 3.0f, 0.0f); - float3 expected_loc1 = make_float3(7.0f / 3.0f, 1.0f / 3.0f, 0.0f); + float3 expected_loc0 = make_float3(1.0f / 3.0f, 1.0f / 3.0f, -0.01f * std::sqrt(10.f)); + float3 expected_loc1 = make_float3(7.0f / 3.0f, 1.0f / 3.0f, -0.01f * std::sqrt(10.f)); std::cout << "Computed location for patch 0: (" << locations2[0].x << ", " << locations2[0].y << ", " << locations2[0].z << ")" << std::endl;