An Improved Position Dynamics Simulation Method Based on Algebraic Multigrid

Abstract

An improved position dynamics simulation method based on algebraic multigrids is presented, belonging to the field of computer graphics physical simulation technology. In the iterative solution of each simulation substep, this method first uses a semi-implicit Euler method to calculate the predicted vertex position, initializing the Lagrange multipliers corresponding to all constraints to zero vectors; it then successively calculates the constraint values ​​and gradients of all constraints, assembling a global system matrix and right-hand vector that preserves the off-diagonal terms corresponding to the adjacency relationships between constraints; it then uses an aggregation-based algebraic multigrid method to solve the linear equation system, obtaining the Lagrange multiplier increments, and synchronously updating the vertex positions and Lagrange multipliers; finally, it uses the dual residual value as the convergence criterion, and after reaching a preset accuracy, calculates the vertex velocity based on the difference between the old and new positions, completing the current substep calculation. This invention designs a matrix-free sparse assembly scheme for spring-mass cloth and tetrahedral soft systems, making full use of the sparsity of the system matrix to reduce memory overhead. It achieves efficient solution through multi-scale algebraic multigrids, improving the convergence speed by one to two orders of magnitude compared to the traditional XPBD method. It effectively solves the typical defects of excessive cloth stretching and soft body stiffness in high-resolution simulation, balancing the accuracy of physical simulation and the running efficiency. It can be widely used in film and animation, virtual reality, industrial mechanics simulation and other fields.

Description

Technical Field

[0001] This invention belongs to the field of physical simulation technology in computer graphics, and specifically relates to an improved position dynamics simulation method based on algebraic multigrids. Background Technology

[0002] In the field of computer graphics, the physical simulation of deformable objects such as cloth and soft bodies is widely used in film and animation, virtual reality, and interactive games. The physical realism and operational efficiency of the simulation directly determine the immersion and interactive experience of the content. Position-Based Dynamics (PBD), which uses vertex positions as the direct solution variable, has the advantages of low computational cost, strong robustness, and ease of implementation. It is currently the mainstream technical solution for real-time physical simulation and can meet the real-time requirements of interactive applications in low-resolution scenarios.

[0003] However, traditional PBD methods have inherent convergence limitations: each simulation substep typically performs only one Gaussian-Seidel or Jacobi iteration, updating vertex positions sequentially for each constraint, without fully considering the coupling effects between constraints, thus failing to iteratively converge to an accurate physical equilibrium state. This deficiency is particularly pronounced in high-resolution simulation scenarios: even with a significant increase in the number of iterations, it remains difficult to eliminate visual imperfections such as excessive fabric stretching and insufficient softness, resulting in insufficient physical accuracy of the simulation results and making it difficult to support the demands of high-precision film and television-level simulations and industrial simulations.

[0004] To improve the convergence performance of Position-Based Dynamics (PBD), Extended Position-Based Dynamics (XPBD) was proposed. This method introduces Lagrange multipliers and compliance parameters, transforming the constraint-solving problem into a system of linear equations about the increments of Lagrange multipliers. Compliance is used instead of stiffness as the physical parameter, improving parameter stability and convergence. However, existing XPBD schemes still haven't broken through the single-step iterative solution mode and haven't designed efficient and accurate solvers for the derived linear systems. When the simulation model resolution increases and the number of constraints grows significantly, the scale of the linear system increases dramatically. Using conventional iterative methods results in slow convergence and numerous iterations, failing to fully utilize the accuracy potential of XPBD and struggling to balance simulation efficiency and physical realism. The problems of overstretching and excessive softness at high resolutions remain unresolved.

[0005] Furthermore, existing linear solution schemes for position dynamics lack optimization for constrained, sparse structures: the system matrices of different simulation objects, such as spring-mass cloth and tetrahedral soft bodies, exhibit typical sparsity and local adjacency characteristics. General solvers cannot fully utilize these structural characteristics to achieve acceleration, further restricting the running efficiency of high-resolution simulations.

[0006] Therefore, there is an urgent need for an improved position dynamics simulation method that can efficiently solve large-scale constrained linear systems, significantly improve convergence speed, and take into account both simulation accuracy and operational efficiency, in order to solve the technical problems of excessive cloth stretching and soft body stiffness distortion in high-resolution scenes. Summary of the Invention

[0007] The technical problem to be solved by this invention is to provide an improved position dynamics simulation method based on algebraic multigrids. This invention designs a matrix-free sparse assembly scheme for spring-mass cloth and tetrahedral soft systems, making full use of the sparsity of the system matrix to reduce memory overhead. It achieves efficient solution through multi-scale algebraic multigrids, improving the convergence speed by one to two orders of magnitude compared with the traditional XPBD method. It effectively solves the typical defects of excessive cloth stretching and soft body stiffness in high-resolution simulation, and balances the accuracy of physical simulation with running efficiency.

[0008] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows: An improved position dynamics simulation method based on algebraic multigrids, comprising the following steps: The solution process for each simulation substep includes the following steps: S1. Use the semi-implicit Euler method to calculate the predicted vertex position and record the old vertex position at the previous moment; S2. Initialize all Lagrange multiplier vectors corresponding to constraints to zero vectors; S3. Determine whether the current iteration result has reached the preset convergence accuracy. If it has, jump to step S8; otherwise, continue to execute the subsequent steps. S4. Calculate the constraint values ​​and constraint gradients of all constraints, and generate constraint vectors and constraint gradient matrices; S5. Assemble the global system matrix and right-hand vector based on the constraint gradient matrix, the mass diagonal matrix and the compliance parameter. The system matrix retains the off-diagonal terms corresponding to the adjacency relationship between constraints. S6. Solve the linear equations corresponding to the system matrix using the algebraic multigrid method to obtain the Lagrange multiplier increments; S7. Calculate the vertex position correction based on the Lagrange multiplier increment, update the Lagrange multipliers and vertex position synchronously, and return to step S3 to enter the next iteration. S8. Calculate the vertex velocity based on the difference between the updated vertex position and the old position, and complete the calculation of the current simulation sub-step.

[0009] Preferably, in step S5, the system matrix and the right-hand vector are derived through constraint gradients, quality parameters, and compliance parameters, and the calculation formulas are as follows: ; ; In the formula: The global system matrix has dimensions of . ; The constraint gradient matrix has dimensions of . ; It is a mass diagonal matrix with dimension . ; For mass diagonal matrix The inverse matrix; Constraint gradient matrix The transpose of the matrix; It is a flexibility diagonal matrix with dimension . ; Let be the right-hand vector of the system, with dimension . ; Let be a constraint value vector with dimension . ; Let be a Lagrange multiplier vector with dimension . ; This represents the total number of constraints in the simulation scenario. This represents the total number of vertices in the simulation scene.

[0010] Preferably, step S5 employs a matrix-free method to quickly assemble the system matrix elements. Leveraging the sparse structure of the system matrix, element calculations can be completed without explicitly storing the complete matrix. This matrix-free method is executed according to the discrete structure type of the simulation object, as follows: For a spring-mass point cloth system, each row of the system matrix corresponds to a spring constraint; the diagonal elements of the row are the sum of the reciprocals of the masses of the two vertices of the spring and the corresponding flexibility parameter; the off-diagonal elements of the row are calculated from the gradient of the adjacent springs that share vertices with the spring and the reciprocal of the mass of the shared vertices, and the number of off-diagonal terms is the same as the number of adjacent springs of the spring. For a tetrahedral soft system, each row of the system matrix corresponds to a tetrahedral constraint; the diagonal elements of the row are the sum of the products of the inverse mass of the four vertices of the tetrahedron and the square of the corresponding gradient; the off-diagonal elements of the row are calculated from the gradients of the adjacent tetrahedrons that share vertices with the tetrahedron and the inverse mass of the shared vertices, and the number of off-diagonal terms is the same as the number of adjacent tetrahedrons of the tetrahedron.

[0011] Preferably, step S6 employs an aggregation-based algebraic multigrid method to solve the linear equation system. The setup phase of the algebraic multigrid method includes the following sub-steps: S6.1 Perform connection strength filtering on the system matrix, remove diagonal terms and weakly connected off-diagonal terms, and generate a strength matrix; S6.2 Determine the connection strength between nodes based on the off-diagonal values ​​of the strength matrix, aggregate nodes with strong mutual influence into non-overlapping and complete aggregate partitions, and generate an aggregate matrix of 0 / 1 form; S6.3. Insert each component of the near-kernel vector into the corresponding position of the aggregation matrix according to the aggregation partition to form a temporary extension matrix; perform QR decomposition on the temporary extension matrix according to the aggregation block, and use the decomposition result as the temporary extension matrix and near-kernel vector of the next level. S6.4 Smooth the temporary extension matrix to obtain the extension operator from the current level to the next level; S6.5 Determine if the coarsest level has been reached. If not, proceed to the next level and repeat steps S6.2 to S6.4 until the operator construction for all levels is completed.

[0012] Preferably, the solution stage of the aggregation-based algebraic multigrid method includes the following sub-steps: S6.6 Perform iterative smoothing on the top-level linear system, calculate the residual of the current layer, and pass the residual to the next layer through the constraint operator to construct the next layer linear system; repeat the smoothing, residual calculation and passing operations until the bottom layer is reached to obtain the residuals of all levels; the constraint operator is the transpose of the corresponding extension operator. S6.7. Solve the lowest level residual linear system using the direct method to obtain the lowest level error vector; S6.8. Pass the bottom-level error vector to the next level through the extension operator, correct the residuals of the corresponding level and perform a smoothing operation; pass and correct it layer by layer until it returns to the top level to obtain the corrected solution vector. S6.9 Repeat steps S6.6 to S6.8 for the V loop operation until the residual meets the preset accuracy requirements, and output the final Lagrange multiplier increment.

[0013] Preferably, in step S3, the dual residual value is used as the criterion for judging the convergence accuracy of the iteration, and the formula for calculating the dual residual value is: ; In the formula: The dual residual value is used to characterize the convergence of the iterative solution; A constraint value vector; This is the ratio of the compliance parameter to the square of the time step, where the compliance parameter and the time step are preset by the simulation scenario. is a Lagrange multiplier vector.

[0014] Preferably, in step S7, the vertex position correction is derived by using the Lagrange multiplier increment and the constraint gradient, and the calculation formula is as follows:

[0015] In the formula: The vertex position correction vector has a dimension of . ; It is the inverse matrix of the mass diagonal matrix; This is the transpose of the constraint gradient matrix; The Lagrange multiplier increment vector has dimension . ; The calculation uses a matrix-free method, which accumulates the product of the gradient of the corresponding constraint and the Lagrange multiplier increment along the vertex dimension, and then multiplies it by the reciprocal of the vertex quality to obtain the single vertex position correction.

[0016] An improved position dynamics simulation system based on algebraic multigrids includes a prediction initialization module, a constraint calculation module, a matrix assembly module, a multigrid solving module, an iterative update module, and a velocity calculation module, which are sequentially connected by communication. The prediction initialization module is used to calculate the predicted vertex positions using the semi-implicit Euler method, record the old positions, and initialize the Lagrange multipliers. The constraint calculation module is used to calculate the constraint values ​​and constraint gradients of all constraints, and generate constraint vectors and constraint gradient matrices. The matrix assembly module is used to assemble the global system matrix and the right-hand vector while retaining adjacent off-diagonal terms; The multigrid solution module is used to solve linear equation systems using the algebraic multigrid method to obtain the Lagrange multiplier increments. The iterative update module is used to calculate the vertex position correction, update the Lagrange multipliers and vertex positions, and determine the convergence of the iteration. The velocity calculation module is used to calculate the vertex velocity based on the difference between the old and new positions after the iteration converges.

[0017] Preferably, the multigrid solving module has built-in hierarchical building units and multi-level solving units; the hierarchical building units are used to generate multi-scale extension operators and constraint operators to complete the setup stage of algebraic multigrids; the multi-level solving units are used to execute the V-loop solving process of residual propagation, low-level solving and error correction to complete the solution stage calculation of algebraic multigrids.

[0018] A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the improved position dynamics simulation method based on algebraic multigrids.

[0019] The present invention can achieve the following beneficial effects: 1. This invention addresses the technical bottlenecks of poor convergence and high iteration steps in traditional PBD and XPBD methods for high-resolution scenarios. It introduces an algebraic multigrid method into the solution process of linear systems in position dynamics, achieving a rapid reduction of dual residuals through multi-scale residual propagation and error correction mechanisms. Compared to the traditional XPBD method, the number of iteration steps required to achieve the same convergence accuracy can be reduced by one to two orders of magnitude, significantly compressing the computational overhead of a single simulation substep.

[0020] 2. In the global system matrix assembly process, this invention fully preserves the off-diagonal terms corresponding to the adjacency relationships between constraints, and fully incorporates the coupling effects between different constraints. This replaces the traditional approximate solution mode that updates constraints independently, making the solution results closer to the real physical equilibrium state. It effectively solves typical visual defects commonly found in traditional methods, such as excessive fabric stretching and soft body stiffness. The simulation physical accuracy under high-resolution models is significantly improved, meeting the quality requirements of high-precision film and television production and industrial simulation.

[0021] 3. This invention designs matrix-free computation schemes for two mainstream simulation objects: spring-mass cloth and tetrahedral soft bodies. Utilizing the sparse structure of the system matrix, it eliminates the need to explicitly store the complete large-scale matrix; matrix element calculations and vector operations can be completed solely through local topological information from adjacency constraints. This significantly reduces memory usage while maintaining solution accuracy, making it more suitable for ultra-large-scale simulation scenarios with hundreds of thousands of vertices.

[0022] 4. This invention adopts an aggregation-based algebraic multigrid solution framework. Through steps such as connection strength filtering, strong connection node aggregation, and QR decomposition optimization extension operator assisted by near kernel vector, it constructs an efficient and stable multi-scale solution hierarchy. Combined with the residual propagation and error correction process of V-loop, it has a stable acceleration effect on constrained linear systems of different scales and types. Moreover, the larger the simulation model, the more prominent the acceleration advantage of the multigrid method becomes.

[0023] 5. This invention uses dual residuals as the quantitative criterion for iterative convergence, allowing for flexible setting of convergence accuracy thresholds according to application scenario requirements, balancing the operational efficiency of real-time interactive scenarios with the quality requirements of offline high-precision simulation. The method is compatible with both spring-mass point placement systems and tetrahedral soft systems, and can be widely applied in multiple technical fields such as game interaction, film and animation, and industrial mechanics simulation. Attached Figure Description

[0024] The present invention will be further described below with reference to the accompanying drawings and embodiments: Figure 1 Overall flowchart of this method; Figure 2 A schematic diagram of the sparse matrix structure in the cloth simulation of this method; Figure 3 A schematic diagram of the sparse matrix structure in the software simulation of this method; Figure 4 The dual residuals and energy convergence curves of this method in cloth simulation; Figure 5 Dual residuals and energy convergence curves of XPBD in cloth simulation; Figure 6 A comparison of the results of this method with XPBD in cloth simulation; Figure 7 The dual residuals and energy convergence curves of this method in software simulation; Figure 8 Dual residuals and energy convergence curves of XPBD in software simulation; Figure 9 The comparison diagram between this method and XPBD in software simulation. Detailed Implementation

[0025] Preferred solutions include Figures 1 to 9 As shown, an improved position dynamics simulation method based on algebraic multigrids is proposed. Figure 1 A flowchart according to this method is shown, and the invention will be further described below with reference to other accompanying drawings and specific embodiments.

[0026] The main steps of the overall algorithm are as follows: Input: Vertex velocity and position at the previous time step, and the simulation mesh.

[0027] Output: The vertex velocity and position at the next moment.

[0028] Step 1: Obtain the predicted position based on the semi-implicit Euler equation and record the old position.

[0029] Step 2: Set the Lagrange multipliers to zero.

[0030] Step 3: Determine if the accuracy meets the requirements. If it does, proceed to step 8; otherwise, proceed to step 4.

[0031] Step 4: Calculate the constraints and constraint gradients.

[0032] Step 5: Assemble the system matrix and the right-hand vector.

[0033] Step 6: Solve the system matrix using the multigrid method to obtain the Lagrange multiplier increments.

[0034] Step 7: Calculate the position correction based on the Lagrange multiplier increment, and correct the Lagrange multipliers and position.

[0035] Step 8: Calculate the speed based on the difference between the new position and the old position.

[0036] The detailed implementation methods for steps 1 to 8 are as follows: Step 1: Calculate the predicted position using the semi-implicit Euler method and record the old position.

[0037] Predicted location for:

[0038] in It is the position at the previous moment. It is the speed of the previous moment. It is the time step. It is a quality matrix. It is an external force calculated based on the position at the previous moment. The predicted position can be directly assigned, that is, the position will be used next... Iterate from the starting point.

[0039] Step 2: Set the Lagrange multipliers to zero.

[0040]

[0041] in The size is m. m is the number of constraints.

[0042] Step 3: Iterate until the maximum number of iterations is reached or exit early due to meeting the accuracy requirements.

[0043] Step 4: Calculate the constraints and constraint gradients.

[0044] This step requires calculations based on the type of constraint. The constraint calculation method is consistent with that in traditional PBD. The following example uses the distance constraint of the fabric element:

[0045] in and These are the two vertices of the distance constraint spring. It is a 2-norm. It is the initial original length of the spring, that is, the length before the simulation begins. .

[0046] The constraint gradient g of a single spring is calculated as follows:

[0047]

[0048] Since in the discrete form of the spring-mass system, a discrete element (i.e., a spring) has two vertices, the gradient has two values. and , respectively corresponding to about (Vertex 1) and about The gradient of (vertex 2). The two are exactly opposites of each other.

[0049] calculate The first equal sign indicates that the constraint gradient is C with respect to... The gradient of C is calculated and normalized. The second equal sign is the gradient calculation formula obtained under the assumption that C is a distance constraint.

[0050] Before discretization, constraint C can be viewed as a scalar function of position x, which becomes an m-dimensional vector after discretization. x, after discretization, is 3n-dimensional. Therefore, the partial derivative of C with respect to x becomes a matrix of size m x 3n. This is matrix G.

[0051] The G matrix is ​​a large sparse matrix of size m x 3n, which has the following form:

[0052] Now, using cloth simulation as an example, we will introduce the sparse structure of the system matrix. The corresponding schematic diagram is shown below. Figure 2 As shown. Suppose we only consider two connected springs in the system, numbered i and ia. ia shares vertex a with i. Spring i has two vertices numbered a and b. Spring ia has two vertices numbered a and c.

[0053] Consider matrix G in a Compact Sparse Row (CSR) manner. G has m rows, corresponding to m springs. The number of non-zero elements in each row of G is determined by the number of vertices of the discrete element. For example, in a spring-mass system, each row has two gradients. In the example, spring i corresponds to the i-th row in G.

[0054] The columns of G correspond to the vertex gradients, so the column positions are determined by the vertex numbers. In the 3D case, the gradient of a single vertex is a 3x1 vector. Therefore, G can be viewed as a block matrix, with each block being a 3x1 vector representing the gradient of a vertex. Thus, in this example, the 3a, 3a+1, and 3a+2 columns of the i-th row correspond to the gradient of spring i with respect to vertex a. It is worth noting that in G, since the two vertices connected by the same spring are not necessarily numbered consecutively, therefore... and They are not necessarily adjacent in G.

[0055] A similar analysis process exists in software simulation, so it will not be elaborated upon here.

[0056] Step 5: Assemble the system matrix and the right-hand vector.

[0057] The method for calculating the system matrix A is as follows:

[0058] Where A is the system matrix of size m x m, and G is the constraint gradient matrix of size m x 3n. It is the transpose of G. M is the mass diagonal matrix. It is the flexibility diagonal matrix divided by the square of time. n is the number of vertices, and m is the number of constraints.

[0059] Sparse structure of system matrix in cloth simulation, such as Figure 2 As shown, the diagonal element in the i-th row of the system matrix that makes up the fabric (corresponding to the spring numbered i) is: Its j-th column (corresponding to the adjacent spring with number j) has the following non-diagonal elements: .

[0060] Similarly, the sparse structure of matrix A in a software simulation system is as follows: Figure 3 As shown. The diagonal element of the i-th row (corresponding to the tetrahedron numbered i) in the system matrix that makes up the software is: Its j-th column (corresponding to the adjacent tetrahedron numbered j) has the following non-diagonal elements: Where k = 0, 1, 2, 3 represents the local numbering of the four vertices of the tetrahedron. This represents the vertex with local number k in tetrahedron i. sv(sharedvertex) represents the number of the shared vertex.

[0061] The method for calculating the right-hand vector b is as follows:

[0062] The constraint vector C, calculated in step 4, is an m x 1 vector. (Lagrange multiplier vector) It is continuously updated in the iteration, accumulating from zero.

[0063] Step 6: Solve for the system matrix to obtain the Lagrange multiplier increments.

[0064] Solve the following linear system.

[0065]

[0066] Obtain the Lagrange multiplier increment .

[0067] The solution methods can include various approaches for solving large-scale sparse linear systems, such as the Gauss-Seidel method and the Jacobi iteration method. Furthermore, since the system matrix A is positive definite and symmetric, methods such as the preconditioned conjugate gradient method, the sparse Cholesky decomposition method, and the algebraic multigrid method can also be employed.

[0068] 1. The linear system solution method using aggregation-based algebraic multigrids described in step 6 mainly consists of a setup phase and a solution phase. The setup phase is as follows: Setting up the stage sub-step 1: Filter the matrix, remove diagonal terms and off-diagonal terms of weak connections, and obtain the strength matrix S.

[0069] Setting up sub-step 2: The connection strength between nodes is determined by the off-diagonal values ​​in matrix S, and nodes with strong mutual influence are aggregated together. This results in non-overlapping and non-omitted aggregate partitions. Each aggregate partition forms a column of the aggregate matrix Agg. The Agg matrix is ​​a 0 / 1 matrix representing the sparse structure.

[0070] Step 3 of the setup phase: For each component of the near kernel vector B, insert it one by one into the corresponding position in the Agg matrix according to the aggregation partition represented in Agg. This forms a temporary extension matrix T. Divide T into blocks according to the aggregation, and use QR decomposition on each block. The decomposed Q and R are used as T and B for the next level, respectively.

[0071] Set up stage sub-step 4: Smooth the temporary extension operator T to obtain the extension operator P from the current layer to a coarser layer.

[0072] Set up stage sub-step 5: If the current level is the coarsest, stop; otherwise, proceed to the next level and return to step 2.

[0073] The solution process is as follows: Solution stage sub-step 1: Smooth the top layer using an iterative method and calculate the residual of the top layer linear system. The residuals are then passed to the next layer via the constraint operator. The limiting operator is the transpose of the extended operator. A new linear system is formed in the next layer. And the residual of the previous layer. As the new right-hand vector Smooth again, calculate the residuals, and propagate the residuals, repeating this process until the maximum number of layers L is reached. Obtain the residuals for all layers. .

[0074] Solution stage sub-step 2: Solve the lowest level residual linear system using the direct method. Where e represents the linear systematic error.

[0075] Solution stage sub-step 3: Convert the lowest level error Through the extension operator Pass to the penultimate level And the residuals of the penultimate layer are corrected using error correction. Then, for the penultimate layer of the linear system... Perform smoothing. Repeat this process until the top layer. Obtain the corrected solution vectors for all layers. .

[0076] Solution stage sub-step 4: Repeat steps 1 to 3 until the residual meets the accuracy requirements.

[0077] Step 7: Calculate the position correction based on the Lagrange multiplier increments, and correct the Lagrange multipliers and position. Position correction amount Calculate using the following formula:

[0078] This formula can be computed using a matrix-free method, resulting in better performance. The matrix multiplications are checked sequentially from right to left. It is the product of a 3n x m matrix and an m x 1 vector. (Vector) Each component corresponds to a spring. Each row corresponds to the gradient of a vertex, and each column corresponds to the gradient of each spring.

[0079]

[0080] A vertex may be connected to multiple springs simultaneously, and therefore has multiple gradients. Therefore, the matrix-vector multiplication is performed by multiplying all gradients of a given vertex (corresponding to different springs) by the Lagrange multiplier increments corresponding to all springs. The result is a vector, where each component corresponds to a vertex. Therefore, it is sufficient to sum the products over each vertex.

[0081]

[0082] Here, 'a' corresponds to different vertices.

[0083] because It is a diagonal matrix, and left-multiplying a diagonal matrix by any matrix is ​​equivalent to scaling that matrix by its rows. For example

[0084] Therefore, the formula for calculating without a matrix is:

[0085] Based on the Lagrange multiplier correction and the position correction, further corrections are made to the Lagrange multiplier and the position.

[0086]

[0087] Step 8: Determine whether the convergence accuracy has been reached based on the dual residual value. If it has, exit the iteration. The dual residual value is exactly the same as the right-hand vector b.

[0088] Here, the constraint vector C is an m x 1 vector. Lagrange multiplier vectors. It is also an mx1 vector. It is softness Divide by the square of the time step, i.e.

[0089] Among them, flexibility and time step All are specified by the user. A reasonable default value for flexibility is... The default reasonable value for the time step is s.

[0090] Step 9: Calculate the speed based on the difference between the new and old positions.

[0091] This invention achieves a faster reduction in residuals compared to XPBD and solves simulation defects such as excessive fabric stretching and overly soft softness. Figure 1 The overall process of this method is demonstrated. Figure 4 The dual residuals and energy convergence curves of this method in cloth simulation are shown. Figure 5 The dual residuals and energy convergence curves of XPBD in the same cloth simulation problem are shown. A comparison reveals that the number of iterations required for this method to achieve the same convergence accuracy is two orders of magnitude less than that of XPBD, resulting in a faster reduction in residuals. Figure 6The image shows a comparison of the results of this method and XPBD in cloth simulation. The left side shows the result of XPBD, and the right side shows the result of this method. It is evident that this method solves the overstretching problem in cloth simulation better than XPBD.

[0092] Figure 7 The dual residuals and energy convergence curves of this method in software simulation are shown. Figure 8 The dual residuals and energy convergence curves of XPBD for the same software simulation problem are shown. Similarly, a comparison shows that the number of iterations required for our method to achieve the same convergence accuracy is about 20 times less than that of XPBD, resulting in a faster decrease in residuals. Figure 9 The image shows a comparison of the results of this method and XPBD in software simulation. The left side shows the result of XPBD, and the right side shows the result of this method. It is evident that this method solves the overly soft problem in software simulation better than XPBD.

[0093] The above embodiments are merely preferred technical solutions of the present invention and should not be considered as limitations on the present invention. The scope of protection of the present invention should be limited to the technical solutions described in the claims, including equivalent substitutions of the technical features described in the claims. That is, equivalent substitutions and improvements within this scope are also within the scope of protection of the present invention.

Claims

1. An improved position dynamics simulation method based on algebraic multigrids, characterized in that, The solution process for each simulation substep includes the following steps: S1. Use the semi-implicit Euler method to calculate the predicted vertex position and record the old vertex position at the previous moment; S2. Initialize all Lagrange multiplier vectors corresponding to constraints to zero vectors; S3. Determine whether the current iteration result has reached the preset convergence accuracy. If it has, jump to step S8; otherwise, continue to execute the subsequent steps. S4. Calculate the constraint values ​​and constraint gradients of all constraints, and generate constraint vectors and constraint gradient matrices; S5. Assemble the global system matrix and right-hand vector based on the constraint gradient matrix, the mass diagonal matrix and the compliance parameter. The system matrix retains the off-diagonal terms corresponding to the adjacency relationship between constraints. S6. Solve the linear equations corresponding to the system matrix using the algebraic multigrid method to obtain the Lagrange multiplier increments; S7. Calculate the vertex position correction based on the Lagrange multiplier increment, update the Lagrange multipliers and vertex position synchronously, and return to step S3 to enter the next iteration. S8. Calculate the vertex velocity based on the difference between the updated vertex position and the old position, and complete the calculation of the current simulation sub-step.

2. The improved position dynamics simulation method based on algebraic multigrids according to claim 1, characterized in that, In step S5, the system matrix and the right-hand vector are derived through constraint gradients, quality parameters, and compliance parameters, and the calculation formulas are as follows: ； ； In the formula: The global system matrix has dimensions of . ; The constraint gradient matrix has dimensions of . ; It is a mass diagonal matrix with dimension . ; For mass diagonal matrix The inverse matrix; Constraint gradient matrix The transpose of the matrix; It is a flexibility diagonal matrix with dimension . ; Let be the right-hand vector of the system, with dimension . ; Let be a constraint value vector with dimension . ; Let be a Lagrange multiplier vector with dimension . ; This represents the total number of constraints in the simulation scenario. This represents the total number of vertices in the simulation scene.

3. The improved position dynamics simulation method based on algebraic multigrids according to claim 2, characterized in that, Step S5 employs a matrix-free method to rapidly assemble the system matrix elements. Leveraging the sparse structure of the system matrix, element calculations can be completed without explicitly storing the complete matrix. This matrix-free method is executed according to the discrete structure type of the simulation object, as detailed below: For a spring-mass point cloth system, each row of the system matrix corresponds to a spring constraint; the diagonal elements of the row are the sum of the reciprocals of the masses of the two vertices of the spring and the corresponding flexibility parameter; the off-diagonal elements of the row are calculated from the gradient of the adjacent springs that share vertices with the spring and the reciprocal of the mass of the shared vertices, and the number of off-diagonal terms is the same as the number of adjacent springs of the spring. For a tetrahedral soft system, each row of the system matrix corresponds to a tetrahedral constraint; the diagonal elements of the row are the sum of the products of the inverse mass of the four vertices of the tetrahedron and the square of the corresponding gradient; the off-diagonal elements of the row are calculated from the gradients of the adjacent tetrahedrons that share vertices with the tetrahedron and the inverse mass of the shared vertices, and the number of off-diagonal terms is the same as the number of adjacent tetrahedrons of the tetrahedron.

4. The improved position dynamics simulation method based on algebraic multigrids according to claim 1, characterized in that, Step S6 employs an aggregation-based algebraic multigrid method to solve the linear equation system. The setup phase of the algebraic multigrid method includes the following sub-steps: S6.1 Perform connection strength filtering on the system matrix, remove diagonal terms and weakly connected off-diagonal terms, and generate a strength matrix; S6.2 Determine the connection strength between nodes based on the off-diagonal values ​​of the strength matrix, aggregate nodes with strong mutual influence into non-overlapping and complete aggregate partitions, and generate an aggregate matrix of 0 / 1 form; S6.

3. Insert each component of the near-kernel vector into the corresponding position of the aggregation matrix according to the aggregation partition to form a temporary extension matrix; perform QR decomposition on the temporary extension matrix according to the aggregation block, and use the decomposition result as the temporary extension matrix and near-kernel vector of the next level. S6.4 Smooth the temporary extension matrix to obtain the extension operator from the current level to the next level; S6.5 Determine if the coarsest level has been reached. If not, proceed to the next level and repeat steps S6.2 to S6.4 until the operator construction for all levels is completed.

5. The improved position dynamics simulation method based on algebraic multigrids according to claim 4, characterized in that, The solution phase of the aggregation-based algebraic multigrid method includes the following sub-steps: S6.6 Perform iterative smoothing on the top-level linear system, calculate the residual of the current layer, and pass the residual to the next layer through the constraint operator to construct the next layer linear system; repeat the smoothing, residual calculation and passing operations until the bottom layer is reached to obtain the residuals of all levels; the constraint operator is the transpose of the corresponding extension operator. S6.

7. Solve the lowest level residual linear system using the direct method to obtain the lowest level error vector; S6.

8. Pass the bottom-level error vector to the next level through the extension operator, correct the residuals of the corresponding level and perform a smoothing operation; pass and correct it layer by layer until it returns to the top level to obtain the corrected solution vector. S6.9 Repeat steps S6.6 to S6.8 for the V loop operation until the residual meets the preset accuracy requirements, and output the final Lagrange multiplier increment.

6. The improved position dynamics simulation method based on algebraic multigrids according to claim 1, characterized in that, In step S3, the dual residual value is used as the criterion for judging the convergence accuracy of the iteration. The formula for calculating the dual residual value is as follows: ； In the formula: The dual residual value is used to characterize the convergence of the iterative solution; A vector of constraint values; This is the ratio of the compliance parameter to the square of the time step, where the compliance parameter and the time step are preset by the simulation scenario; is a Lagrange multiplier vector.

7. The improved position dynamics simulation method based on algebraic multigrids according to claim 1, characterized in that, In step S7, the vertex position correction is derived using the Lagrange multiplier increment and the constraint gradient, and the calculation formula is as follows: ； In the formula: The vertex position correction vector has a dimension of . ; It is the inverse matrix of the mass diagonal matrix; This is the transpose of the constraint gradient matrix; The Lagrange multiplier increment vector has dimension . ; The calculation uses a matrix-free method, which accumulates the product of the gradient of the corresponding constraint and the Lagrange multiplier increment along the vertex dimension, and then multiplies it by the reciprocal of the vertex quality to obtain the single vertex position correction.

8. An improved position dynamics simulation system based on algebraic multigrids, characterized in that, The improved position dynamics simulation method based on algebraic multigrids as described in any one of claims 1-7 is adopted. The system comprises a prediction initialization module, a constraint calculation module, a matrix assembly module, a multigrid solving module, an iterative update module, and a velocity calculation module, which are sequentially connected in communication: wherein: The prediction initialization module is used to calculate the predicted vertex positions using the semi-implicit Euler method, record the old positions, and initialize the Lagrange multipliers. The constraint calculation module is used to calculate the constraint values ​​and constraint gradients of all constraints, and generate constraint vectors and constraint gradient matrices. The matrix assembly module is used to assemble the global system matrix and the right-hand vector while retaining adjacent off-diagonal terms; The multigrid solution module is used to solve linear equation systems using the algebraic multigrid method to obtain the Lagrange multiplier increments. The iterative update module is used to calculate the vertex position correction, update the Lagrange multipliers and vertex positions, and determine the convergence of the iteration. The velocity calculation module is used to calculate the vertex velocity based on the difference between the old and new positions after the iteration converges.

9. The improved position dynamics simulation system based on algebraic multigrids according to claim 8, characterized in that, The multigrid solving module has built-in hierarchical building units and multilevel solving units. The hierarchical building units are used to generate multi-scale extension operators and constraint operators to complete the setup stage of algebraic multigrids. The multilevel solving units are used to execute the V-loop solving process of residual propagation, low-level solving and error correction to complete the solution stage of algebraic multigrids.

10. A computer-readable storage medium, characterized in that, A computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the improved position dynamics simulation method based on algebraic multigrids as described in any one of claims 1 to 7.

Images