A Parallel Contact Detection Method Based on Dual-Domain Activation FDEM

Abstract

This invention discloses a parallel contact detection method for FDEM based on dual-domain activation, comprising: constructing a hash domain and a contact domain; dividing the problem domain into hash units; marking hash units containing TET4 units as active states; and constructing a contact domain centered on the host unit according to spatial neighborhood rules; sorting the TET4 units according to hash keys; traversing the sorted sequence to count the number of units in each hash unit, generating an index interval, marking the active states, and constructing an active unit list; based on the active unit list, sequentially selecting the host hash units, determining the visitor hash unit set according to spatial neighborhood rules, and performing pair detection on the TET4 units in the host and visitor units after activation state filtering; and generating a global contact pair set through two rounds of scanning for subsequent contact mechanics calculations. This invention, through a dual-domain activation mechanism and parallel sorting using GPU, limits the contact search to local active regions, avoiding global traversal and dynamic hash reconstruction.

Description

Technical Field

[0001] This invention belongs to the field of engineering calculation and geomechanical simulation technology, and particularly relates to a parallel contact detection method based on dual-domain activation FDEM. Background Technology

[0002] The Finite-Discrete Element Method (FDEM) is widely used in problems involving large deformation and crack propagation, such as rock fracturing, blasting, cutting, and impact dynamics. As cracks propagate, elements separate, and the block moves, previously adjacent TET4 elements continuously generate new potential contact pairs over multiple time steps. To ensure computational accuracy, the contact state between all TET4 elements needs to be detected at each time step.

[0003] Contact detection is typically divided into coarse contact search (Broad Phase) and fine contact determination (Narrow Phase). The coarse search, used to search for potential TET4 elements from a massive number of elements, is far more computationally intensive than the fine contact determination, making it one of the most computationally demanding and time-consuming parts of FDEM numerical simulations. At the million-element scale, insufficient efficiency of the coarse search can severely slow down the overall computational efficiency of the FDEM solver.

[0004] To improve contact search efficiency, various spatial indexing or neighborhood management methods have been proposed, including bounding volume hierarchy (BVH), sweep-and-Prune (SAP) method, verlet list, uniform mesh partitioning and spatial hashing, etc. However, in the process of three-dimensional FDEM cracking simulation, the computational object often exhibits highly dynamic evolutionary characteristics, and its geometry, topological relationship and local contact state may change significantly between adjacent time steps.

[0005] In this case, the index structure used in existing technologies may need to be frequently updated or reconstructed in practical applications, which may introduce additional computational overhead or cause uneven load distribution during parallel execution, thus adversely affecting the overall execution efficiency. Summary of the Invention

[0006] The purpose of this invention is to provide a parallel contact detection method based on dual-domain activation FDEM. This model constructs an efficient contact detection and collision response calculation mechanism suitable for GPU architecture by parallelizing and reconstructing the contact search and contact determination process, thereby significantly improving the overall computational efficiency of large-scale contact problems while ensuring numerical calculation accuracy and stability.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A parallel contact detection method based on dual-domain activation using FDEM includes the following steps: S1. Based on the discrete model, TET4 units are used to construct a problem domain covering the computation space, and a hash domain is constructed according to the problem domain. Based on the spatial hash grid, hash units containing TET4 units are marked as active to serve as host units. With the host units as the center, a limited number of hash units are selected from the spatial hash grid according to the spatial neighborhood rule to form a contact domain. S2. Obtain the hash key corresponding to each TET4 unit in the hash field, and sort the TET4 units using the hash key as the key, and obtain the hash key sequence. S3. Traverse the sorted hash key sequence, count the number of TET4 units in each hash unit, generate the storage index range corresponding to each hash unit, and mark the hash units containing at least one TET4 unit as active states, and construct a list of active units consisting of all active hash units. S4. Select the active hash units from the list of active units in sequence as host hash units, determine their corresponding visitor hash units according to the spatial neighborhood rules, and filter the visitor hash units based on their active status; perform contact detection on each TET4 unit in the host hash unit and the filtered visitor hash units to obtain contact pairs between the TET4 units in the active visitor hash units and the TET4 units in the host hash units.

[0008] S5. Using a two-round scanning method, generate all candidate contact pairs that can make contact based on the contact pairs, and store them in the global contact pair set to obtain the global contact pair set.

[0009] Furthermore, step S1 is detailed as follows: S101. Construct TET4 elements: Import the geometric model of the rock strata working area or engineering object into the simulation system, and use the mesh generation tool to perform three-dimensional volume meshing of the geometric model to generate a discrete model composed of TET4 elements. S102. Constructing the problem domain: Statistically determine the node coordinates of the TET4 elements to identify the minimum cuboid space covering all TET4 elements; and introduce a preset safety margin outside the space to construct a problem domain for unified management of the effective space for simulation calculations. S103. Hash Unit Construction: Divide the problem domain into cubic hash units of the same size along the three spatial dimensions of x, y, and z. S104. Construction of the hash field: Define the spatial range for limiting the participation in contact retrieval as the hash field, wherein the hash field is a fixed spatial range that covers the problem domain or is equivalent to the problem domain, and its size and structure remain unchanged during the simulation. Mapping between S105, TET4 units, and hash units: This is achieved through the coordinates of the four vertices of each TET4 unit in the Euler coordinate system. ( Calculate the geometric center of the TET4 element. Using the geometric center as the reference point for the spatial location of the TET4 cell, the TET4 cell is assigned to the hash cell, and a hash key is constructed.

[0010] S106. Constructing the contact domain: Based on the spatial hash grid, the activation state of the hash units is determined. Hash units that contain at least one TET4 unit and have contact interactions are marked as active and serve as host units. Centered on a single host unit, 13 hash units selected from the spatial hash grid according to predefined spatial neighborhood rules constitute the contact domain.

[0011] Furthermore, the side length of the hash unit in step S03 for: Based on the equivalent characteristic length of a single TET4 unit Calculate the side length of the hash cell The equivalent feature length is based on the volume of the TET4 unit. Defined as: , It is a proportionality coefficient greater than or equal to 1.

[0012] Furthermore, the hash key construction method in step S105 is as follows: According to the Euler coordinate system, along , , The number of hash units in each direction is expressed as follows: , , For those located at coordinates For the TET4 cell, calculate its hash coordinates based on the side length of the hash cell:

[0013] in, The range of values ​​for are respectively , , The hash key HK is calculated using the following linear mapping: This enables the indexing of hash cells within the hash grid.

[0014] Furthermore, the specific method for obtaining the contact domain in step S106 is as follows: a. Select the research object from the host unit as the host hash unit, and let the hash coordinate of the host hash unit be... ; b. In a spatial hash grid, adjacent hash cells are defined as hash cells whose hash coordinates satisfy the spatial neighborhood rule with the host hash cell; The spatial neighborhood rule is as follows: When a hash unit is at the same layer as the host hash unit, an adjacent hash unit is defined as a hash unit whose index difference with the host hash unit in the x or y direction is ±1, and whose index difference in the other direction is 0, satisfying the following condition: or The hash unit; In the layer adjacent to the host hash unit, the adjacent hash units are defined as having an index difference that satisfies the following condition within the x–y plane projection range: The hash unit satisfies The hash unit; c. The contact domain is formed by the host hash unit itself, the four adjacent hash units in the same layer, and the nine adjacent hash units in the adjacent layers.

[0015] Furthermore, step S2 specifically involves: S201. Each TET4 cell contains basic parameter information, including a unique and unchanging cell identifier TET4ID and a hash key HK related to its spatial location. S202. In parallel, determine whether the geometric center of each TET4 unit is within the range of the hash domain and the contact domain. If so, record its unit identifier and the corresponding hash key; otherwise, mark it as an invalid unit. S203. Using the hash key HK of each valid TET4 unit as the key and its unit identifier TET4ID as the value, perform key-value sorting using the GPU's CUDA parallel sorting interface to generate the sorted hash key sequence Key_sorted and its corresponding value sequence Val_sorted.

[0016] Furthermore, step S3 specifically includes: S301. Traverse the sorted hash key sequence Key_sorted and count the number of TET4 units corresponding to each hash key k: Count[k]. S302. Generate the index range for each hash unit based on the number of TET4 units, Count[k]. The starting position of the index interval beg[k] is less than The count of all hash keys is summed, and the ending position end[k] = beg[k] + Count[k] is obtained. S303. When Count[k] > 0, the hash unit is marked as active; when Count[k] = 0, its index range is empty. S304. Based on the activation state of the hash units, rearrange all hash units, store all activated hash units consecutively at the beginning of the data structure, and place inactive hash units at the end to form a list of activated units.

[0017] Furthermore, step S4 specifically involves: S401. Select the hash units in the active state from the list of active units in turn as the host hash units; S402. For the current host hash unit, determine its corresponding visitor hash unit according to the spatial neighborhood rule. The visitor hash unit set includes the host unit itself, 4 adjacent units in the same layer, and 9 adjacent units in the upper and lower layers, for a total of 13 hash units. S403. Before performing contact detection, filter the visitor hash unit set based on the activation state array to obtain visitor hash units with an activation state of 1. S404. Perform contact detection on each TET4 unit in the host hash unit and the filtered visitor hash unit to obtain contact pairs.

[0018] Furthermore, the specific method for obtaining the contact pair in step S404 is as follows: Contact detection is performed on all TET4 units within the host hash unit except itself, and contact detection is performed on each TET4 unit contained in all visitor hash units to form contact pairs.

[0019] Furthermore, step S5 specifically involves: S501, First round of scanning phase: Parallel traversal of each active hash unit and its corresponding visitor hash unit, counting of unit pairs that meet the contact conditions, generating the number of potential contact pairs corresponding to each active hash unit, and calculating the prefix and index information of each active hash unit based on the number of contact pairs, where the prefix and index information is used to determine the write start position of each active hash unit in the global contact pairs. S502, Second round of scanning phase: Based on the prefix and index information of the hash unit of each active state, a non-overlapping global contact pair set writing interval is allocated for each active state hash unit. Each computing thread writes all the identified potential contact pairs into the global contact pair set in parallel within its corresponding writing interval. S503. The global contact pair set is used as input data to drive subsequent contact force calculations or collision response solutions to update the motion state of the TET4 unit.

[0020] The beneficial effects of this invention are as follows: (1) Significantly improve contact calculation efficiency: This invention maps the contact search and collision determination process to the GPU parallel computing platform and utilizes the characteristics of GPU large-scale thread concurrent execution to achieve parallel processing of multi-body contact relationships. Compared with the traditional contact model based on CPU serial or limited parallel methods, it can significantly reduce the time complexity of contact calculation and improve the overall calculation efficiency.

[0021] (2) Enhanced adaptability to large-scale contact problems: The contact collision model of the present invention is designed in parallel for the contact characteristics of high-density, multi-body systems, which can effectively support contact detection and collision response calculation of large-scale particle or block systems, avoiding the problem of significant decrease in computational efficiency when the number of contacts increases sharply in the existing technology, and has good scalability.

[0022] (3) Fully leverage the advantages of GPU hardware computing resources: The algorithm structure of this invention is matched with the GPU thread model and storage access characteristics. By parallelizing the contact search and collision determination process, unnecessary serial operations and thread divergences are reduced, improving the utilization rate of GPU computing resources, thereby achieving high-performance numerical simulation.

[0023] (4) Ensuring computational accuracy and numerical stability: While improving computational efficiency, the contact collision model of the present invention maintains the same physical constraints as the traditional contact model in the process of contact determination and mechanical response calculation, ensuring the accuracy of contact force calculation results and the stability of numerical simulation process, and is suitable for engineering analysis and scientific computing scenarios. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of hash unit partitioning provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the defined contact area provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of spatial hash partitioning and TET4 cell mapping provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of the host unit and the visitor unit provided in an embodiment of the present invention; Figure 5 This is a flowchart of the contact detection process provided in an embodiment of the present invention; Figure 6 This is a flowchart of parallel contact detection with dual-domain activation of hash domain and contact domain provided in an embodiment of the present invention. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0026] The application principle of the present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0027] Example 1: A parallel contact detection method based on dual-domain activation using FDEM includes the following steps: S1. Construct the hash field and contact field.

[0028] S101, Constructing the TET4 unit First, the geometric model of the rock strata working area or engineering object is imported into the simulation system. Then, a three-dimensional volume mesh is generated from the geometric model using a mesh generation tool (such as Gmsh), producing a discrete model composed of multiple four-node tetrahedral elements (TET4 elements). The TET4 element is the smallest mechanical and geometric discrete element of the simulation object. Each TET4 element is uniquely determined by its node coordinates, and its spatial position and volume parameters can be directly obtained from the mesh data of the three-dimensional volume mesh.

[0029] S102. Constructing the Problem Domain The problem domain is the spatial computational region used for numerical solution and contact retrieval. It defines the effective spatial range of the simulation computation and serves as the reference boundary for spatial hash grid construction and coordinate mapping. The problem domain is constructed during the simulation initialization phase, as follows:

[0030] After obtaining all TET4 elements, spatial statistics are performed on their node coordinates. The minimum and maximum coordinate values ​​in the x, y, and z spatial directions are calculated in the Euler coordinate system to determine the minimum cuboid spatial range covering all TET4 elements. A preset safety margin is introduced outside this spatial range to construct the spatial management boundary for numerical updates and contact retrieval, i.e., the problem domain. Therefore, the problem domain is obtained by inverse calculation based on the spatial distribution range of the TET4 elements and is used to uniformly manage the effective space for simulation calculations.

[0031] S103, Hash Unit Construction During the simulation, when a TET4 element flies out of the hash domain used for contact retrieval due to breakage, collision, or inertial motion, it is still allowed to continue to update its velocity and position kinematics within the problem domain. When the TET4 element goes further out of the problem domain boundary, it is marked as a delocalized element, and its subsequent calculations are frozen. That is, for delocalized TET4 elements that have gone out of the problem domain, the updating of their velocity, position, and other kinematic parameters is stopped, and they are no longer included in the scope of contact detection and contact force calculation, so as to ensure the stability of the simulation calculation process and control the calculation scale.

[0032] See Figure 1 In the Euler coordinate system, the problem domain is divided into cubic hash cells of equal size along the x, y, and z spatial dimensions. The side length of the hash cell in each dimension is determined based on the geometric characteristics of the TET4 cell. Specifically, it is based on the volume of a single TET4 cell. Define its equivalent feature length ,in

[0033] The equivalent characteristic length is used to characterize the average spatial interaction scale of TET4 cells when contact occurs. Based on this, the side length of the hash cell is set. It satisfies a proportional relationship with the equivalent feature length:

[0034] in, A scaling factor greater than or equal to 1 is used to ensure that any TET4 cells that may come into contact are spatially mapped to the same or adjacent hash cells, thereby ensuring the completeness of the contact search.

[0035] S104. Construction of Hash Domain In the above description, for ease of explanation, the hash domain is implicitly set to be consistent with the problem domain. However, in practical applications, the spatial range of the hash domain does not necessarily cover the entire problem domain; its size can be selected according to the specific research object and simulation requirements. The hash domain is mainly used to limit the spatial range of participants in contact retrieval, and its setting principle is to cover areas where potential contact behaviors are relatively concentrated.

[0036] In embodiments of the present invention, the hash domain is preferably set to be consistent with the physical size of the object under study, so that the spatial range of the hash domain covers the main computational object involved in the contact and fracturing behavior. For example, in a simulation scenario of PDC drill bit cutting rock, the object under study is usually the rock sample body to be fractured, so the overall volume of the rock sample can be used as the spatial range of the hash domain to construct a regular spatial hash grid.

[0037] Based on this, the hash field can also be set to cover the entire problem domain or be equivalent to the problem domain, so as to avoid frequently adjusting the hash field range or rebuilding the hash structure during the simulation process. This helps to maintain the consistency of the space partitioning method and simplify the organization and implementation of the calculation process.

[0038] By spatially covering the problem domain with a hash field and maintaining its fixed size and structure throughout the simulation, the determinism and predictability of hash index computation can be guaranteed, thus adapting to the implementation requirements of parallel contact detection. Furthermore, the hash field covering the problem domain does not mean that all hash cells participate in contact computation; in the actual computation process, contact detection is performed only on hash cells containing active TET4 cells, while empty hash cells located far from the contact region do not introduce additional computational burden.

[0039] Taking the contact between the cutting tooth and the rock interface as an example, the contact behavior is mainly concentrated in a local area near the cutting tooth and the rock surface, while rock units far from this area have little impact on the contact results. Even if the material undergoes severe fragmentation and some rock fragments detach from the main contact area, as long as they remain within the hash domain coverage, they can be managed through a unified spatial hash index. When the fragments further move away from the physical range corresponding to the research object, they can be excluded from subsequent contact detection according to preset rules. In this way, efficient calculation of local contact behavior is achieved while avoiding dynamic hash reconstruction.

[0040] Mapping of S105, TET4 units and hash units To facilitate the explanation of the specific process of hash mapping and index construction in this embodiment, Figure 6 The diagram illustrates the correspondence between some TET4 units and hash units. It should be noted that this correspondence is not arbitrarily set, but determined according to a unified space partitioning rule and mapping process.

[0041] In this embodiment, for each TET4 element participating in the calculation, the coordinates of the four vertices of the TET4 element in the Euler coordinate system are first used. ( Calculate the geometric center of the TET4 element. The formula for calculating the geometric center is as follows: The geometric center is a unique reference point used to characterize the spatial location of the TET4 cell. Based on the location of the reference point in space, the TET4 cell is assigned to a hash cell. Through the above assignment process, it can be guaranteed that any TET4 cell corresponds to only one specific hash cell at any given time.

[0042] After determining the hash cell to which the TET4 cell belongs, a unique hash key is assigned to the TET4 cell based on its spatial index position in the hash field. The hash key is used to identify the position of the TET4 cell in the spatial hash structure and serves as the index for subsequent contact pair generation and collision detection processes.

[0043] It should be understood that the TET4 unit numbers, hash key values, and their arrangement shown in the figure are only for illustrating the basic process of hash mapping and index construction, and do not constitute a limitation on the specific numbering rules, storage order, data structure form, or implementation method of this invention. In practical applications, the specific storage method, index order, and data organization form of the relevant mapping relationships can all be equivalently adjusted according to the model size, computing platform, and parallel implementation requirements without changing the above mapping logic. In the Euler coordinate system, along , , The number of hash units in each direction is expressed as follows: , , For those located at coordinates For the TET4 cell, calculate its hash coordinates based on the side length of the hash cell:

[0044] in, The range of values ​​for are respectively , , The hash key HK is further calculated using the following linear mapping relationship:

[0045] This enables a unique index for each hash unit within the hash grid.

[0046] S106, Constructing the contact domain (1) In this invention, the activation state of hash units is first determined based on the spatial hash grid, and hash units that contain at least one TET4 unit and have the possibility of contact interaction are marked as active. The active hash units are used as host units in the contact detection process to trigger subsequent local contact search.

[0047] (2) such as Figure 4 As shown, a contact domain refers to a local region centered on a single host cell and composed of a finite number of hash cells selected from a spatial hash grid according to predefined spatial neighborhood rules. The spatial neighborhood rules are used to limit the search range for contact detection, ensuring that the calculation of contact interactions is performed only within the host cell and its adjacent hash cells, as detailed below:

[0048] a. Select the research object from the host unit as the host hash unit, and let the hash coordinate of the host hash unit be... .

[0049] b. In a spatial hash grid, adjacent hash cells are defined as hash cells whose hash coordinates satisfy the spatial neighborhood rule with the host hash cell.

[0050] The spatial neighborhood rule is as follows: When a hash unit is at the same layer as the host hash unit, adjacent hash units are defined as hash units whose index difference with the host hash unit in the x or y direction is ±1, and whose index difference in the other direction is 0, i.e., satisfying... or The hash units are thus excluded from the diagonal direction of the same layer.

[0051] In the layer adjacent to the host hash unit, the adjacent hash units are defined as having an index difference that satisfies the following condition within the x–y plane projection range: The hash unit that satisfies The hash unit.

[0052] c. Thus, a single contact domain consists of the host hash unit itself, four adjacent hash units in the same layer, and nine adjacent hash units in adjacent layers, totaling 13 hash units. The aforementioned adjacency relationships are entirely determined by the discrete offset rule (spatial neighborhood rule) of the hash index, and are independent of specific geometric distances.

[0053] Unlike hash domains, contact domains are dynamic, with their spatial locations and constituent units updating as the activated hash units change during the simulation. For example, under external loads or during tool movement, as the position of the PDC cutting teeth changes, different hash units will be activated sequentially and become new host units, thus forming new contact domains. Furthermore, multiple independent contact domains can exist simultaneously in the same simulation, each used for contact detection calculations within different local areas. In contrast, hash domains, as a unified spatial partitioning structure, maintain a relatively fixed spatial extent and mesh topology during the simulation.

[0054] S2. Sort the TET4 units according to the hash key.

[0055] (1) A data organization structure is constructed using TET4 cells as the basic objects, and its length is consistent with the total number of TET4 cells participating in the calculation in the hash domain. For each TET4 cell, there is a set of basic parameter information, including the cell identifier and the hash mapping structure related to its spatial location. The cell identifier is used to uniquely identify the TET4 cell and remains unchanged throughout the simulation process to distinguish different TET4 cells; the cell identifier can be denoted as TET4ID. The hash mapping structure related to the spatial location includes its corresponding hash key. The hash key is used to characterize the spatial belonging relationship of the TET4 cell in the hash domain and can be dynamically updated as the TET4 cell moves in Eulerian space.

[0056] Furthermore, the basic parameter information associated with each TET4 element may also include the element's geometric center position, vertex coordinates, velocity information, and other state parameters used for mechanical calculations or contact determination. By separating the stable element identifier from the dynamically updated spatial mapping information, it is beneficial to efficiently manage and access TET4 elements in parallel during simulation.

[0057] (2) Subsequently, for each TET4 cell, it is determined in parallel whether it is within the predefined hash domain and contact domain. The geometric center position of the TET4 cell in the Euler coordinate system is used as the spatial determination criterion. When the geometric center is within the spatial range of the hash domain or contact domain, the corresponding cell identifier and the hash key obtained by mapping are recorded and marked as a valid cell. For TET4 cells whose geometric center does not fall within the range of the hash domain and contact domain, they are uniformly marked as invalid cells and are automatically ignored in subsequent steps such as sorting and contact detection, and do not participate in the calculation.

[0058] (3) Then, using the hash key corresponding to each valid TET4 unit Using the CUDA (Compute Unified Device Architecture, a parallel computing platform and programming model from NVIDIA that leverages the parallel computing capabilities of graphics processing units (GPUs) to accelerate general-purpose computing, the GPU's CUDA parallel sorting interface is used to perform key-value sorting on all TET4 cells. This allows for key-value sorting of millions of TET4 cells, adapting to the needs of large-scale data sorting. During the sorting process, the cell identifier of the TET4 cell (TET4ID or its index in the cell array) is used as the value, based on the hash key. The TET4 cells are sorted by comparing their numerical values, regardless of the cell number. For example... Figure 6As shown in Figure B, after sorting, TET4 units with the same hash key are distributed continuously in memory, thereby aggregating local units that fall into the same hash unit in space in the data structure, providing a basis for subsequent enumeration of neighborhood candidate units and generation of contact pairs.

[0059] In its implementation, the sorting method includes the following steps: The hash key of each TET4 cell is computed in parallel. Generate a key array Key[i] = AND value array Value[i] = ; The GPU parallel sorting operator is invoked to sort the key array Key in ascending order and simultaneously rearrange the values ​​to obtain the sorted Key_sorted and Val_sorted. Val_sorted represents the sequence of TET4 unit identifiers corresponding to each position after sorting, and Key_sorted represents the sequence of hash keys that correspond to each of them.

[0060] The sorting operation described above allows TET4 units corresponding to the same hash unit to form a contiguous storage segment. This enables the candidate unit set within the hash unit to be directly obtained through a single index location in subsequent contact detection, avoiding global traversal, improving neighborhood search efficiency, and enhancing the continuity and consistency of GPU parallel memory access.

[0061] S3, Hash Unit Activation Management.

[0062] Traverse the sorted hash key sequence, count the number of TET4s in each hash unit, as shown in C of Figure 6, and generate the index range corresponding to each hash unit. Specifically, for each hash key (Corresponding hash unit), let Count[k] be the number of TET4 units it contains; then its index range starting position beg[k] is less than The result is obtained by summing the counts of all hash keys, and the ending position is end[k] = beg[k] + Count[k].

[0063] When Count[k] > 0, it corresponds to a non-zero continuous storage interval in the sorted sequence, indicating that at least one TET4 unit is mapped to this hash unit. This type of hash unit has practical significance in subsequent contact detection, so it is marked as an active state. The hash units of the active states are counted to construct an active state array. When Count[k] = 0, beg[k] = end[k], and its index interval degenerates into an empty interval. It does not participate in any contact candidate unit enumeration in subsequent calculations, thus avoiding unnecessary computational overhead.

[0064] Therefore, the storage range of all TET4 units within each hash unit in the sorted sequence can be determined by... The only certainty.

[0065] like Figure 6 As shown in Figure D, the hash units are rearranged according to their activation states. The original hash unit activation state array is traversed, and the activation state identifier corresponding to each hash key is collected. An auxiliary list is constructed based on the activation state identifiers, where the auxiliary list contains all hash keys with an activation state of 1. The hash keys in the auxiliary list are sorted in ascending order to form an ordered sequence of activation keys. This ordered sequence of activation keys is then sequentially filled into the starting position of the rearranged data structure. Hash keys with an activation state of 0 from the original hash unit activation array are appended to subsequent positions in the rearranged data structure in their original order, so that activated hash units are continuously clustered in the array, while inactive hash units are placed at the end, forming an activation unit list.

[0066] S4. Neighborhood contact detection based on activation domain.

[0067] From the list of active units constructed in step S3, hash units in active states are selected sequentially as the host hash units for the current contact detection. Specifically, by counting the number of TET4 units within a hash unit and generating an index range, an active state array Active(k) can be obtained. When Active(k) = 1, it indicates that the k-th hash unit contains at least one TET4 unit and has the potential to undergo contact interaction.

[0068] For each host hash unit, a set of visitor hash units is constructed in its spatial neighborhood according to the spatial neighborhood rules in step S106. This set includes four hash units adjacent to each other in the same layer along the x and y directions, and nine hash units in the lower layer corresponding to the 3×3 grid centered on the projection of the host hash unit onto the x–y plane, thus forming a total of 13 hash units. These 13 hash units are then used as visitor hash units.

[0069] Before performing contact detection, the visitor hash cells are first filtered based on the active state array. Only when the visitor hash cell satisfies Active(k)=1, that is, it contains at least one TET4 cell, will the index range corresponding to the hash cell be accessed and contact detection of the TET4 cell within it be performed.

[0070] Contact detection is performed on each TET4 unit in both the host hash unit and the filtered visitor hash units. Two types of contact detection operations are executed: First, contact detection is performed on all other TET4 units within the host hash unit except itself, meaning TET4 units within the same hash unit are considered contact pairs by default. Second, contact detection is performed on each TET4 unit contained in all visitor hash units, meaning that TET4 units in all active visitor hash units are considered contact pairs with TET4 units in the host hash unit by default. For inactive visitor hash units, contact detection is not performed on their TET4 units, thus avoiding invalid computation of spatial regions where contact interaction is impossible in the current time step. This activation state filtering mechanism significantly reduces the number of units participating in computation while maintaining the integrity of contact detection.

[0071] S5. Two rounds of scanning are used to construct a contact candidate list.

[0072] In the first round of scanning, only the number of potential contact pairs within each active visitor hash unit and host hash unit is counted; no actual data write operations are performed. Specifically, each computation thread counts the units (active visitor hash units and host hash units) that meet the coarse contact conditions for the corresponding active visitor hash unit and host hash unit and their neighborhood range, generating statistical results reflecting the number of potential contact pairs for each hash unit. Based on these statistical results, the prefix and index information of each active hash unit is further calculated to determine the write start position and interval length of each hash unit in the global contact pair storage structure.

[0073] In the second round of scanning, based on the prefix and index information generated in the first round, a unique and non-overlapping write interval is assigned to each active hash unit. Subsequently, each computation thread writes all potential contact pairs identified during the contact detection process sequentially into the global contact pair set only within its corresponding write interval. Since the write intervals have been logically divided, there is no competition for the same storage location between different threads, thus enabling parallel writing without the need for synchronization or atomic operations.

[0074] After the second round of scanning is completed, the generated global contact pair set is used to characterize all TET4 element pairs that may have contact interactions within the current time step. The global contact pair set serves as input data for subsequent physics calculation modules, driving contact force calculation, collision response solving, and element motion state updates.

[0075] In practical applications, the global set of contact pairs can be input into the contact force model or constraint solution model to calculate the normal force, tangential force or contact constraint conditions between each contact pair, and update the velocity, displacement and force state of the corresponding TET4 element accordingly, thereby ensuring the accuracy and stability of the physical simulation process when contact behavior occurs.

[0076] In the above manner, the global contact pair set constructed by this invention provides an effective candidate unit pair set for the calculation of contact-related physical quantities, enabling contact calculation to be limited to the unit range related to contact behavior, thereby reducing the participation of irrelevant units in the calculation, which is beneficial to improving the efficiency of simulation calculation and enhancing the stability of the calculation process and the rationality of the results.

[0077] Example 2: (1) Structural comparison This comparison, independent of specific computational examples or experimental time, is used to illustrate the inherent differences between different contact detection methods in a parallel environment from the perspective of algorithm organization and GPU execution mechanism, with a focus on control flow consistency, memory access patterns, and concurrent write behavior.

[0078] Table 1. Comparison of Spatial Organization Characteristics of Different Contact Detection Methods As shown in Table 1 above, different contact detection methods have fundamental differences in spatial organization and activation strategies. These differences directly affect their execution characteristics and scalability in parallel computing environments, especially on GPU platforms.

[0079] The LBVH method, based on a hierarchical bounding volume structure, achieves contact clipping by constructing a hierarchical spatial bounding volume, effectively reducing invalid detections in scenes with relatively stable geometric distributions. However, the activation process of this method relies on tree structure traversal, and whether spatial units participate in the computation is implicitly determined by the traversal path. When the geometry changes drastically, the tree structure needs frequent reconstruction, making it difficult to maintain stable activation behavior, thus affecting execution consistency in strongly fragmented FDEM simulations.

[0080] The neighborhood list approach primarily relies on the sorting of objects' projection intervals on one or more coordinate axes to identify potential contact objects. Its effectiveness is based on the assumption that the relative order of objects changes only slightly between adjacent time steps. During the fission process, a large number of new units are generated, accompanied by rapid changes in their positions and sorting relationships. This necessitates frequent updates to the sorting results and neighborhood information, resulting in a global activation granularity that makes it difficult to form a stable local computational scope.

[0081] Uniform grid or spatial hashing methods limit the contact search range through regular spatial discretization. Their spatial organization is relatively simple and does not rely on the continuity of geometric topology, thus adapting to a certain degree of dynamic change. However, these methods typically employ a global activation strategy, meaning that all spatial units participate in the contact detection process at every time step, even if a large number of units do not have potential contact objects at the current stage. Structurally, this makes it difficult to avoid the introduction of invalid computations.

[0082] In contrast, the dual-domain activation FDEM contact detection method proposed in this invention explicitly distinguishes the global spatial index structure (hash domain) from the local computational regions (contact domains) where actual contact interactions occur at the spatial organization level. Through dynamic management of the contact domain, only spatial regions containing potential contact units are marked as active, thus spatially confining the contact detection process to locally active regions. This activation mechanism does not rely on the temporal continuity of geometric structures or object ordering relationships, and can decouple from the topological evolution during the fission process, thereby maintaining good structural stability under strong fission scenarios.

[0083] Therefore, although the above methods can all perform coarse contact search at the functional level, their structural differences in spatial activation methods determine that their applicability and execution characteristics in three-dimensional FDEM strong dynamic cracking simulation are significantly different. This invention, by introducing an explicit, local, and dynamically updatable activation mechanism, provides a more stable structural foundation for subsequent parallel contact detection and candidate pair construction.

[0084] (0) Comparison of parallel execution mechanisms Table 2 Comparison of Parallel Execution Mechanisms of Different Contact Detection Methods As shown in Table 2 above, different contact detection methods have significant differences in parallel execution mechanisms. These differences are mainly reflected in concurrent write methods, thread control flow consistency, and data layout methods, which directly affect their execution efficiency and scalability on the GPU platform.

[0085] For the LBVH method based on the linear bounding volume hierarchy, the traversal process is usually implemented non-recursively, avoiding explicit function recursive calls. However, in actual traversal, the node paths visited by different threads are greatly affected by the geometric distribution, resulting in significant differences in execution paths between threads and easily leading to branch divergence. Furthermore, the node access order is determined by the tree structure, making it difficult to guarantee continuous data access patterns during parallel execution.

[0086] The neighborhood list approach primarily relies on sorting and interval scanning operations to filter potential contact objects. While this method exhibits high execution efficiency when object distribution is relatively stable, in parallel environments, the number of comparisons and neighborhood update frequencies differ between different objects, leading to an imbalance in thread execution paths and computational load. In scenarios involving fragmentation or drastic rearrangement, this imbalance is further exacerbated, negatively impacting control flow consistency.

[0087] In parallel implementations, uniform grid or spatial hashing methods typically construct data structures by appending object indices or contact pairs to spatial cells. When multiple threads simultaneously write data to the same spatial cell, atomic operations are often required to ensure data consistency. When objects are densely packed or contacts are highly concentrated within a spatial cell, this type of concurrent write behavior can easily generate additional synchronization overhead, thus impacting overall parallel efficiency. Furthermore, the data distribution between spatial cells is often quite discrete, which is detrimental to contiguous memory access.

[0088] In contrast, the FDEM parallel contact detection method based on dual-domain activation proposed in this invention explicitly avoids concurrent write behavior in its parallel execution mechanism. By pre-allocating non-overlapping write intervals for different spatial units during the contact pair construction phase, each thread performs sequential write operations only within its corresponding interval, structurally avoiding contention for global counters or shared memory locations. Simultaneously, the contact detection process employs a consistent execution path within the activated hash unit range, ensuring a high degree of consistency in thread control flow. Combined with a continuous data layout, this improves GPU pipeline execution efficiency.

[0089] Therefore, the differences in parallel execution mechanisms among various methods do not stem from the choice of specific implementation details, but are inherently determined by their algorithm organization and data structure design. This invention, through systematic design in concurrent write control and data layout, enables higher execution consistency and stability in a GPU parallel environment, providing a favorable parallel computing foundation for large-scale 3D FDEM contact detection.

[0090] (3) Scale reduction ratio To facilitate the analysis of the differences in computational scale among different contact detection methods, a unified symbol definition is first introduced: For the number of global discrete units, The number of units participating in the coarse search, The number of spatial units (grids or hashes) in the entire domain. To activate the number of spatial units, This represents the number of candidate contact pairs.

[0091] Table 3 Comparison of computational scale reduction for different contact detection methods Based on the above symbol definitions, a comparative analysis of the participation scale of different contact detection methods in the coarse contact search stage is conducted, as shown in Table 3. Different methods show significant differences in the participation unit scale and spatial unit scale.

[0092] For the LBVH method based on bounding volume hierarchy, its coarse contact search process relies on tree structure traversal. Whether a cell participates in the search is dynamically determined by the traversal path, and is greatly affected by geometric distribution and hierarchical pruning effects. When the geometry is relatively stable, some spatial regions can be effectively pruned; however, during the 3D FDEM cracking process, as cracks propagate and cells separate, the tree structure changes frequently, and the traversal path fluctuates significantly with time steps. This makes it difficult to keep the size of the cells actually participating in the search stable, resulting in a participation characteristic that varies with the traversal path.

[0093] The neighborhood list approach primarily filters potential contact objects based on sorting the objects' projection intervals along the coordinate axes. This method implicitly assumes that the relative order between objects changes only slightly between adjacent time steps. In strongly fragmented scenarios, a large number of new units are generated, accompanied by significant changes in spatial location. The sorting results and neighborhood information need frequent updates, causing the vast majority of units to participate in the computation during the coarse search phase. The number of participating units is statistically close to the global unit count. It is difficult to achieve effective scaling down.

[0094] Uniform grid or spatial hashing methods define the search range by discretizing a regular space. They do not rely on the continuity of geometric topology and therefore can adapt to a certain degree of dynamic change. However, these methods typically employ a global activation strategy, meaning that at each time step, all spatial units are considered potential search areas, corresponding to a certain number of participating spatial units. Approximately equal to the total number of spatial units in the entire domain The number of units participating in the search simultaneously Also close to In scenarios where the element distribution is highly uneven or the contact is highly localized, this structural characteristic makes it difficult to avoid a large amount of invalid computation.

[0095] In contrast, the FDEM parallel contact detection method based on dual-domain activation proposed in this invention introduces a direct constraint on the scale of participation by explicitly distinguishing between the global spatial index structure and the local computational regions where actual contact interactions occur. On the one hand, it includes only units located within the contact domain in the coarse contact search, thus limiting the scale of participating units. Significantly smaller than the number of global units On the other hand, only spatial units containing potential contact units are activated, thus increasing the number of activated spatial units. Much smaller than the total number of spatial units in the entire domain This scale constraint mechanism does not depend on the temporal continuity of geometric structure or object ordering relationships, and can remain relatively stable during the fission evolution process.

[0096] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A parallel contact detection method based on dual-domain activation using FDEM, characterized in that, Includes the following steps: S1. Construct a problem domain covering the computation space using TET4 units based on the discrete model, and construct a hash domain based on the problem domain; Based on the spatial hash grid, hash units containing TET4 units are marked as active to serve as host units, and a limited number of hash units are selected from the spatial hash grid to form a contact domain, centered on the host unit and according to the spatial neighborhood rule. S2. Obtain the hash key corresponding to each TET4 unit in the hash field, and sort the TET4 units using the hash key as the key, and obtain the hash key sequence. S3. Traverse the sorted hash key sequence, count the number of TET4 units in each hash unit, generate the storage index range corresponding to each hash unit, and mark the hash units containing at least one TET4 unit as active states, and construct a list of active units consisting of all active hash units. S4. Select the active hash units from the list of active units in sequence as host hash units, determine the corresponding visitor hash units according to the spatial neighborhood rules, and filter the visitor hash units based on their active status; perform contact detection on each TET4 unit in the host hash unit and the filtered visitor hash units to obtain contact pairs between the TET4 units in the active visitor hash units and the TET4 units in the host hash units. S5. Using a two-round scanning method, generate all candidate contact pairs that can make contact based on the contact pairs, and store them in the global contact pair set to obtain the global contact pair set.

2. The FDEM parallel contact detection method based on dual-domain activation according to claim 1, characterized in that, Step S1 is as follows: S101. Construct TET4 elements: Import the geometric model of the rock strata working area or engineering object into the simulation system, and use the mesh generation tool to perform three-dimensional volume meshing of the geometric model to generate a discrete model composed of TET4 elements. S102. Constructing the problem domain: Statistically determine the node coordinates of the TET4 elements to identify the minimum cuboid space covering all TET4 elements; and introduce a preset safety margin outside the space to construct a problem domain for unified management of the effective space for simulation calculations. S103. Hash Unit Construction: Divide the problem domain into cubic hash units of the same size along the three spatial dimensions of x, y, and z. S104. Construction of the hash field: Define the spatial range for limiting the participation in contact retrieval as the hash field, wherein the hash field is a fixed spatial range that covers the problem domain or is equivalent to the problem domain, and its size and structure remain unchanged during the simulation. Mapping between S105, TET4 units, and hash units: This is achieved through the coordinates of the four vertices of each TET4 unit in the Euler coordinate system. ( Calculate the geometric center of the TET4 element. Using the geometric center as the reference point for the spatial location of the TET4 cell, the TET4 cell is assigned to a hash cell, and a hash key is constructed. S106. Constructing a contact domain: Based on the spatial hash grid, the activation state of the hash units is determined. Hash units that contain at least one TET4 unit and have contact interactions are marked as active and serve as host units. Centered on a single host unit, 13 hash units selected from the spatial hash grid according to predefined spatial neighborhood rules constitute a contact domain.

3. The FDEM parallel contact detection method based on dual-domain activation according to claim 2, characterized in that, The side length of the hash cell in step S103 for: Based on the equivalent characteristic length of a single TET4 unit Calculate the side length of the hash cell The equivalent feature length is based on the volume of the TET4 unit. Defined as: It is a proportionality coefficient greater than or equal to 1.

4. The FDEM parallel contact detection method based on dual-domain activation according to claim 2, characterized in that, The hash key construction method in step S105 is as follows: According to the Euler coordinate system, along , , The number of hash units in each direction is expressed as follows: , , For those located at coordinates For the TET4 cell, calculate its hash coordinates based on the side length of the hash cell: in, The range of values ​​for are respectively , , ; and calculate the hash key using a linear mapping relationship: This enables the indexing of hash cells within the hash grid.

5. The FDEM parallel contact detection method based on dual-domain activation according to claim 2, characterized in that, The specific method for obtaining the contact domain in step S106 is as follows: a. Select the research object from the host unit as the host hash unit, and let the hash coordinate of the host hash unit be... ; b. In a spatial hash grid, adjacent hash cells are defined as hash cells whose hash coordinates satisfy the spatial neighborhood rule with the host hash cell; The spatial neighborhood rule is as follows: When a hash unit is at the same layer as the host hash unit, an adjacent hash unit is defined as a hash unit whose index difference with the host hash unit in the x or y direction is ±1, and whose index difference in the other direction is 0, satisfying the following condition: or The hash unit; In the layer adjacent to the host hash unit, the adjacent hash units are defined as having an index difference that satisfies the following condition within the x–y plane projection range: The hash unit satisfies The hash unit; c. The contact domain is formed by the host hash unit itself, the four adjacent hash units in the same layer, and the nine adjacent hash units in the adjacent layers.

6. The FDEM parallel contact detection method based on dual-domain activation according to claim 1, characterized in that, Step S2 is as follows: S201. Each TET4 cell contains basic parameter information, including a unique and unchanging cell identifier TET4ID, and a hash key related to its spatial location. S202. In parallel, determine whether the geometric center of each TET4 unit is within the range of the hash domain and the contact domain. If so, record its unit identifier and the corresponding hash key; otherwise, mark it as an invalid unit. S203. Using the hash key of each valid TET4 unit as the key and its unit identifier TET4ID as the value, perform key-value sorting using the GPU's CUDA parallel sorting interface to generate the sorted hash key sequence Key_sorted and its corresponding value sequence Val_sorted.

7. The FDEM parallel contact detection method based on dual-domain activation according to claim 1, characterized in that, Step S3 is as follows: S301. Traverse the sorted hash key sequence Key_sorted and count the number of TET4 units corresponding to each hash key Count[k]. S302. Generate the index range for each hash unit based on the number of TET4 units, Count[k]. The starting position of the index interval beg[k] is less than The count of all hash keys is summed, and the ending position end[k] = beg[k] + Count[k] is obtained. S303. When Count[k] > 0, the hash unit is marked as active; when Count[k] = 0, its index range is empty. S304. Based on the activation state of the hash units, rearrange all hash units, store all activated hash units consecutively at the beginning of the data structure, and place inactive hash units at the end to form a list of activated units.

8. The FDEM parallel contact detection method based on dual-domain activation according to claim 1, characterized in that, Step S4 is as follows: S401. Select the hash units in the active state from the list of active units in turn as the host hash units; S402. For the current host hash unit, determine its corresponding visitor hash unit according to the spatial neighborhood rule. The visitor hash unit set includes the host unit itself, 4 adjacent units in the same layer, and 9 adjacent units in the upper and lower layers, for a total of 13 hash units. S403. Before performing contact detection, filter the visitor hash unit set based on the activation state array to obtain visitor hash units with an activation state of 1. S404. Perform contact detection on each TET4 unit in the host hash unit and the filtered visitor hash unit to obtain contact pairs.

9. The FDEM parallel contact detection method based on dual-domain activation according to claim 8, characterized in that, The specific method for obtaining the contact pair in step S404 is as follows: Contact detection is performed on all TET4 units within the host hash unit except itself, and contact detection is performed on each TET4 unit contained in all visitor hash units to form contact pairs.

10. The FDEM parallel contact detection method based on dual-domain activation according to claim 1, characterized in that, Step S5 is as follows: S501, First round of scanning phase: Parallel traversal of each active hash unit and its corresponding visitor hash unit, counting of unit pairs that meet the contact conditions, generating the number of potential contact pairs corresponding to each active hash unit, and calculating the prefix and index information of each active hash unit based on the number of contact pairs, where the prefix and index information is used to determine the write start position of each active hash unit in the global contact pairs. S502, Second round of scanning phase: Based on the prefix and index information of the hash unit of each active state, a non-overlapping global contact pair set writing interval is allocated for each active state hash unit. Each computing thread writes all the identified potential contact pairs into the global contact pair set in parallel within its corresponding writing interval. S503. The global contact pair set is used as input data to drive subsequent contact force calculations or collision response solutions to update the motion state of the TET4 unit.

Images