Improve mesh patch center derivation

VERSION

@@ -1 +1 @@
-3.0.18
+3.0.19

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,

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) {

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,

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<float3> ComputePatchLocations() const;
+    std::vector<float3> 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 };
 

src/DEM/dT.cpp

@@ -1108,7 +1108,8 @@ void DEMDynamicThread::populateEntityArrays(const std::vector<std::shared_ptr<DE
         if (input_mesh_objs.at(i)->patch_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;
 }

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;

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"))

src/DEM/utils/MeshDeformation.hpp

@@ -4,6 +4,7 @@
 #define DEME_MESH_DEFORMATION_HPP
 
 #include <memory>
+#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<DEMMesh> 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<size_t> patch_offsets;
     DualArray<size_t> patch_triangles;
     DualArray<unsigned int> invalid;
+    // One scalar computed on the same stream before patch reduction; never downloaded during deformation.
+    DualArray<double> diagonal;
 };
 }  // namespace deme
 #endif

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<float3> DEMMesh::ComputePatchLocations() const {
-    std::vector<float3> 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<int> 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<float3>(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<patchID_t>(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<float3> 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<float3> 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<std::array<double, 7>> sums(nPatches, std::array<double, 7>{});
+    std::vector<size_t> 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<patchID_t>(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).

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) ###