An isogeometric stiffness matrix cpu / gpu heterogeneous parallel computing method

Abstract

The application belongs to the technical field of isogeometric analysis, and discloses an isogeometric stiffness matrix CPU / GPU heterogeneous parallel computing method.The method comprises the following steps: S1, obtaining an isogeometric body grid model of an object to be processed, dividing the isogeometric body grid model to obtain a plurality of body units, control points and Gauss points, and numbering the body units, control points and Gauss points respectively; S2, obtaining corresponding information vector groups of all interactive control point pairs in all body units; S3, obtaining a shape function partial derivative matrix at each Gauss point; S4, calculating a global stiffness matrix of the interactive control point pairs, thereby obtaining global stiffness matrices of all interactive control point pairs, and assembling the global stiffness matrices of all interactive control point pairs, so as to realize heterogeneous parallel assembly of the stiffness matrix.Through the application, the problems of low calculation efficiency, large memory consumption and incapability of isogeometric analysis on large-scale high-resolution models in the existing isogeometric analysis method are solved.

Description

A CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix Technical Field

[0001] This invention belongs to the technical field of isogeometric analysis, and more specifically, relates to a CPU / GPU heterogeneous parallel computing method for isogeometric stiffness matrices. Background Technology

[0002] Isogeometric Analysis (IGA), as a novel finite element analysis method, aims to seamlessly integrate Computer-Aided Design (CAD) and Computer-Aided Engineering (CAE). In traditional CAD and CAE, models designed during the CAD process are typically represented using non-uniform rational B-splines (NURBS), while CAE analysis uses Finite Element Analysis (FEA). FEA uses a mesh as its geometric representation, meaning that the NURBS model needs to be converted to a mesh before FEA analysis. This conversion process has two drawbacks: first, it is extremely time-consuming, accounting for 80% of the entire FEA time; second, the conversion process uses a mesh to approximate the CAD model, introducing errors. To overcome these shortcomings, IGA directly uses NURBS as its geometric representation, eliminating the time-consuming conversion process and avoiding the errors introduced during conversion, thus achieving seamless integration of CAD and CAE.

[0003] Although IGA (Integrated Geometry Analysis) is promising compared to FEA (Functional Elements in Calculation), possessing advantages such as accurate geometric discretization, higher-order continuity, and unified data representation for CAD and CAE models, the computation of the stiffness matrix is ​​more laborious, increasing the cost of IGA in practical applications. This is because IGA has higher element continuity; for the same degrees of freedom, IGA generates significantly more elements than FEA. More elements lead to an increased number of Gaussian points, thus increasing the computational cost of assembling the stiffness matrix. This drawback significantly increases the computational cost of IGA, especially in 3D analysis. In recent years, many scholars have studied how to accelerate the computation of stiffness matrices for isogeometric analysis to achieve isogeometric analysis of large-scale high-resolution models. Parallelizing the computation process is one of the mainstream methods. Parallel computing is computation performed on parallel computers. Parallel algorithms divide the computational task into multiple subtasks, which are executed simultaneously by different GPUs and CPUs to achieve higher computational efficiency. Based on CPU / GPU heterogeneous architectures, existing parallel computation methods for stiffness matrices often neglect the parallel computing capabilities of CPUs. The CPU is responsible for data interaction with the GPU, and all computational tasks are undertaken by the GPU. This results in the ineffective utilization of local computing and memory resources, making the isogeometric analysis method unable to process large-scale high-resolution models.

[0004] Based on a CPU / GPU heterogeneous architecture, this method utilizes the collaborative parallel computation of the CPU and GPU to compute the isogeometric stiffness matrix. This approach effectively leverages local computing and memory resources, enabling isogeometric analysis methods to handle large-scale, high-resolution models. Therefore, designing a CPU / GPU heterogeneous parallel computation method for the isogeometric stiffness matrix is ​​of great significance for isogeometric analysis. Summary of the Invention

[0005] To address the aforementioned shortcomings or improvement needs of existing technologies, this invention provides a CPU / GPU heterogeneous parallel computing method for isogeometric stiffness matrices, solving the problems of low efficiency in stiffness matrix calculation, high memory consumption, and inability to perform isogeometric analysis on large-scale high-resolution models in existing isogeometric analysis methods.

[0006] To achieve the above objectives, according to the present invention, a CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix is ​​provided, the method comprising the following steps:

[0007] S1 performs isogeometric meshing on the object to be processed to obtain an isogeometric mesh model, divides the isogeometric mesh model to obtain multiple volume elements and corresponding control points and Gaussian points in each volume element, and numbers the volume elements, control points and Gaussian points respectively.

[0008] S2 For any pair of interactive control points in the body unit, the information vector group of the pair of interactive control points is formed by using the number of the body unit corresponding to the pair of interactive control points and the number of each control point, so as to obtain the information vector group corresponding to all pairs of interactive control points in all body units.

[0009] S3 performs heterogeneous parallel computation of the partial derivatives of the shape function at each Gaussian point of the volume element to obtain the partial derivative matrix of the shape function at each Gaussian point.

[0010] S4 uses the number of the corresponding volume element in the information vector group of the interactive control point pair to determine the shape function partial derivative matrix of the Gaussian point corresponding to the volume element of the interactive control point pair, and uses the shape function partial derivative matrix to calculate the strain matrix, thereby establishing the correspondence between the interactive control point and the strain matrix; uses the correspondence to calculate the global stiffness coefficient submatrix of the interactive control point pair, thereby obtaining the global stiffness coefficient submatrix of all interactive control point pairs, and assembling the global stiffness coefficient submatrix of all interactive control point pairs to realize the heterogeneous parallel assembly of the stiffness matrix.

[0011] More preferably, in step S3, the heterogeneous parallel computation of the partial derivatives of the shape functions is performed in the following manner:

[0012] S31 takes the partial derivative of the shape function at all Gaussian points in each volume element as a subtask, thus forming multiple subtasks;

[0013] S32 allocates the multiple tasks between the CPU and GPU and starts CPU and GPU multi-threading, wherein the number of CPU threads is the number of CPU cores, and one thread will be responsible for executing multiple sub-tasks; the number of GPU threads is the number of sub-tasks allocated to the GPU, and one thread will be responsible for executing one sub-task.

[0014] S33 describes the partial derivatives of shape functions calculated in parallel by multi-threaded CPU and GPU.

[0015] More preferably, in step S31, the two-dimensional shape function is determined according to the following relationship:

[0016]

[0017] in, It is a shape function in a two-dimensional parametric coordinate space, where p and q are NURBS orders, i and j are control point labels, ξ and η are coordinate values ​​in the parametric coordinate space, and W i It is a weighting factor with the property of the total tensor product: W i,j =W i W j n and m are the number of parameter coordinates. These are shape functions on coordinates ξ and η in a one-dimensional parameter space, respectively.

[0018] More preferably, in step S33, the partial derivatives of the shape functions are calculated according to the following relationship:

[0019]

[0020] in, It is a shape function in a two-dimensional parametric coordinate space, where p and q are NURBS orders, i and j are control point labels, ξ and η are coordinate values ​​in the parametric coordinate space, and W i It is a weighting factor with the property of the total tensor product, W i,j =W i W j n and m are the number of parameter coordinates. These are shape functions on coordinates ξ and η in a one-dimensional parameter space, respectively.

[0021] More preferably, in step S4, the calculation of the global stiffness matrix of the interactive control point pair is performed according to the following steps:

[0022] S41 divides the solution of the global stiffness submatrix of all interactive control point pairs into multiple subtasks, thereby obtaining multiple subtasks.

[0023] S42 distributes subtasks between the CPU and GPU and enables CPU and GPU multithreading. The CPU thread is responsible for calculating the stiffness coefficient submatrix of multiple interactive control point pairs; the GPU thread is only responsible for calculating the stiffness coefficient submatrix of one interactive control point pair.

[0024] The S43 CPU and GPU perform multi-threaded parallel computation of the global stiffness matrix.

[0025] More preferably, in step S42, the global stiffness matrix is ​​calculated in the following manner:

[0026] S431 determines the shared unit of the interactive control point pair;

[0027] S432 calculates the stiffness coefficient sub-element for each of the shared elements;

[0028] S433 constructs the relationship between the stiffness coefficient sub-element of each shared element and the global stiffness coefficient sub-matrix, thereby obtaining the global stiffness coefficient sub-matrix of the interactive control point pair;

[0029] S434 assembles the global stiffness coefficient submatrices of all interactive control points to obtain the global stiffness matrix of the object to be processed.

[0030] More preferably, in step S432, for an interactive control point pair ij, the stiffness coefficient submatrix H of the shared element to which the interactive control point pair belongs... ij Calculate according to the following formula:

[0031]

[0032] Where D is the elasticity matrix, B i B j These are the strain matrices corresponding to control points i and j, respectively, and strain matrix B. i and B j The elements in the equation are the partial derivatives of the shape function at the Gaussian point of the volume element, representing the geometric coordinates of the derivatives.

[0033] More preferably, in step S433, the global stiffness coefficient submatrix is ​​determined according to the following relationship:

[0034]

[0035] Among them, K ij It is the global stiffness coefficient submatrix, w G It is the weighting coefficient, H ij It is the stiffness coefficient submatrix.

[0036] More preferably, in step S434, the global stiffness matrix is ​​a sparse matrix, wherein the elements are indexed in the following manner:

[0037]

[0038] Where x and y are the coordinate directions of the horizontal and vertical axes in geometric space, respectively, and edof(i) and edof(j) represent the degree-of-freedom label vectors of control points i and j in the x and y directions. This represents the operation of converting the degree-of-freedom labels to row and column indices in the stiffness matrix K, iK ij It is a vector of row index numbers, jK ij It is a vector of column index numbers.

[0039] In summary, the technical solutions conceived by this invention have the following beneficial effects compared with the prior art:

[0040] 1. The present invention proposes a CPU / GPU heterogeneous parallel computing method for isogeometric stiffness matrices. This method assembles stiffness matrices based on interactive control point pairs. The designed parallel algorithm avoids memory write competition between parallel threads and ensures parallel computing efficiency without data synchronization. Threads directly calculate the non-zero stiffness coefficients at the corresponding positions of the stiffness matrix according to the control point pairs, making efficient heterogeneous parallel computing of stiffness matrices possible.

[0041] 2. The heterogeneous parallel algorithm proposed in this invention enables the CPU to start multiple threads to participate in parallel computing. Compared with the existing technology in which the CPU is only responsible for data interaction with the GPU, it makes full use of the local computing power and memory resources, making it possible to efficiently process large-scale high-resolution models in isogeometric analysis.

[0042] 3. This invention proposes a dynamic load balancing method for heterogeneous parallel computing. The algorithm is based on the general formula method but with improvements. This method adjusts the computing load on the CPU and GPU according to the real-time running parameters of the local machine, so that the computing time of the CPU and GPU is close to the same. At this time, the total computing time is the shortest and the heterogeneous parallel computing efficiency is the highest. Attached Figure Description

[0043] Figure 1 is a flowchart of a CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix constructed according to a preferred embodiment of the present invention;

[0044] Figure 2 is a flowchart of interactive control point pairs in the identification geometric model constructed according to a preferred embodiment of the present invention;

[0045] Figure 3 is a schematic diagram of the control point vectors of the volume element in the isogeometric model mesh topology constructed according to a preferred embodiment of the present invention;

[0046] Figure 4 is a schematic diagram of the relevant volume element vectors of the control points in the isogeometric model mesh topology constructed according to a preferred embodiment of the present invention;

[0047] Figure 5 is a schematic diagram of the search radius of the interactive control point under different NURBS basis function orders in the label space, constructed according to a preferred embodiment of the present invention;

[0048] Figure 6 is a schematic diagram of a CPU / GPU heterogeneous parallel computing model constructed according to a preferred embodiment of the present invention;

[0049] Figure 7 is a flowchart of heterogeneous parallel computation of partial derivatives of shape functions constructed according to a preferred embodiment of the present invention;

[0050] Figure 8 is a schematic diagram of a heterogeneous parallel computing thread model of partial derivatives of shape functions constructed according to a preferred embodiment of the present invention;

[0051] Figure 9 is a schematic diagram of the heterogeneous parallel assembly process of stiffness matrix constructed according to a preferred embodiment of the present invention;

[0052] Figure 10 is a schematic diagram of a heterogeneous parallel assembly thread model of stiffness matrix constructed according to a preferred embodiment of the present invention;

[0053] Figure 11 is a schematic diagram of the COO format used in the compressed storage sparse matrix K constructed according to a preferred embodiment of the present invention;

[0054] Figure 12 is a flowchart of a dynamic workload balancing algorithm for heterogeneous parallel computing constructed according to a preferred embodiment of the present invention. Detailed Implementation

[0055] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0056] A heterogeneous parallel computing method for CPU / GPU with equal geometric stiffness matrix is ​​shown in Figure 1. Its main steps include:

[0057] Step 1: Identify interactive control point pairs in the geometric model. Based on the volume element where the control point pair is located and the local position of the control point in the volume element, generate a control point pair information vector group.

[0058] In this invention, an interactive control point pair is defined as two control points that jointly participate in forming at least one isogeometric element, and the two control points may be the same. Multiple isogeometric elements formed by a pair of interactive control points are defined as shared elements, which interact within the shared elements, and the resulting element stiffness coefficients contribute to the formation of the global stiffness coefficient.

[0059] As shown in Figure 2, the identification of interactive control point pairs includes:

[0060] 1) Initialize the mesh topology information. As shown in Figure 3, traverse the volume elements sequentially according to the isogeometric mesh model, record the control points contained in the volume elements to generate the volume element control point vector; traverse the control points sequentially, and determine the isogeometric elements in which the control points participate based on the element control point vectors, generating the control point element vector as shown in Figure 4;

[0061] 2) Search for interactive control point pairs. Map control points to the label space based on global labels. In 3D and other geometric models, the mapping relationship between control point labels and coordinates can be expressed as:

[0062]

[0063] In the formula, x, y, z are the coordinates of the control point in the label space, idx is the global label of the control point, and n x n y n z Let be the total number of control points along each coordinate axis. When searching for the interactive control points of control point i in the label space, the search radius is determined by the degree of the NURBS basis function of the isogeometric model, as shown in Figure 5. When the NURBS degree p = 2, 3, and 4, the search radii on the coordinate axes are 2, 3, and 4 respectively. Therefore, the coordinate relationship between interactive control points i and j can be expressed as:

[0064]

[0065] 3) Generate control point pair information vectors. In this invention, the information vector set for control point pairs ij includes two types: ① Shared unit vectors. In the control point unit vectors, the shared unit vectors are generated by searching for the units jointly formed by control points i and j; ② Local label vectors. In the control point unit vectors, the local labels of control points i and j within the shared units are searched.

[0066] Step 2 employs a heterogeneous parallel algorithm for the partial derivatives of shape functions to perform heterogeneous parallel computation on the partial derivatives of shape functions at Gaussian points of equal geometric units, in order to generate the partial derivative matrix of shape functions.

[0067] In this invention, heterogeneous parallel computing is defined as follows: CPU and GPU, based on a heterogeneous computer architecture, collaboratively solve computational problems using multi-threading technology. This invention employs CPU / GPU heterogeneous parallel computing based on OpenMP (Open Multi-Processing) + CUDA (Compute Unified Device Architecture). As shown in Figure 6, OpenMP is used to enable multi-threading on the CPU, where the main thread is responsible for interacting with the GPU for data, and other threads are responsible for computation. The CPU and GPU work collaboratively simultaneously under specific workloads, with the total workload distributed across both the CPU and GPU. The CPU operates in a "single-threaded multi-node" mode, where each thread iterates through multiple tasks in a loop; the GPU operates in a "single-threaded single-node" mode, where each thread executes only one task.

[0068] As shown in Figure 7, the heterogeneous parallel computation of shape function partial derivatives includes: 1) Assigning tasks and enabling CPU and GPU multithreading. As shown in Figure 8, using isogeometric units as units, the total task is divided into subtasks of the number of units using a divide-and-conquer approach. Each subtask solves for the partial derivatives of the shape function at all Gaussian points in each unit. Subtasks are assigned between the CPU and GPU, and CPU and GPU multithreading is enabled. The number of CPU threads is equal to the number of CPU cores, with one thread responsible for executing multiple subtasks; the number of GPU threads is equal to the number of subtasks allocated to the GPU, with one thread responsible for executing one subtask; 2) Parallel computation of shape function partial derivatives using CPU and GPU multithreading. In this invention, the partial derivative of the shape function is defined as the derivative of the shape function at a Gaussian point of the unit with respect to the parameter coordinates. In isogeometric analysis, NURBS curve basis functions are used as shape functions. In the 2D example, the shape function can be expressed as:

[0069]

[0070] Among them, W i These are weighting factors that possess the properties of the total tensor product:

[0071] W i,j =W i W j (4)

[0072] According to the above formula, the derivative of the shape function with respect to the parametric coordinate ξ can be expressed as:

[0073]

[0074] The CPU and GPU threads will traverse all Gaussian points in the cells corresponding to their subtasks, calculate the shape functions based on the parameter coordinates of the Gaussian points, and then sequentially complete the differentiation of the shape functions with respect to each parameter coordinate, finally completing the partial derivative calculation of the shape functions for all cells.

[0075] Step 3: Based on steps 1 and 2, calculate the stiffness coefficient matrix between control point pairs and its index vector in the global stiffness matrix, completing the heterogeneous parallel assembly of the stiffness matrix. The shared element vector and local label vector generated in step 1, as well as the shape function partial derivative matrix generated in step 2, will be stored in the GPU global memory and CPU shared memory as the raw material for calculating the stiffness coefficient matrix in step 3. As shown in Figure 9, the heterogeneous parallel assembly method for stiffness matrix has the following steps:

[0076] 1) Assign tasks and enable CPU and GPU multithreading. Divide the overall task into control point pairs, and distribute subtasks between the CPU and GPU, enabling CPU and GPU multithreading. As shown in Figure 10, the CPU thread will be responsible for calculating the stiffness coefficients of multiple interactive control point pairs; the GPU thread will only be responsible for calculating the stiffness coefficients of one interactive control point pair.

[0077] 2) The stiffness matrix can be assembled in parallel using CPU and GPU multi-threading, which can be achieved through the following steps:

[0078] ① Calculate the global stiffness coefficient submatrix based on the interactive control point pairs ij. ij The element stiffness coefficient submatrix for each control point pair ij in each shared element needs to be calculated sequentially, and the calculation formula is as follows:

[0079]

[0080] Where D is the elasticity matrix, B i B j This is the strain matrix corresponding to the control point pair ij. Strain matrix B i B j The elements in the matrix are the partial derivatives of the shape function in geometric coordinates at the Gaussian point of the element, which can be generated by calculating the partial derivatives of the shape function in parametric coordinates and the Jacobian matrix. The stiffness coefficient submatrix H of all shared elements is then used. ij After multiplying by the weighting coefficients, the results are summed to generate the final global stiffness coefficient submatrix for the interactive control points, as shown in the following formula:

[0081]

[0082] Specifically, the CPU and GPU threads determine all shared elements for the control point pair ij based on the shared element vector, and calculate the element stiffness submatrix H using the local label vector and the partial derivative matrix of the element Gaussian point shape function. ij The summation produces the global stiffness coefficient submatrix K. ij When calculating the global stiffness coefficient submatrix based on control point pairs, the elements in the submatrix are the final values ​​and can be stored in different memory addresses. This avoids memory contention and ensures that threads perform safe concurrent read and write operations on global video memory or shared memory throughout the entire calculation process.

[0083] ② Calculate the stiffness coefficient index vector. In the global stiffness matrix K, the number of non-zero elements is much smaller than the number of zero elements, making the matrix sparse. To save memory resources, this invention uses the COO (Coordinate) format to compress and store the sparse matrix K. This compression format only records the information of non-zero elements, thus saving memory overhead for storing matrix K. As shown in Figure 11, the COO format uses C / C++ arrays to store three vectors: row and column index vectors (iK, jK) and a vector of non-zero element values ​​(vK). The global stiffness coefficient submatrix K for control point pairs ij. ij K is a constituent unit of vector vK. ij The row and column index pairs corresponding to the elements in matrix K are unique and non-repeating. ij The index of an element is determined by the degree of freedom labels of control points i and j. In the 2D example, the index can be calculated as follows:

[0084]

[0085] Where edof(i) and edof(j) represent the degree-of-freedom label vectors of control points i and j in the x and y directions, respectively. This represents the operation of converting the degree-of-freedom labels to row and column indices in the stiffness matrix K. Concurrent threads perform this operation after calculating all submatrices K. ij After obtaining the corresponding index vector, the assembly of the sparse matrix K can be completed.

[0086] In steps 2 and 3 above, balancing the workload between the CPU and GPU is crucial for achieving efficient computing in heterogeneous parallel computing. If the task allocation scheme is determined directly by quantifying hardware computing power, although this power can be obtained from APIs, the actual computing efficiency will be affected by the parallel algorithm and hardware operation, resulting in the failure to achieve workload balancing. Therefore, this invention uses a dynamic workload balancing method, combining pre-running with theoretical analysis to obtain the workload balancing scheme.

[0087] Assume there are N x ×N y ×N z For each independent task, the CPU time to complete the task can be expressed as:

[0088] τ CPU =(1-α)(N) x ×N y ×N z )t CPU (9)

[0089] Among them, t CPUα is the time required for the CPU to execute a subtask, and α is the workload allocation rate. Considering the time consumed by data transfer between the CPU and GPU, the time for the GPU to complete a task can be expressed as:

[0090] τ GPU =τ DT +τ G (10)

[0091] Where, τ DT It is the time for data transmission, τ G This is the time spent on GPU computation. Given a workload allocation rate α, τ DT and τ GPU It can be represented as:

[0092]

[0093] Where k represents the space complexity coefficient, S val t is the number of bytes per data unit, and v is the bandwidth capacity of the PCI-E bus connecting the CPU and GPU for data transfer. dt t represents the average data transfer time of a task. GPU This represents the computation time for the GPU to execute a task. When τ CPU =τ GPU When the total computation time τ is minimized, as shown below:

[0094] (1-α)(N x ×N y ×N z )t CPU =α(N) x ×N y ×N z )(t dt +t GPU (12)

[0095] Therefore, the workload allocation rate α can be expressed as:

[0096]

[0097] In the dynamic approach, the pre-run phase aims to obtain the actual data transmission time t. dt and calculation time t GPU t CPU As shown in Figure 12. Workload during the pre-operation phase Π pre The number of subtasks is greater than (n-1), ensuring that each CPU core has a subtask to execute. After the pre-running phase, the actual execution times Γ1 and Γ2 of the CPU and GPU will be obtained. Therefore, the data transfer time and computation time can be evaluated as follows:

[0098]

[0099] Finally, the workload allocation rate α can be expressed as follows:

[0100]

[0101] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix, characterized in that, The method includes the following steps: S1. The object to be processed is meshed into an isogeometric mesh model. The isogeometric mesh model is then divided to obtain multiple volume elements and corresponding control points and Gaussian points in each volume element. The volume elements, control points, and Gaussian points are numbered respectively. S2. For any pair of interactive control points in the volume element, an information vector group for the pair of interactive control points is formed using the number of the volume element corresponding to the pair of interactive control points and the number of each control point. This yields the information vector groups corresponding to all pairs of interactive control points in all volume elements. S3. The partial derivative of the shape function at each Gaussian point of the volume element is calculated. Heterogeneous parallel computation is performed to obtain the shape function partial derivative matrix at each Gaussian point; S4 uses the number of the corresponding volume element in the information vector group of the interactive control point pair to determine the shape function partial derivative matrix of the Gaussian point corresponding to the volume element of the interactive control point pair, and uses the shape function partial derivative matrix to calculate the strain matrix, thereby establishing the correspondence between the interactive control point and the strain matrix; the global stiffness coefficient submatrix of the interactive control point pair is calculated using the correspondence, thereby obtaining the global stiffness coefficient submatrix of all interactive control point pairs, and the global stiffness coefficient submatrix of all interactive control point pairs is assembled to realize the heterogeneous parallel assembly of the stiffness matrix.

2. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 1, characterized in that, In step S3, the heterogeneous parallel computation of the partial derivatives of the shape functions is performed as follows: S31, the partial derivatives of the shape functions at all Gaussian points in each volume element are taken as a subtask, thus forming multiple subtasks; S32, the multiple tasks are allocated between the CPU and GPU and CPU and GPU multithreading is started, wherein the number of threads of the CPU is the number of CPU cores, and one thread will be responsible for executing multiple subtasks; the number of threads of the GPU is the number of subtasks allocated to the GPU, and one thread will be responsible for executing one subtask; S33, the partial derivatives of the shape functions are computed in parallel by CPU and GPU multithreading.

3. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 2, characterized in that, In step S31, the two-dimensional shape function is determined according to the following relationship: in, It is a shape function in a two-dimensional parametric coordinate space, where p and q are NURBS orders, i and j are control point labels, ξ and η are coordinate values ​​in the parametric coordinate space, and W i It is a weighting factor with the property of the total tensor product: W i,j =W i W j n and m are the number of parameter coordinates. These are shape functions on coordinates ξ and η in a one-dimensional parameter space, respectively.

4. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 3, characterized in that, In step S33, the partial derivatives of the shape functions are calculated according to the following relationship: in, It is a shape function in a two-dimensional parametric coordinate space, where p and q are NURBS orders, i and j are control point labels, ξ and η are coordinate values ​​in the parametric coordinate space, and W i It is a weighting factor with the property of the total tensor product, W i,j =W i W j n and m are the number of parameter coordinates. These are shape functions on coordinates ξ and η in a one-dimensional parameter space, respectively.

5. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 1, characterized in that, In step S4, the calculation of the global stiffness matrix of the interactive control point pairs is performed according to the following steps: S41 Taking the interactive control point pairs as the unit, the solution of the global stiffness sub-matrix of all interactive control point pairs is divided into multiple sub-tasks to obtain multiple sub-tasks. S42 distributes the subtasks between the CPU and GPU and starts CPU and GPU multi-threading. The CPU thread is responsible for calculating the stiffness coefficient submatrix of multiple interactive control point pairs; the GPU thread is only responsible for calculating the stiffness coefficient submatrix of one interactive control point pair; S43 CPU and GPU multi-threading calculate the global stiffness matrix in parallel.

6. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 5, characterized in that, In step S43, the global stiffness matrix is ​​calculated as follows: S431, the shared units of the interactive control point pairs are determined; S432, the stiffness coefficient sub-unit of each shared unit is calculated; S433, the relationship between the stiffness coefficient sub-unit of each shared unit and the global stiffness coefficient sub-matrix is ​​constructed to obtain the global stiffness coefficient sub-matrix of the interactive control point pairs; S434, the global stiffness coefficient sub-matrixes of all interactive control point pairs are assembled to obtain the global stiffness matrix of the object to be processed.

7. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 6, characterized in that, In step S432, for an interactive control point pair ij, the stiffness coefficient submatrix H of the shared element to which the interactive control point pair belongs... ij Calculate according to the following formula: Where D is the elasticity matrix, B i B j These are the strain matrices corresponding to control points i and j, respectively, and strain matrix B. i and B j The elements in the equation are the partial derivatives of the shape function at the Gaussian point of the volume element, representing the geometric coordinates of the derivatives.

8. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 6, characterized in that, In step S433, the global stiffness coefficient submatrix is ​​determined according to the following relationship: Among them, K ij It is the global stiffness coefficient submatrix, w G It is the weighting coefficient, H ij It is the stiffness coefficient submatrix.

9. The CPU / GPU heterogeneous parallel computing method with equal geometric stiffness matrix as described in claim 6, characterized in that, In step S434, the global stiffness matrix is ​​a sparse matrix, and its elements are indexed in the following manner: edof(i) = [ix, iy] and edof(j) = [jx, jy], where x and y are the coordinate directions of the horizontal and vertical coordinates in the geometric space, respectively, and edof(i) and edof(j) represent the degree-of-freedom label vectors of control points i and j in the x and y directions. This represents the operation of converting the degree-of-freedom labels to row and column indices in the stiffness matrix K, iK ij It is a vector of row index numbers, jK ij It is a vector of column index numbers.

Images