Energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters

Abstract

This invention discloses an energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters. First, the 3D deformable body model is discretized into a particle system, and then divided into multiple overlapping neighboring particle clusters through voxelization. For each cluster, the local deformation gradient is fitted using the initial and current positions of the particles using a weighted least squares method. The volume and shear constraint functions of the Neo-Hookean model are combined into coupling constraints, and the gradient is calculated. At the current time step, a global energy function is constructed by combining the system momentum term and the potential energy term of all cluster coupling constraints. An energy descent block coordinate optimization strategy is adopted, iteratively solving the local problem particle by particle until the termination condition is met. Simultaneously, an index mapping between vertices and clusters is established. A graph is constructed and colored according to cluster dependencies, and vertices of the same color are assigned to parallel computing units for synchronous solution and update. This invention can alleviate the problems of slow convergence and poor stability of traditional methods under strong nonlinearity, achieving faster energy descent and higher simulation quality under the same computing resources.

Description

Technical Field

[0001] This invention relates to the field of computer graphics technology, and in particular to a simulation method for optimizing the coordinates of energy-decreasing blocks based on coupling constraints and overlapping clusters. Background Technology

[0002] With the increasing demands for physical realism in fields such as virtual reality, computer games, and surgical simulations, the need for real-time software simulation is growing. While high-precision physical models (such as the finite element method) can accurately simulate the nonlinear behavior of materials, their high computational cost makes it difficult to meet the frame rate requirements of real-time interaction. Simplified models designed to improve performance (such as traditional volume-based dynamics (PBD) or position-based extended dynamics (XPBD)) often encounter problems like convergence difficulties, volume collapse, surface wrinkling, and energy oscillations under strongly nonlinear and large deformation conditions, sacrificing simulation stability and physical realism. One improvement approach in existing technologies is to use meshless cluster structures to organize deformable bodies, decomposing the global system into local subproblems to improve computational efficiency. However, when dealing with nonlinear material models such as Neo-Hookean models, these methods typically involve decoupled, sequential projection updates of multiple constraints constituting the material's constitutive relations (such as volume constraints and shear constraints). This decoupling process ignores the inherent coupling relationships between different physical constraints, which can easily lead to gradient direction conflicts and energy disturbances during iterative solutions, resulting in slow convergence and insufficient stability under limited iterative budgets. Therefore, designing a simulation method in a meshless framework that can maintain numerical stability and physical consistency while significantly improving the efficiency of solving nonlinear problems has become a pressing technical problem in this field. Summary of the Invention

[0003] The purpose of this invention is to provide a simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters. This invention effectively alleviates the problems of slow convergence and poor stability of traditional methods under strong nonlinear conditions, achieving faster energy descent speed and higher simulation quality under the same computational budget.

[0004] The technical solution of this invention: a simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, comprising the following steps:

[0005] Step 1: Discretize the three-dimensional deformable body model to be simulated into a particle system. Based on the voxelization method, divide the particle system into a set of overlapping clusters, each cluster containing multiple neighboring particles.

[0006] Step 2: For each cluster, based on the position of the particles within the cluster in the initial reference configuration and the current configuration, fit the local deformation gradient using the weighted least squares method; based on the local deformation gradient, define the combination of the volume constraint function and the shear constraint function in the Neo-Hookean model as a coupling constraint, and calculate its gradient to obtain the local deformation potential energy term of the cluster containing each particle within the cluster.

[0007] Step 3: At the current time step, for each particle in the particle system, sum the local variable kinetic energy terms of all clusters to which the particle belongs, and then add the current kinetic energy term of the particle to obtain the global energy function corresponding to the particle; adopt the block coordinate optimization strategy based on energy descent to iteratively perform local solution for each particle vertex in the system until the preset iteration termination condition is met. The object of the local solution is the local system composed of a single particle.

[0008] Step 4: Establish an index mapping between vertices and their respective clusters, construct a graph based on the cluster dependencies between vertices and perform color grouping, and assign vertices of the same color group to parallel computing units to synchronously execute the local solution and update in step S3.

[0009] In the above-mentioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, step 1 involves dividing the particle system into a set of overlapping clusters using a voxelization method, including:

[0010] A three-dimensional voxel mesh is constructed for the model. With each voxel as the center, all particles within the neighboring voxels that are less than or equal to the preset radius from the center voxel are collected to form a cluster.

[0011] The aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters also includes:

[0012] Calculate the volume of each cluster. For voxels that are completely inside the model, the volume is the cube of the voxel's side length. For surface voxels that intersect the model surface, the volume is approximated by subdividing the surface voxel into smaller sub-voxels, counting the number of sub-voxels inside the model, and summing them up.

[0013] In the aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, step 2, the local deformation gradient is calculated using the following formula:

[0014] ;

[0015] In the formula, Indicates the number of particles in the cluster. For particles in the current configuration The offset vector relative to the cluster centroid Particles in reference configuration The offset vector relative to the cluster centroid For particles The quality; It is the transpose symbol; This represents the edge vector matrix of the current configuration; This represents the inverse of the edge vector matrix of the reference configuration.

[0016] In the aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, step S2 defines the coupling constraints and their gradients as follows:

[0017] ;

[0018] ;

[0019] In the formula, For coupling constraints, The gradient represents the coupling constraint. Due to volume constraints, It is a shear constraint; The gradient is the volume constraint. It is the gradient of the shear constraint.

[0020] In the aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, in step S3, the global energy function corresponding to the particle is defined as:

[0021] ;

[0022] In the formula, It is the global target energy corresponding to the particle; It is the particle's position vector; The particle is a predicted position vector; It is the transpose symbol. It is a quality matrix; It is a constraint function vector; It is a diagonal matrix composed of volume constraint compliance parameters and shear constraint compliance parameters.

[0023] In the aforementioned simulation method for energy-descent block coordinate optimization based on coupling constraints and overlapping clusters, step S3, which involves using an energy-descent-based block coordinate optimization strategy to locally solve for each particle vertex, includes:

[0024] For a local system consisting of a single particle, fix the positions of all other particles in the system and construct the local energy function corresponding to the particle; calculate the first-order gradient and Hessian matrix of the local energy function at the current particle position; solve the linear system to obtain the particle's vertex position increment and update the particle's vertex position.

[0025] In the aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, the first-order gradient and Hessian matrix are calculated using the following formulas:

[0026] ;

[0027] In the formula, Global energy For particle position The first gradient; Global energy For particle position The Hessian matrix; For time step; For coupling constraints; The gradient represents the coupling constraint. A diagonal matrix specifically designed for Neo-Hookean materials.

[0028] In the aforementioned simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, step S4, specifically the construction of a graph based on the cluster dependency relationship between vertices and the coloring and grouping, involves:

[0029] If two vertices belong to at least one of the same clusters, then the two vertices are determined to have a dependency relationship, and a connecting edge is established between the corresponding nodes in the graph; the graph is then colored so that vertices with a dependency relationship are assigned different colors, and vertices without a dependency relationship are assigned the same color.

[0030] Compared with the prior art, the present invention has the following beneficial effects:

[0031] 1. This invention solves the Neo-Hookean model by coupling the volume and shear constraints into a unified whole, fundamentally avoiding gradient direction conflicts and energy oscillations caused by traditional decoupled projection methods. Experiments show that, with the same number of iterations, the method of this invention achieves faster energy reduction, effectively alleviates volume collapse and surface wrinkling phenomena, and exhibits superior stability under strongly nonlinear and large deformation conditions.

[0032] 2. Under the unified momentum-constrained potential energy framework, this invention uses block coordinate optimization based on energy decrease to solve the problem. Essentially, this method finds the direction of the fastest energy decrease in each local solution step, which is more rigorous in physical and energy terms than the traditional explicit constraint projection.

[0033] 3. The meshless overlapping cluster structure adopted in this invention reduces the coupling complexity of the global system. Combined with a specially designed graph coloring parallel strategy based on vertex (particle) dependencies, it can be efficiently mapped to parallel hardware such as GPUs, significantly improving the simulation frame rate while ensuring accuracy, and meeting the needs of real-time interactive applications.

[0034] 4. The coupling constraint framework proposed in this invention is not limited to the Neo-Hookean model, but can be extended to other material models or combinations of geometric constraints composed of multiple associated constraints, and has good versatility. Attached Figure Description

[0035] Figure 1 It is a cluster partitioning diagram of a meshless model;

[0036] Figure 2 This is a schematic diagram of surface voxels and internal voxels;

[0037] Figure 3 This is a diagram showing the coloring of cluster vertices;

[0038] Figure 4 This is a diagram showing the model perturbation recovery results;

[0039] Figure 5 It is a diagram of total elastic energy;

[0040] Figure 6 It is an energy difference graph;

[0041] Figure 7 It is a relative lifting graph;

[0042] Figure 8 The simulation results of the liver model using the coupling formula (a) and the decoupling formula (b) are shown in comparison under the same parameters.

[0043] Figure 9 The simulation results of the armadillo model using the coupling formula (a) and the decoupling formula (b) are shown in comparison under the same parameters.

[0044] Figure 10 The global volume change curve of the model is shown throughout the simulation process. Detailed Implementation

[0045] The present invention will be further described below with reference to the embodiments and accompanying drawings, but this should not be construed as limiting the present invention.

[0046] Example: An energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters is applicable to real-time software simulation of 3D deformable bodies in the field of computer graphics. It can effectively solve the problems of convergence difficulty, volume collapse, and energy oscillation of traditional methods under strong nonlinearity and large deformation conditions, and achieve high-precision and high-frame-rate simulation results. The specific steps are as follows:

[0047] Step 1: Discretize the 3D deformable body model to be simulated into a particle system, using particles as the basic computational unit to represent the geometric shape and physical properties of the deformable body; then, based on the voxelization method, perform overlapping cluster partitioning on the particle system, specifically:

[0048] 1.1 Particle System Discretization: The three-dimensional deformable model to be simulated is discretized into a discrete system composed of a large number of particles. The particles are used as the basic computational units, and the position of each particle under the initial reference configuration is recorded. ,quality Physical properties such as particle spatial distribution and density are set according to simulation accuracy requirements.

[0049] 1.2 3D Voxel Mesh Construction: Construct a 3D voxel mesh covering the entire deformable model, setting the voxel side length to... A uniform voxel space is generated based on the bounding box of the model, and the particle system is mapped onto the voxel grid, recording the voxel number to which each particle belongs.

[0050] 1.3 Overlapping Cluster Partitioning: Using each voxel as the center, define a cluster neighborhood radius. Collect particles from all neighboring voxels whose Euclidean distance from the central voxel is ≤ the cluster neighborhood radius, forming an independent cluster, such as... Figure 1 As shown, each cluster contains multiple neighboring particles, and there is particle overlap between different clusters (i.e., a particle can belong to multiple clusters), thereby maintaining the overall continuity of the deformable body and avoiding model breakage caused by local cluster division.

[0051] 1.4 Approximate Calculation of Cluster Volume: Cluster volume is crucial for subsequent compliance parameter calculations. Based on the positional relationship between voxels and the model surface, voxels are divided into internal voxels and surface voxels, and their volumes are calculated separately. Figure 2 As shown:

[0052] For a voxel that is entirely located inside the model, its volume is directly taken as the cube of the voxel's side length:

[0053] ;

[0054] For surface voxels intersecting the model surface, they are subdivided into smaller sub-voxels. The center of each sub-voxel is checked to see if it is located inside the model. The number of sub-voxels located inside the model is counted, and their volumes are summed to approximate the effective volume of the surface voxels.

[0055] ;

[0056] In the formula, The ratio for further division, This indicates the number of internal daughter cells.

[0057] The final volume of a cluster is the sum of the effective volumes of all the voxels it contains:

[0058] ;

[0059] in, The set of voxels contained in this cluster. The effective volume for each voxel.

[0060] Step 2: For each cluster, based on the positions of particles within the cluster in the initial reference configuration and the current configuration, fit the local deformation gradient using the weighted least squares method; based on the local deformation gradient, define the combination of the volume constraint function and the shear constraint function in the Neo-Hookean model as a coupling constraint, and calculate its gradient to obtain the local deformation potential energy term of the cluster containing all particles within the cluster; the specific operations are as follows:

[0061] 2.1 Cluster Centroid Calculation: Calculate the centroid of each particle within the cluster under both the initial reference configuration and the current configuration. The centroid is the mass-weighted average position of the particles. This calculation is performed during the preprocessing stage to reduce runtime overhead. The formula is:

[0062] ;

[0063] In the formula, Indicates the number of particles in the cluster. For particles quality Particles in reference configuration The location.

[0064] 2.2 Calculation of Particle Relative Offset Vector: Calculate the offset vector of a particle relative to the cluster centroid, reflecting the particle's local position within the cluster. The formula is:

[0065] Current configuration offset: ;

[0066] Reference configuration offset: .

[0067] 2.3 Local Deformation Gradient Fitting: The local deformation gradient (deformation matrix) of the cluster is fitted using the weighted least squares method to minimize the deformation error of all particles within the cluster under this transformation. Specifically:

[0068] ;

[0069] In the formula, Indicates the number of particles in the cluster. For particles in the current configuration The offset vector relative to the cluster centroid Particles in reference configuration The offset vector relative to the cluster centroid For particles The quality; It is the transpose symbol; This represents the edge vector matrix of the current configuration; This represents the inverse of the edge vector matrix of the reference configuration.

[0070] 2.4 Constrained Representation of Neo-Hookean Model: The stable-Neo-Hookean model is used to describe the elastic behavior of deformable bodies. Its strain energy density is decomposed into volume terms and shear terms, which are then transformed into volume constraint functions and shear constraint functions to reflect the volume change and shear deformation characteristics of particles within the cluster.

[0071] 2.5 Coupling Constraints and Gradient Definitions: Based on the fitted local deformation gradient, the volume constraint function (hydrostatic constraint) and shear constraint function (partial energy constraint) in the Neo-Hookean model are combined and defined as coupling constraints. Their gradients are calculated using the chain rule and matrix differential relations. The specific formula is as follows:

[0072] ;

[0073] ;

[0074] In the formula, For coupling constraints, The gradient represents the coupling constraint. Due to volume constraints, It is a shear constraint; The gradient is the volume constraint. It is the gradient of the shear constraint.

[0075] The values ​​of volume constraints and shear constraints are calculated by combining the local deformation gradient with the Lamé parameters and stability terms of the Neo-Hookean model, reflecting the volume changes and shear deformation characteristics of particles within the cluster.

[0076] Step 3: At the current time step, for each particle in the particle system, sum the local variable kinetic energy terms of all clusters to which the particle belongs, and then add the current kinetic energy term of the particle to obtain the global energy function corresponding to the particle; using a block coordinate optimization strategy based on energy descent, iteratively perform local solutions for each particle vertex in the system (e.g., when the number of iterations reaches a threshold, the energy change is less than a set value, etc.) until the preset iteration termination condition is met. The object of the local solution is the local system composed of a single particle; the specific operation is as follows:

[0077] 3.1 Establishment of Global Energy Function: The global target energy consists of the system momentum term and the coupled constraint potential energy term, defined by the following formula:

[0078] ;

[0079] In the formula, It is the global target energy; It is the particle's position vector; It is the predicted position vector; It is the transpose symbol. It is a quality matrix; It is a constraint function vector; It is a diagonal matrix composed of volume constraint compliance parameters and shear constraint compliance parameters.

[0080] 3.2 Local solution for a single particle vertex: With the positions of all other particles in the system fixed, the optimization solution is performed only on the position of the vertex of this particle. The specific process is as follows:

[0081] 3.2.1 Construct the local energy function corresponding to the vertex of the particle. This local energy function is a simplified form of the global energy function when the positions of other particles are fixed.

[0082] 3.2.2 Calculate the first-order gradient and Hessian matrix of the local energy function at the current particle position. The calculation formula is as follows:

[0083] ;

[0084] In the formula, Global energy For particle position The first gradient; Global energy For particle position The Hessian matrix; For time step; For coupling constraints; The gradient represents the coupling constraint. Diagonal matrix specifically for Neo-Hookean materials:

[0085] .

[0086] 3.2.3 Solving linear systems This gives the position increment of the vertex and updates the vertex position.

[0087] 3.2.4 Repeat steps 3.2.1-3.2.3 above for the updated vertex position until the preset iteration termination condition is met, thus completing the local solution for the particle vertex.

[0088] Step 4: Establish an index mapping between vertices and their respective clusters. Construct a graph based on the cluster dependencies between vertices and perform color grouping. Assign vertices in the same color group to parallel computing units to synchronously execute the local solution and update in step S3. This effectively avoids write conflicts in parallel computing and improves simulation efficiency. The specific operations are as follows:

[0089] 4.1 Establishing a vertex-cluster index mapping: An index table is established for each particle vertex, recording the numbers of all clusters to which the vertex belongs, to achieve fast lookup and association between vertices and clusters.

[0090] 4.2 Constructing a vertex dependency graph: Using each particle vertex as a node in the graph, if two vertices belong to at least one of the same clusters, then it is determined that these two vertices have a dependency relationship. An undirected connection edge is established between the corresponding two nodes in the graph to complete the construction of the vertex dependency graph.

[0091] 4.3 Graph Coloring Grouping: A graph coloring algorithm is used to color the constructed vertex dependency graph. The coloring rule is: vertices with dependencies are assigned different colors, and vertices without dependencies can be assigned the same color; ultimately, each color corresponds to a group of vertices, and there are no cluster dependencies between vertices within the same group, such as... Figure 3 As shown.

[0092] 4.4 Parallel Solving and Updating: The vertex set of each color group is assigned to an independent parallel computing unit (such as different sub-threads of the GPU). Each parallel computing unit synchronously performs the local solving and position update operations in step 3 on the vertices within the group. After the parallel computing of all color groups is completed, all particle vertices of the current time step are updated, and the simulation calculation of the next time step is started.

[0093] The method in this embodiment avoids gradient conflicts and energy oscillations in traditional decoupled projection by solving the volume and shear constraints of the Neo-Hookean model in a coupled manner; the combination of block coordinate optimization strategy with energy descent ensures the rigor of physical and energy meanings; and the combination of overlapping cluster structure and graph coloring parallel strategy greatly improves computational efficiency while ensuring simulation accuracy.

[0094] Furthermore, the method described in this embodiment was verified through experiments. All CPU experiments were run on a single thread on a 2.30GHz Core i5-8300H processor with 16GB of memory. GPU experiments were run on an NVIDIA GeForce GTX 1050 Ti graphics card, and the simulation results were rendered using Polyscope. Table 1 shows the data such as the number of model vertices, clusters, number of vertex shading colors, and model memory size (MB) used in the experiments.

[0095] Table 1

[0096]

[0097] To evaluate the actual acceleration performance of the parallel graph coloring strategy, 200 simulated frames were used with the same material parameters. Under the conditions of 4 sub-steps per frame and 8 iterations per sub-step, the average time consumption (milliseconds / frame) per frame for the baseline solver and the coupled constraint method on the GPU was measured, and various models were compared. Table 2 summarizes the time consumption per frame for different models under the two partitioning methods, quantifying the actual performance improvement brought by cluster partitioning and graph coloring.

[0098] Table 2

[0099]

[0100] As can be seen from Table 2, under the same experimental conditions, the average computation time (milliseconds / frame) of the coupled constraint solver (i.e. the method described in this invention) for all test models is much lower than that of the baseline solver, achieving significant computational acceleration. Moreover, the more vertices / clusters the model has (the higher the complexity), the more significant the acceleration effect.

[0101] To further verify the robustness and stability of the proposed solution system under extreme conditions, a large-scale vertex perturbation experiment was designed and conducted. In this experiment, all vertices of the model were randomly perturbed from their initial equilibrium positions to within ±80%. Such a large random offset would cause the initial geometry to completely lose its original structural characteristics, resulting in a highly chaotic and almost unrecognizable mesh configuration. This setup aims to test whether the solution system can still maintain stable iteration and gradually recover to a reasonable geometry under strongly nonlinear and large deformation conditions.

[0102] After vertex perturbation, the model is dynamically reconstructed under uniform simulation parameter settings to eliminate the influence of factors such as time step, substep size, or number of iterations on stability. The Neo-Hookean material model is used as the primary type of elastic constraint in the simulation. This model exhibits good volume preservation characteristics under large deformation conditions, providing a physically consistent recovery trend for the system. Furthermore, to ensure simulation stability and numerical accuracy, the time step is fixed at s, with four substeps executed per frame, and 16 iterations of constraint solving within each substep. This solution strategy balances computational efficiency and convergence characteristics, effectively avoiding unstable behaviors such as divergence, jitter, or mesh flipping under large perturbation conditions.

[0103] Under the above conditions, driven by Neo-Hookean energy, the model gradually converges from a completely random chaotic state to the original smooth geometry. Figure 4Several representative simulation frames from the experiment are presented. It can be seen that the recovery process, from the initial disordered deformation to the final structured geometry, is smooth and without significant oscillations, demonstrating strong stability and tolerance to extreme inputs. This result intuitively shows that the proposed method can maintain stable solutions even when faced with large-scale vertex perturbations, exhibiting good robustness and numerical stability.

[0104] To evaluate the convergence of the solver, an experimental configuration exhibiting the highest total elastic potential energy during the simulation was selected in the beam example, and a solution process of 100 sub-steps was performed on it. Subsequently, based on the cluster-partitioned model structure, a constrained objective function was used and Newton's method was employed for solution. To verify the impact of the constraint coupling strategy on convergence, the solver was run under conditions of no constraint coupling and full constraint coupling, and the elastic energy decreasing gradually with the number of iterations was compared to evaluate the convergence improvement brought by the coupling strategy. It should be noted that the beam model in this experiment consists of 3460 clusters, and the material Lamé coefficient is... , .

[0105] The result is as follows Figures 5-7 As shown. Figure 5 The diagram shows the total elastic energy. Figure 6 This is an energy difference plot; Figure 7 This is a relative lifting graph. During the computation of a single sub-step, from... Figure 5 It is evident that the potential energy reduction of the model obtained using the coupled constraint processing method is greater than that obtained using the decoupling method. This difference is not only reflected in the numerical convergence properties but also directly in the model's visualization. The difference is particularly significant when the baseline solver and the coupled solver are compared under the same computational budget, material parameters, and simulation parameters. Specifically, under a limited iteration budget, decoupling and updating each constraint often makes it difficult to coordinate the interrelated physical constraints in a timely manner, which may introduce high-frequency numerical errors with local inconsistencies, manifesting as pseudo-deformations such as surface wrinkles or jitter. In contrast, the model using coupled updates can synchronously correct related variables in a single iteration, effectively suppressing oscillations and thus exhibiting more continuous and stable surface deformation. To further illustrate the improvement effect of this method, energy difference curves are given respectively ( Figure 6 ) and relative lifting curve ( Figure 7 ). Figure 6 This indicates that the coupling method reduces energy more significantly in a single sub-step than the traditional method; Figure 7 Further evidence confirms that, with the same number of iterations, the coupling method significantly outperforms both XPBD and the basic method.

[0106] Figure 8The simulation results of the liver model using the coupling formula (a) and the decoupling formula (b) are shown in comparison under the same parameters. Figure 9 The simulation results of the armadillo model using the coupling formula (a) and the decoupling formula (b) are compared under the same parameters. It can be seen that, under the same parameter settings, the model using the coupling formula exhibits better stability and convergence behavior across multiple time steps, while the model using the decoupling formula struggles to maintain stability under the same computational conditions.

[0107] In the preceding experiments, it was verified that the constraint coupling strategy can effectively improve the convergence of the cluster-partition-based solver during the elastic energy descent process. The following section further examines its applicability to other types of constraint combinations. Since volume conservation constraints possess strong globality and high nonlinearity, they are very suitable as an effective use case for testing the generalization ability of constraint coupling methods. Therefore, this experiment evaluates the effectiveness of the constraint coupling strategy in preserving volume constraints in cluster structure models by comparing the volume stability during simulation with the baseline solver, the strategy coupled with Neo-Hookean constraints, and the strategy coupled with volume constraints and edge distance constraints, thereby verifying the generalizability of this strategy. Figure 10 The global volume change curve of the model throughout the simulation is presented. It can be seen that the coupling of volume constraints and distance constraints also converges stably in this example, eventually approaching a volume level similar to that of the Neo-Hookean coupling. This experiment demonstrates that the proposed coupling framework does not depend on a specific physical energy form and has the potential to be extended to combinations of geometric constraints and even non-physical constraints.

[0108] In summary, this invention proposes an energy descent block coordinate optimization method based on coupling constraints and overlapping clusters. This method integrates the constraint modeling advantages of XPBD with the energy iteration framework of VBD, significantly improving solution stability and convergence performance while maintaining energy consistency. By dividing the object into overlapping cluster structures and using weighted least squares fitting to approximate the deformation of the continuum, this method can construct local deformation gradients without relying on mesh topology. Furthermore, the graph coloring parallel strategy proposed in this invention effectively avoids write conflict problems, allowing the vertex-based local energy accumulation process to be efficiently mapped to the GPU parallel architecture, further improving computational efficiency while ensuring accuracy.

[0109] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters, characterized in that, Includes the following steps: Step 1: Discretize the three-dimensional deformable body model to be simulated into a particle system. Based on the voxelization method, divide the particle system into a set of overlapping clusters, each cluster containing multiple neighboring particles. Step 2: For each cluster, based on the position of the particles within the cluster in the initial reference configuration and the current configuration, fit the local deformation gradient using the weighted least squares method; based on the local deformation gradient, define the combination of the volume constraint function and the shear constraint function in the Neo-Hookean model as a coupling constraint, and calculate its gradient to obtain the local deformation potential energy term of the cluster containing each particle within the cluster. Step 3: At the current time step, for each particle in the particle system, sum the local variable kinetic energy terms of all clusters to which the particle belongs, and then add the current kinetic energy term of the particle to obtain the global energy function corresponding to the particle; adopt the block coordinate optimization strategy based on energy descent to iteratively perform local solution for each particle vertex in the system until the preset iteration termination condition is met. The object of the local solution is the local system composed of a single particle. Step 4: Establish an index mapping between vertices and their respective clusters, construct a graph based on the cluster dependencies between vertices and perform color grouping, and assign vertices of the same color group to parallel computing units to synchronously execute the local solution and update in step S3.

2. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 1, characterized in that, In step 1, the voxelization-based method divides the particle system into a set of overlapping clusters, including: A three-dimensional voxel mesh is constructed for the model. With each voxel as the center, all particles within the neighboring voxels that are less than or equal to the preset radius from the center voxel are collected to form a cluster.

3. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 2, characterized in that, Also includes: Calculate the volume of each cluster, where the volume of a voxel that is completely inside the model is the cube of the voxel's side length. For a surface voxel that intersects with the model surface, its volume is approximated by subdividing the surface voxel into smaller sub-voxels, counting the number of sub-voxels located within the model, and summing them up.

4. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 1, characterized in that, In step 2, the local deformation gradient is calculated using the following formula: ； In the formula, Indicates the number of particles in the cluster. For particles in the current configuration The offset vector relative to the cluster centroid Particles in reference configuration The offset vector relative to the cluster centroid For particles The quality; It is the transpose symbol; This represents the edge vector matrix of the current configuration; This represents the inverse of the edge vector matrix of the reference configuration.

5. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 1, characterized in that, In step S2, the coupling constraint and its gradient are defined as follows: ； ； In the formula, For coupling constraints, The gradient represents the coupling constraint. Due to volume constraints, It is a shear constraint; The gradient is the volume constraint. It is the gradient of the shear constraint.

6. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 5, characterized in that, In step S3, the global energy function corresponding to the particle is defined as: ； In the formula, It is the global target energy corresponding to the particle; It is the particle's position vector; The particle is the predicted position vector; It is the transpose symbol. It is a quality matrix; It is a constraint function vector; It is a diagonal matrix composed of volume constraint compliance parameters and shear constraint compliance parameters.

7. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 6, characterized in that, In step S3, the step of using an energy-decreasing block coordinate optimization strategy to locally solve for each particle vertex includes: For a local system consisting of a single particle, fix the positions of all other particles in the system and construct the local energy function corresponding to the particle; calculate the first-order gradient and Hessian matrix of the local energy function at the current particle position; solve the linear system to obtain the particle's vertex position increment and update the particle's vertex position.

8. The energy descent block coordinate optimization simulation method based on coupling constraints and overlapping clusters according to claim 7, characterized in that, The first-order gradient and the Hessian matrix are calculated using the following formula: ； In the formula, Global energy For particle position The first gradient; Global energy For particle position The Hessian matrix; For time step; For coupling constraints; The gradient represents the coupling constraint. A diagonal matrix specifically designed for Neo-Hookean materials.

9. The simulation method for energy descent block coordinate optimization based on coupling constraints and overlapping clusters according to claim 1, characterized in that, In step S4, the specific steps of constructing a graph based on the cluster dependency relationship between vertices and performing coloring and grouping are as follows: If two vertices belong to at least one of the same clusters, then the two vertices are determined to have a dependency relationship, and a connecting edge is established between the corresponding nodes in the graph; the graph is then colored so that vertices with a dependency relationship are assigned different colors, and vertices without a dependency relationship are assigned the same color.

Images