An efficient automatic multi-substructure out-of-core parallel computing method

Abstract

The present application relates to an efficient automatic multi-substructure out-of-core parallel computing method, which divides a finite element simulation analysis model into a binary elimination tree composed of substructures; initializes the memory space; calculates the memory size required for node assembly, compares it with the allocated assembly memory, and updates the assembly memory; completes the matrix assembly of the node based on the updated assembly memory and updated data; writes the constraint matrix and eigenvector obtained by the solution to the disk as the first type of data; calculates the update data of the node to its ancestor node based on the constraint matrix and eigenvector solved in the node, and writes it to the disk as the second type of data; solves the generalized characteristic equation of dimensionality reduction in each node to obtain the characteristic pair; and obtains the formation vector based on the characteristic pair of each node and the constraint matrix back-substitution calculation. The method improves the memory utilization rate of the computer, and avoids performance loss as much as possible while reducing the peak memory in the calculation process.

Description

Technical Field

[0001] The present application relates to the technical field of out-of-core parallel computing, and in particular to an efficient automatic multiple substructure out-of-core parallel computing method. Background Art

[0002] The automatic multiple substructure method has a wide range of applications in engineering fields such as structural dynamics analysis, structural optimization design, vibration control, electromagnetics, fluid-solid interaction problems, etc., and can provide efficient and accurate solutions for dynamic analysis and optimization design. However, as the complexity of structural design increases, the accuracy and number of nodes in the model grid also increase, which leads to the huge scale of simulation problems. The memory capacity of a single computer is increasingly unable to meet the needs of actual calculations. At this time, neither expanding the number of memory bars nor using cluster computing is an economical option. In comparison, low-priced and high-capacity hard disks are undoubtedly more cost-effective, but there is currently no technical solution to achieve out-of-core parallelism of this algorithm.

[0003] The introduction and improvement of the automatic multiple substructure solving algorithm in the existing solutions are all based on in-core multithreading, aiming to achieve better load balancing or higher calculation accuracy. These technologies will inevitably encounter the problem of memory limitation in large-scale simulation calculations, and it is difficult to meet the actual needs of large-scale model simulation calculations on a single computer. Summary of the invention

[0004] Based on this, it is necessary to provide an efficient automatic multi-substructure out-of-core parallel computing method, which includes:

[0005] S1: Divide the finite element simulation analysis model into a binary elimination tree composed of substructures; and initialize the memory space to allocate memory for the nodes in the binary elimination tree;

[0006] S2: Calculate the memory size required for node assembly, compare it with the allocated assembly memory, and update the assembly memory; read update data from the disk according to the type of the node, and complete the matrix assembly of the node based on the updated assembly memory and the update data;

[0007] S3: Solve the constraint matrix and generalized characteristic equation in the node, and write the solved constraint matrix and characteristic vector to the disk as the first type of data; calculate the update data of the node to its ancestor node based on the constraint matrix and characteristic vector solved in the node; and write the updated data to the disk as the second type of data;

[0008] S4: Release the degree of freedom update position memory, assembly memory, and thread calculation memory allocated for the node;

[0009] S5: traversing the binary elimination tree, solving the reduced-dimensional generalized characteristic equation in each node, obtaining the characteristic pair, and releasing the allocated global memory;

[0010] S6: Based on the feature pairs of each node and the constraint matrix back-substitution calculation, the formation vector is obtained, the file stream of the thread is closed and the allocated constraint matrix memory is released.

[0011] Beneficial effects: This method solves the problem of insufficient memory on a single computer due to excessive data volume during large-scale model calculations, improves computer memory utilization, and avoids performance loss as much as possible while reducing peak memory during the calculation process. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0013] Figure 1 This is a flowchart of the automatic multi-substructure out-of-core parallel computing method according to an embodiment of the present application.

[0014] Figure 2 A schematic diagram of a substructure binary elimination tree according to an embodiment of the present application.

[0015] Figure 3 This is a flowchart of the out-of-core parallel substructure calculation of an embodiment of the present application.

[0016] Figure 4 A flowchart of node assembly according to an embodiment of the present application.

[0017] Figure 5 A schematic diagram of a structural subtree according to an embodiment of the present application.

[0018] Figure 6 This is a flowchart of node calculation according to an embodiment of the present application. DETAILED DESCRIPTION

[0019] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the specific implementation methods of the present application are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present application. However, the present application can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present application, so the present application is not limited by the specific embodiments disclosed below.

[0020] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In the description of this application, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0021] This embodiment divides and analyzes the calculation process of the substructure algorithm in detail, and adopts targeted processing strategies for the various intermediate data generated according to the subsequent use time nodes and methods. At the same time, in order to achieve efficient multi-threaded parallelism and avoid memory fragmentation caused by frequent allocation and release of memory, a complete set of data structures is designed as a buffer between memory and disk.

[0022] like Figure 1 As shown, this embodiment provides an efficient automatic multi-substructure out-of-core parallel computing method, the method comprising:

[0023] S1: Divide the finite element simulation analysis model into a binary elimination tree composed of substructures; and initialize the memory space to allocate memory for the nodes in the binary elimination tree.

[0024] In this embodiment, the structure of the substructure binary elimination tree is as follows Figure 2 As shown; the nodes in the binary elimination tree are divided into three categories: leaf nodes, nodes whose child nodes are leaf nodes, and other nodes. Different node assembly modes and corresponding node update data transmission modes are designed.

[0025] The partitioning method adopts the binary nested partitioning algorithm in computer graphics. Currently, the METIS software package is generally used in the field to complete the partitioning.

[0026] Specifically, the initializing memory space to allocate memory for nodes in the binary elimination tree includes:

[0027] Allocate a block of size N dof The global memory, N dof Represents the overall degrees of freedom of the finite element simulation analysis model; the global memory is used to store an array composed of the eigenvalues;

[0028] Allocate a second global memory for all non-root nodes in the binary elimination tree, wherein the second global memory includes a plurality of degree of freedom update position memories, each of which is used to record the degree of freedom update position of each non-root node to all its ancestor nodes; Figure 2 Node 1 in needs to determine which degrees of freedom of nodes 3, 7, and 15 will be updated.

[0029] Calculate the size of the constraint matrix of each node according to the number of degrees of freedom of each node, and select the largest constraint matrix size S max , allocate a block of size S t The third global memory; S t =N t *S max , N t Indicates the number of CPU threads used; each CPU thread of each node in the third global memory occupies a block of size S max Constraint matrix memory;

[0030] A piece of the assembly memory (1M in size) and a piece of the thread calculation memory (500Kb in size) are allocated to each CPU thread participating in the calculation; both the assembly memory and the thread calculation memory can be released separately and grow dynamically.

[0031] In this embodiment, each thread has assembled memory for temporary data generated during node calculation. These memories grow dynamically during node calculation and are released after all nodes are calculated. No memory is allocated independently during single node calculation.

[0032] The overall process of steps S2 and S3 is the calculation process of the out-of-core parallel substructure, and its flow chart is as follows: Figure 3 shown.

[0033] S2: Calculate the memory size required for node assembly, compare it with the allocated assembly memory, and update the assembly memory; read update data from the disk according to the type of node, and complete the matrix assembly of the node based on the updated assembly memory and the update data.

[0034] The node assembly process is as follows Figure 4 As shown, the specific process includes the following:

[0035] Allocating corresponding space for each node before assembling the node and releasing it after the end can easily lead to memory fragmentation and performance loss, so first you need to calculate the memory size required for node assembly;

[0036] The calculation formula for the memory size required for node assembly is:

[0037] ;

[0038] in, Indicates the memory size required for node assembly; Indicates the degree of freedom of the i-th node; represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree; Indicates the number of degrees of freedom that node i generates for node j as recorded in the degree of freedom update position memory;

[0039] When the node's assembled memory size is less than , immediately release the current assembly memory and apply for a block of size The memory is used as the new assembled memory of the node.

[0040] Further, the updating assembly memory includes:

[0041] When the node's assembled memory size is less than , immediately release the current assembly memory and apply for a block of size The memory is used as the new assembled memory of the node.

[0042] Furthermore, reading update data from the disk according to the type of the node includes:

[0043] If the current node has two or more levels of child nodes, the update data from its left child node is read from the disk; the second update data calculated by the right child node of the current node will be immediately used by the current node;

[0044] For example: When a thread calculates Figure 5 When the structural subtree shown is used, it will be calculated in the subsequent order of 1, 2, 3, 4, 5, 6, and 7. We notice that the updated data calculated by the right child node will be used by the parent node immediately, such as substructure 3 will be calculated immediately after substructure 2 is calculated. Therefore, the updated data of the right child node will not be written to the disk but will be released immediately after the parent node is used. On the other hand, the updated data of the left child, which is a leaf node, is only one leaf node away from the parent node, such as substructure 1 and substructure 3, and substructure 4 and substructure 6. Considering that the time in between is extremely short and the updated data of leaf nodes is generally small, the updated data of such nodes will not be written to the disk either.

[0045] In the updated assembly memory, the matrix assembly of the current node is completed based on the updated data.

[0046] In this embodiment, the assembled matrix includes: a stiffness matrix, a mass matrix, a coupling stiffness matrix and a coupling mass matrix; the calculation formulas of the stiffness matrix and the coupling stiffness matrix are as follows:

[0047] ;

[0048] , ;

[0049] in, represents the stiffness matrix of the i-th node; represents the coupling stiffness matrix between node i and node j; Represents the original data of the i-th node; Represents the updated data of the left child node of the i-th node; Represents the updated data of the right child node of the i-th node; represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree; Represents the left child node pair matrix of the i-th node Updated data; Represents the right child node pair matrix of the i-th node Updated data;

[0050] If the i-th node is a leaf node, then , All are 0;

[0051] The calculation process of the mass matrix is ​​consistent with that of the stiffness matrix;

[0052] The calculation process of the coupling mass matrix is ​​consistent with that of the coupling stiffness matrix.

[0053] S3: Solve the constraint matrix and generalized characteristic equation in the node, and write the solved constraint matrix and characteristic vector to the disk as the first type of data; calculate the update data of the node to its ancestor node based on the constraint matrix and characteristic vector solved in the node; and write the updated data to the disk as the second type of data.

[0054] The constraint matrix in the node is calculated based on the stiffness matrix and the coupled stiffness matrix. The calculation formula is:

[0055] ;

[0056] in, Represents the constraint matrix of node i on node j; represents the stiffness matrix of the i-th node; Represents the coupling stiffness matrix between node i and node j.

[0057] Calculating the updated data includes calculating the updated coupling stiffness matrix and coupling mass matrix for each ancestor node j of node i; the calculation formula is:

[0058] ;

[0059] ;

[0060] in, represents the coupling mass matrix between node j and node k after update; represents the coupling mass matrix between node j and node k; represents the constraint matrix of node i on node j; T represents transpose; Represents the constraint matrix of node i on node k; represents the coupling mass matrix between node i and node k; represents the coupling mass matrix between node i and node j; represents the coupling mass matrix between node i and node i; represents the updated coupling stiffness matrix between node j and node k; represents the coupling stiffness matrix between node j and node k; represents the coupling stiffness matrix between node i and node k; j represents the jth ancestor node of the ith node; k represents the node numbered k in the binary elimination tree.

[0061] like Figure 6 As shown, in the process of solving the constraint matrix and modal equation in the node, the constraint matrix in the node is written into the constraint matrix memory, and the size of the constraint matrix is ​​recorded; after the modal equation is calculated, if the equation has a solution, the obtained eigenvalue is written into the array in the global memory, and the eigennumber in the array will correspond to the position of the degree of freedom of the node. For example, if the current substructure calculates 10 eigenpairs, and its degree of freedom represents the 100th to 200th degrees of freedom of the global structure, the eigenvalue is written into the 100th to 109th position of the array, and the eigenvector size Svec is recorded (Svec = 0 when the equation has no eigensolution).

[0062] Afterwards, in order to reduce the number of disk files, the thread calculating the substructure will immediately merge the feature vector (if any) with the constraint matrix and write it into the thread file of the current thread. The size is Smat plus Svec, and the corresponding file name, data offset when writing and data size are recorded in the node.

[0063] When calculating the generalized characteristic equation of a node, if there is an eigensolution within the given calculation range, the eigenvalue is written to the global array opened up during the elimination tree construction phase, and the eigenvector is merged with the constraint matrix of the node and written to disk.

[0064] In this embodiment, after the update data is calculated, if the node is the left child node of its parent node and the node is a non-leaf node, a binary file identified by the node id is created in the current thread on disk, and the update data is written into the file.

[0065] Each thread creates a separate binary file on disk to store the node constraint matrix. The node constraint matrix data calculated by the same thread is written to the same file in an appended manner. At the same time, the file name, offset, and data size are recorded inside the node.

[0066] S4: Release the degree of freedom update position memory, assembly memory, and thread calculation memory allocated for the node.

[0067] Before calculating the dimension reduction generalized characteristic equation, release the memory of the global array, record the memory occupied by the substructure to update the position of the parent node, and assemble the memory.

[0068] S5: traverse the binary elimination tree, solve the reduced-dimensional generalized characteristic equation in each node, obtain the characteristic pair, and release the allocated global memory. The reduced-dimensional generalized characteristic equation is composed of the reduced-dimensional stiffness matrix and the reduced-dimensional mass matrix assembled from the node data.

[0069] Calculate the nodes in the entire binary elimination tree in the order of post-order traversal, read the corresponding data according to the offset recorded in the node and the number of bytes written to the disk, and release the thread constraint space memory and close the file stream after the calculation is completed.

[0070] Specifically, solving the reduced-dimensional generalized characteristic equation, the characteristic solutions include:

[0071] The calculation formula of the reduced-dimensional generalized characteristic equation is:

[0072] ;

[0073] in, represents the reduced-dimensional stiffness matrix; The quality matrix representing the reduced dimension; x represents the feature vector; represents the eigenvalue;

[0074] Solving the reduced-dimensional generalized characteristic equation is to calculate the reduced-dimensional generalized characteristic equation in a given interval A feature pair within, the feature pair comprising the feature vector and the feature value;

[0075] The diagonal elements in the reduced-dimensional stiffness matrix are the eigenvalues ​​stored in the array;

[0076] In this embodiment, the characteristic equation of dimensionality reduction adopts different solving methods according to the size of the dimension, as follows:

[0077] when When the dimension is not greater than 2000, the mass matrix of the node is extracted in the form of a dense matrix to obtain the reduced-dimensional mass matrix, and the LAPACKE_dgeev function of the Lapack module in OpenBLAS is used to solve all the eigenvalues ​​of the generalized characteristic equation, and then the eigenvalues ​​in the given interval are extracted. The feature pairs within;

[0078] when When the dimension is greater than 2000, the reduced-dimensional quality matrix is ​​assembled in the form of a CSR sparse matrix, and the generalized characteristic equation is solved using the Block Lanczos method, and the global memory is released after the solution.

[0079] S6: Based on the characteristic solution of each node and the constraint matrix back-substitution calculation, the formation vector is obtained, the file stream of the thread is closed and the allocated constraint matrix memory is released.

[0080] The calculation formula of the formation vector is:

[0081] , ;

[0082] in, represents the formation vector of the i-th node; The matrix representing the eigenvector of the i-th node solved by S3; represents the characteristic solution of the i-th node; Represents the constraint matrix of node i on node j; represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree.

[0083] In this embodiment, the difference from the parallel computing in the core is that the matrix and All are stored on disk, so for each child node, the data needs to be read according to the file name, data offset and size stored during the node calculation process. And each substructure will The data is read into the thread's constraint matrix memory to further avoid frequent memory allocation and release.

[0084] The calculated formation vector can be used to obtain the actual response of the structure under a specific load, thereby evaluating whether the dynamic characteristics of the structure meet the requirements. It can also be used to predict dangerous parts of the structure and verify the accuracy of the finite element model.

[0085] The efficient automatic multi-substructure out-of-core parallel computing method provided in this embodiment divides the finite element simulation analysis model into a plurality of substructures (and records the division order in the form of a binary tree); initializes the memory space; calculates the memory size required for node assembly, and compares it with the allocated assembly memory, and updates the assembly memory; completes the matrix assembly of the node based on the updated assembly memory and the updated data; writes the constraint matrix and the eigenvector obtained by the solution to the corresponding position of the disk as the first type of data in a specific order; calculates the updated data of its ancestor node based on the constraint matrix and the eigenpair obtained by the solution in the node; and writes the updated data to the disk as the second type of data; assembles the reduced-dimensional stiffness matrix and mass matrix to solve the generalized characteristic equation to obtain the reduced-dimensional characteristic solution; obtains the formation vector based on the characteristic solution of each node, the reduced-dimensional characteristic solution and the constraint matrix back-substitution calculation. This method improves the memory utilization rate of the computer, and avoids performance loss as much as possible while reducing the peak memory in the calculation process.

[0086] The method provided in this embodiment effectively solves the problem of excessive memory required by the automatic multiple substructure algorithm in the process of large-scale model CAE simulation calculation. Before calculating the substructure model conversion, the corresponding mapping array is constructed for the nodes in the elimination tree. During the calculation process, different data processing strategies are adopted for the assembly data, coupling update matrix, and internal matrix of the nodes, which significantly reduces the amount of data remaining in the memory and realizes an efficient out-of-core parallel calculation method for automatic multiple substructures of large-scale sparse generalized characteristic equations. At the same time, the present invention ensures that the calculation results are consistent with the parallel accuracy in the core, and can effectively meet the solution requirements of large-scale models in a number of fields such as finite element simulation analysis, computational fluid dynamics, and electromagnetism.

[0087] The technical features of the above-described embodiments may be arbitrarily combined. To make the description concise, not all possible combinations of the technical features in the above-described embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0088] The above-described embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be construed as limiting the scope of the patent application. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent application shall be subject to the attached claims.

Claims

1. An efficient automatic multi-substructure out-of-core parallel computing method, characterized in that: include: S1: Divide the finite element simulation analysis model into a binary elimination tree consisting of substructures; And initialize the memory space and allocate memory for the nodes in the binary elimination tree; Initializing the memory space to allocate memory for the nodes in the binary elimination tree includes: Allocate a block of size N dof Global memory, N dof Represents the overall degree of freedom of the finite element simulation analysis model; the global memory is used to store an array composed of eigenvalues; Allocating a second global memory for all non-root nodes in the binary elimination tree, wherein the second global memory includes a plurality of degree-of-freedom update position memories, each of which is used to record the degree-of-freedom update position of each non-root node to all its ancestor nodes; Calculate the size of the constraint matrix of each node according to the number of degrees of freedom of each node, and select the largest constraint matrix size S max , allocate a block of size S t The third global memory; S t =N t *S max , N t Indicates the number of CPU threads used; each CPU thread of each node in the third global memory occupies a block of size S max Constraint matrix memory; Allocate a piece of assembly memory and a piece of thread calculation memory to each CPU thread involved in the calculation; both the assembly memory and the thread calculation memory can be released separately and grow dynamically; S2: Calculate the memory size required for node assembly, compare it with the allocated assembly memory, and update the assembly memory; read update data from the disk according to the type of the node, and complete the matrix assembly of the node based on the updated assembly memory and the update data; S3: Solve the constraint matrix and generalized characteristic equation in the node, and write the solved constraint matrix and characteristic vector to the disk as the first type of data; calculate the update data of the node to its ancestor node based on the constraint matrix and characteristic vector solved in the node; and write the updated data to the disk as the second type of data; S4: Release the degree of freedom update position memory, assembly memory, and thread calculation memory allocated for the node; S5: traverse the binary elimination tree, solve the reduced-dimensional generalized characteristic equation in each node, obtain characteristic pairs, and release the allocated global memory; S6: Based on the feature pairs of each node and the constraint matrix back-substitution calculation, the formation vector is obtained, the file stream of the thread is closed and the allocated constraint matrix memory is released; The calculation formula of the formation vector is: in, V represents the formation vector of the i-th node; i The matrix representing the eigenvector of the i-th node solved by S3; x i represents the eigenvector of the i-th node; Ψ ij A represents the constraint matrix of node i on node j; i represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree.

2. The automatic multi-substructure out-of-core parallel computing method according to claim 1, characterized in that: The calculation formula for the memory size required for node assembly is: Among them, M total Indicates the memory size required for node assembly; d i A represents the degree of freedom of the i-th node; i represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree; d i,j Indicates the number of degrees of freedom that node i generates for node j as recorded in the degree of freedom update position memory; When the size of the node's assembled memory is less than M total When the current assembly memory is released immediately, a block of size M is requested. total The memory is used as the new assembled memory of the node.

3. The automatic multi-substructure out-of-core parallel computing method according to claim 1, characterized in that: Reading update data from the disk according to the type of the node includes: If the current node has two or more levels of child nodes, the first update data from its left child node is read from disk; the second update data calculated by the right child node of the current node will be used immediately by the current node; In the updated assembly memory, the matrix assembly of the current node is completed based on the updated data.

4. The automatic multi-substructure out-of-core parallel computing method according to claim 1, characterized in that: In the process of solving the constraint matrix and modal equations in the node, the constraint matrix in the node is written into the constraint matrix memory, and the size of the constraint matrix is ​​recorded; after the modal equation is calculated, if the equation has a solution, the obtained eigenvalue is written into the array in the global memory, and the eigenvalue in the array will correspond to the degree of freedom position of the node.

5. The automatic multi-substructure out-of-core parallel computing method according to claim 1, characterized in that: After calculating and obtaining the update data, if the node is the left child node of its parent node and the node is a non-leaf node, the current thread creates a binary file in the disk using the id of the node as identification; Write the updated data into it.

6. The automatic multi-substructure out-of-core parallel computing method according to claim 4, characterized in that: Solving the reduced-dimensional generalized characteristic equation, the characteristic vectors obtained include: The calculation formula of the reduced-dimensional generalized characteristic equation is: in, represents the reduced-dimensional stiffness matrix; represents the quality matrix of dimensionality reduction; x represents the eigenvector; λ represents the eigenvalue; Solving the reduced-dimensional generalized characteristic equation is to calculate the reduced-dimensional generalized characteristic equation in a given interval [σ min ,σ max ] within a feature pair, the feature pair comprising the feature vector and the feature value; The diagonal elements in the reduced-dimensional stiffness matrix are the eigenvalues ​​stored in the array; when When the dimension is not greater than 2000, the mass matrix of the node is extracted in the form of a dense matrix to obtain the reduced-dimensional mass matrix, and the LAPACKE_dgeev function of the Lapack module in OpenBLAS is used to solve all the characteristic pairs of the generalized characteristic equation, and then the characteristic solution is extracted in the given interval [σ min ,σ max ]; when When the dimension is greater than 2000, the reduced-dimensional quality matrix is ​​assembled in the form of a CSR sparse matrix, and the generalized characteristic equation is solved using the BlockLanczos method, and the global memory is released after the solution.

7. The automatic multi-substructure out-of-core parallel computing method according to claim 6, characterized in that: The assembled matrices include: stiffness matrix, mass matrix, coupled stiffness matrix and coupled mass matrix; the calculation formulas of stiffness matrix and coupled stiffness matrix are as follows: Among them, K i represents the stiffness matrix of the i-th node; K i,j represents the coupling stiffness matrix between node i and node j; Represents the original data of the i-th node; Represents the left child node pair matrix K of the i-th node i Updated data; Represents the right child node pair matrix K of the i-th node i Update data of A i represents the set of all ancestor nodes of the i-th node; j represents the node numbered j in the binary elimination tree; Represents the left child node pair matrix K of the i-th node i,j Updated data; Represents the right child node pair matrix K of the i-th node i,j Updated data; If the i-th node is a leaf node, then All are 0; The calculation process of the mass matrix is ​​consistent with that of the stiffness matrix; The calculation process of the coupling mass matrix is ​​consistent with that of the coupling stiffness matrix.

Images