A method, system, medium, and apparatus for numerical solution of finite elements of geometric structures
By employing a solution method based on matrix-free operators and hybrid preconditioners, the problems of high memory consumption and high computational cost in complex geometric structures of the finite element method are solved, achieving efficient mechanical property analysis and response prediction, and applicable to complex working conditions such as multi-material composites and nonlinear deformations.
Patent Information
- Application Number
- CN202610318108.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-03-16
- Publication Date
- 2026-07-10
AI Technical Summary
Existing finite element methods face problems such as high memory consumption, wasted computational resources, and increased computational costs when dealing with complex geometric structures. In particular, high-order finite element methods struggle to effectively handle the fine mechanical characteristics and nonlinear problems of large geometric structures.
The matrix-vector product calculation is performed directly without matrix operators. The solution method is constructed by grid discretization and hybrid preconditioners, avoiding the construction and storage of the global system matrix. The solution is performed using matrix-free operators and hybrid preconditioners, and memory usage and computational cost are reduced by summation, decomposition and combination.
It effectively reduces memory usage, improves solution efficiency, and reduces computational costs. It can efficiently handle the mechanical performance analysis and response prediction of complex geometric structures and is suitable for complex working conditions such as multi-material composites and nonlinear deformation.
Smart Images
Figure CN122365979A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of numerical solution technology for geometric structures, and in particular to a finite element numerical solution method, system, medium, and equipment for geometric structures. Background Technology
[0002] Geometric structures in engineering scenarios often have complex characteristics, such as the curved contours of aerospace components, the splicing of multiple components in building structures, and the irregular cavities of mechanical parts. These structures cannot be described by simple geometric equations corresponding to spheres and cubes. Moreover, geometric structures usually face complex working conditions such as multi-material composites, nonlinear deformation, and dynamic loads. Their mechanical performance analysis and working condition response prediction are the core requirements of engineering design and safety verification.
[0003] The Finite Element Method (FEM), as a numerical method for solving geometrically adaptable problems, effectively handles complex geometries, non-uniform material properties, and nonlinear physical problems by decomposing the solution domain into a finite number of elements and constructing local approximate solutions on these elements. In the finite element analysis of general geometric structures, the solution process typically includes two stages: system matrix construction and linear system solution. First, based on the mathematical description of the mechanical properties, boundary constraints, and material properties of the geometric structure, a discretized global system matrix A and a right-hand side vector f are constructed to obtain a linear system Au=f representing the force-displacement relationship of the geometric structure. Then, this linear system is solved to obtain the finite element solution u.
[0004] Traditional finite element analysis methods for general geometric structures employ explicit matrix methods, which involve first explicitly constructing and storing the global system matrix A, and then solving it using numerical linear algebra methods. However, this method suffers from the following problems: 1) As the demand for geometric structure analysis increases, such as using high-order finite element methods to obtain the fine mechanical characteristics of complex geometric structures, or performing ultra-fine meshing on large geometric structures, the dimension of the global system matrix A will increase significantly, which may exceed the carrying capacity of computer memory. Even if a sparse matrix storage format is used, memory consumption may still become a bottleneck. 2) In the iterative solution of geometric linear systems, sparse matrix-vector multiplication is the core computational operation. However, due to the irregularity of sparse matrix data access patterns, its theoretical computational intensity is far lower than the ideal computational intensity that modern hardware can provide, resulting in a waste of computational resources and significantly limiting the efficiency of solving geometric structures. 3) When dealing with nonlinear or time-dependent problems, such as geometric structures with large deformation or plastic deformation of materials, the global system matrix A may need to change with the deformation and stress state of the geometric structure and be reassembled in each iteration or time step, which significantly increases the computational cost. This is especially true in high-order finite element methods, where the complexity of element matrix generation and global assembly is high.
[0005] Therefore, for geometric structures, there is an urgent need for a finite element numerical solution method that can reduce memory usage, improve solution efficiency, and reduce computational costs. Summary of the Invention
[0006] In order to overcome the shortcomings of the prior art, the present invention aims to provide a finite element numerical solution method, system, medium and device for geometric structures, which can reduce memory usage, improve solution efficiency and reduce computational cost.
[0007] This invention is implemented according to the following scheme: A finite element numerical solution method for geometric structures is provided, including: Obtain the structural parameters and boundary constraint parameters of the geometric structure, wherein the boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under a displacement field or a force field; Based on the structural parameters and the boundary constraint parameters, construct the geometric model of the geometric structure; The geometric model is discretized to obtain the discretization parameters of the geometric model and the geometric mesh sequence constituting the geometric model. The geometric mesh sequence includes geometric meshes of multiple levels. Based on the discretization parameters, a matrix-free operator for the geometric mesh is constructed based on summation decomposition. The matrix-free operator is used to directly perform matrix-vector product calculations in the geometric model solution process. Based on the matrix-free operator and the geometric mesh sequence, a hybrid preconditioner is constructed; The geometric model is solved based on the hybrid preconditioner to obtain the solution result of the geometric structure. The solution result includes the structural data of the geometric structure under the constraints of displacement field or force field.
[0008] Optionally, the geometric model is discretized into a mesh to obtain the geometric mesh sequence, including: The solution domain of the geometric model is partitioned using hexahedral elements to obtain an initial discrete mesh, wherein the hexahedral elements are the basic geometric elements constituting the initial discrete mesh; The initial discrete mesh is refined multiple times to generate multiple geometric meshes ordered by their coarseness, where the coarseness is determined by the order of the refinement process.
[0009] Optionally, the discretization parameters include: The subdivision density parameter of the hexahedral element, which is used to indicate the mesh accuracy of the geometric mesh sequence at different levels; p-order tensor product Lagrange basis function, used to discretize the weak form equations of the geometric model on the hexahedral elements; The geometric mapping Jacobian matrix is used to convert the physical domain of the geometric model into the reference domain of a preset reference cell.
[0010] Optionally, based on the discretization parameters, constructing a matrix-free operator for the geometric mesh using a summation decomposition includes: Based on the discretization parameters, the three-dimensional summation process of the geometric mesh at the integration points is decomposed into three-dimensional tensor contraction calculations in three directions. The one-dimensional tensor contraction calculations in the three directions are performed respectively to construct the matrix-free operator, including: By forward transformation, the displacement gradient at the unit node of the geometric mesh is obtained, which maps the displacement at the integration point to the displacement gradient at the integration point. Based on the displacement gradient, a preset geometry and material operator is applied at the integration point to perform a physical field evaluation, and the stress of the geometric mesh at the integration point is obtained. The stress is mapped back to the element node by a backward transformation to obtain the matrix-free operator.
[0011] Optionally, before constructing the matrix-free operator by calculating the one-dimensional tensor contraction in the three directions respectively, the solution method further includes: The cell node data of the geometric mesh is rearranged from the node-first format to the component-first format. The node-first format is used to indicate that multiple physical components of the same cell node are stored consecutively, and the component-first format is used to indicate that the same physical component of different cell nodes is stored consecutively. Based on forward transformation, physical field evaluation, and backward transformation, a target computation kernel is constructed, which is used to perform one-dimensional tensor contraction calculations in three directions; Based on the one-dimensional tensor contraction in three directions, the computational order used to indicate the calculation of the one-dimensional tensor contraction in the three directions is determined.
[0012] Optionally, the hybrid preconditioner is constructed based on the matrix-free operator and the geometric mesh sequence, including: Determine the diagonal inverse matrix based on the matrix-free operator; Based on the matrix-free operator and the diagonal inverse matrix, a smoother is constructed, which is used to solve the remaining geometric grids in the geometric grid sequence except for the coarsest geometric grid; Based on the algebraic multigrid algorithm, a coarse grid solver is constructed to solve the coarsest geometric grid. Based on the multigrid method, an inter-layer transport operator is constructed, which is used to transmit data between geometric grids of adjacent layers. The hybrid preconditioner is constructed based on the smoother, the coarse mesh solver, and the interlayer transport operator.
[0013] Optionally, the geometric model is solved based on the hybrid preconditioner to obtain the solution result of the geometric structure, including: The hybrid preconditioner is embedded into the conjugate gradient method to construct a solver; The solver iteratively processes the algebraic equations of the geometric model to obtain the solution result.
[0014] A finite element numerical solution system for geometric structures is also provided, which applies the aforementioned finite element numerical solution method for geometric structures, including: The model building module is used for: Obtain the structural parameters and boundary constraint parameters of the geometric structure, wherein the boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under a displacement field or a force field; Based on the structural parameters and the boundary constraint parameters, construct the geometric model of the geometric structure; The processing module is used to perform mesh discretization processing on the geometric model to obtain the discretization parameters of the geometric model and the geometric mesh sequence constituting the geometric model, wherein the geometric mesh sequence includes multiple levels of geometric mesh; Based on the discretization parameters, a matrix-free operator for the geometric mesh is constructed using summation decomposition. Based on the matrix-free operator and the geometric mesh sequence, a hybrid preconditioner is constructed; The solver module is used to solve the geometric model according to the hybrid preconditioner to obtain the solution result of the geometric structure. The solution result includes the structural data of the geometric structure under the constraints of displacement field or force field.
[0015] A computer-readable storage medium is also provided, which is a computer-readable storage medium storing a computer program thereon, which, when executed, implements the finite element numerical solution method for a geometric structure.
[0016] A computer device is also provided, including a processor and a memory, wherein the memory stores at least one instruction, at least one program, code set or instruction set, and the at least one instruction, at least one program, code set or instruction set is loaded and executed by the processor to implement the finite element numerical solution method for a geometric structure.
[0017] Compared with the prior art, the beneficial effects of the finite element data solution method for geometric structures of the present invention are as follows: by directly executing the matrix-vector product calculation in the solution process of the geometric model without matrix operators, the matrix-vector product calculation can be performed without the construction and storage process of the global system matrix A in the traditional finite element numerical solution process. It does not require memory to store the global system matrix A, nor does it require reassembling the global system matrix A in each iteration, which effectively reduces memory usage, improves solution efficiency and reduces computational cost. Attached Figure Description
[0018] Figure 1 This is a flowchart of the solution method of the present invention; Figure 2 This is a schematic diagram illustrating the discretization process of the geometric model using the solution method of this invention; Figure 3 This is a schematic diagram of a traditional assembly method and a linear system; Figure 4 This is a schematic diagram of the linear system corresponding to the solution method and geometric model of this invention; Figure 5 This is a schematic diagram illustrating the data rearrangement and operator fusion performed by the solution method of the present invention before constructing the matrix-free operator. Detailed Implementation
[0019] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.
[0020] It should be noted that, for the sake of simplicity, the method embodiments are all described as a series of actions. However, those skilled in the art should understand that the embodiments of this application are not limited to the described order of actions, because according to the embodiments of this application, some steps can be performed in other orders or simultaneously. Secondly, those skilled in the art should also understand that the embodiments described in the specification are all preferred embodiments, and the actions involved are not necessarily required by the embodiments of this application.
[0021] In the following description, when referring to the accompanying drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this application as detailed in the appended claims. In the description of this application, it should be understood that the terms "first," "second," "third," etc., are used only to distinguish similar objects and are not necessarily used to describe a specific order or sequence, nor should they be construed as indicating or implying relative importance. Those skilled in the art can understand the specific meaning of the above terms in this application according to the specific circumstances.
[0022] The solution method of this invention is aimed at structures with complex geometric features in engineering, such as the curved contours of aerospace components, multi-component splices of building structures, and irregular cavities of mechanical parts. Under given displacement field constraints (Dirichlet boundary conditions) or force field constraints (Neumann boundary conditions), the displacement field distribution u(x) of the geometric structure is obtained through the finite element numerical method. Then, the key mechanical response data such as strain field ε(x) and stress field σ(x) are derived from the displacement field, thereby obtaining complete mechanical characteristic information of the geometric structure under actual working conditions.
[0023] See Figure 1 As shown, the present invention provides a finite element numerical solution method for geometric structures, comprising: S1: Obtain the structural parameters and boundary constraint parameters of the geometric structure. The boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under displacement or force fields. In one embodiment of the present invention, the structural parameters include geometric feature parameters, material property parameters, and body force parameters. The geometric feature parameters represent the three-dimensional dimensional parameters of the geometric structure and the topological parameters of the solution domain Ω, where the solution domain Ω is the physical solution domain range of the geometric structure. The material property parameters represent the Lamé constant λ and shear modulus μ of the isotropic linear elastic material, and are used to construct the constitutive relationship between stress and strain. These parameters change with spatial position. The body force parameters are body force vectors f acting inside the geometric structure, and are used to represent the stress state inside the geometric structure.
[0024] In one embodiment of the present invention, the boundary constraint parameters include displacement constraint parameters and force constraint parameters. The displacement constraint parameters are Dirichlet boundary conditions, which are used to represent the displacement limit value of the fixed end or specified boundary of the geometric structure, such as limiting the displacement u=0 of the fixed end of the geometric structure to define the displacement boundary state of the geometric structure. The force constraint parameters are Neumann boundary conditions, which are used to represent the surface force vector of the force-bearing end or free boundary of the geometric structure to represent the force situation of the geometric structure boundary, such as the geometric structure being subjected to a constant tensile force f.
[0025] S2: Based on the structural parameters and boundary constraint parameters, construct the geometric model of the geometric structure, including: determining the three-dimensional physical solution domain occupied by the geometric structure based on the geometric feature parameters in the structural parameters. To determine the overall topology, spatial dimensions, and subdomain partitioning characteristics of the solution domain Ω, a basic geometric outline of the geometric structure is formed. If the geometric structure is a multi-material structure, the solution domain is divided according to the materials to obtain sub-solution domains corresponding to different materials. Next, the static equilibrium of the linear elastic problem in structural mechanics is constructed into a strong-form governing equation, including: under the assumption of small deformation, the static equilibrium equation satisfied by the displacement field u is constructed, and its expression is as follows: σ+f=0 Where σ is the Cauchy stress tensor, and f is the body force vector in the structural parameters. σ is the divergence of the stress tensor, used to represent the variation and equilibrium relationship of stress in space; Based on the material property parameters in the structural parameters, a linear elastic constitutive model is constructed, including: using an isotropic linear elastic constitutive model, establishing the relationship between stress and strain, the expression of which is as follows: σ=λ( u)I+2με(u) Where λ and μ are Lamé constants, ε(u) = ( u+( u) T ) / 2 is the linear strain tensor. The geometric structures of multiple materials are matched with the corresponding Lamé constants according to the sub-solution domains. u is the divergence of the displacement field, and its value is equal to the trace of the linear strain tensor. , used to characterize the volumetric deformation of the geometric structure, u is the displacement field vector, I is the unit tensor, ( u) T It is the transpose of the displacement gradient tensor; Based on the variational principle, the strong-form governing equations are transformed into weak-form equations requiring solutions. This involves multiplying both sides of the static equilibrium equations by the test function and integrating by parts, then eliminating higher-order derivatives to reduce the regularity requirement of the solution function, resulting in the weak-form equations, expressed as follows:
[0026] Where Ω represents the three-dimensional physical solution domain occupied by the geometric structure, and v is the test function used to perform weighted integration on the strong-form governing equations to derive the weak-form equations. For the strain tensor of the test function, The force constraint parameter is the Neumann boundary condition at the force-bearing end of the geometric structure. This weak form equation can be adapted to the finite element solution of complex geometric structures and boundary conditions. At the same time, this weak form equation reduces the requirement for the continuity of the solution and is suitable for finite element discrete solution. Finally, based on the boundary constraint parameters, the boundaries of the geometric model are determined, including marking the boundary properties such as displacement constraints (Dirichlet boundaries, such as fixed end u=0), force constraints (Neumann boundaries, such as force g on the stressed end face) and free surfaces at the corresponding boundary positions of the geometric model, so as to embed the boundary constraint parameters into the boundary integral terms of the weak form equations and clarify the numerical constraint conditions of each boundary. Based on the basic geometric profile, linear elastic constitutive model, weak formal equations and boundaries obtained above, a complete geometric model is constructed. This geometric model can simultaneously represent the geometric shape, material distribution, physical and mechanical characteristics and full boundary constraints of the geometric structure.
[0027] S3: Perform mesh discretization on the geometric model to obtain the discretization parameters of the geometric model and the geometric mesh sequence that constitutes the geometric model. The geometric mesh sequence includes geometric meshes of multiple levels. In one embodiment of the present invention, the geometric model is discretized to obtain a geometric mesh sequence, including: using hexahedral elements to divide the solution domain of the geometric model to obtain an initial discrete mesh, wherein the hexahedral elements are the basic geometric elements constituting the initial discrete mesh; and performing multiple refinement processes on the initial discrete mesh to generate multiple geometric meshes sorted by coarseness, wherein the coarseness is the order of the refinement processes.
[0028] In one embodiment of the present invention, hexahedral elements are used to partition the solution domain of the geometric model. Specifically, the continuous three-dimensional physical solution domain Ω corresponding to the geometric model is divided into subdivisions. R 3 A set T divided into a series of hexahedral physical units hThe initial discrete mesh exists as the coarsest layer (Level 0) in the geometric mesh sequence. The multiple refinement processes performed on the initial discrete mesh are all uniform refinement processes. Specifically, after the initial discrete mesh is uniformly refined once, mesh Level 1 is obtained. After mesh Level 1 is uniformly refined once, mesh Level 2 is obtained, and so on. After L uniform refinement processes, fine-layer meshes from mesh Level 1 to mesh Level L are generated in the refinement order to obtain a complete geometric mesh sequence. Mesh Level L, which exists as the finest layer, is the highest level of this geometric mesh sequence. Moreover, the geometric meshes at each level maintain the tensor product structure of hexahedral elements, which is precisely matched with the topology of the geometric model, the multi-material subdomain partitioning, and the boundary constraint positions.
[0029] In one embodiment of the present invention, each level of the geometric mesh sequence has the core features of the geometric model. All hexahedral elements are labeled with material properties according to the material distribution of the geometric model. The hexahedral elements of the multi-material geometric structure are matched with the corresponding Lamé constants according to the sub-solution domain corresponding to the material. At the same time, the hexahedral elements at the boundary are precisely matched with the displacement constraint parameters (Dirichlet boundary conditions) and force constraint parameters (Neumann boundary conditions) of the geometric model to completely preserve the boundary constraint markings of the geometric model. This ensures that the boundary constraint conditions after discretization of the geometric model can be effectively transferred to the geometric meshes at each level. The geometric meshes at each level together constitute the discretized solution domain of the geometric model, providing a discretization carrier for subsequent solutions.
[0030] In practical applications, the solution domain of a geometric model can be divided using elements such as triangles, quadrilaterals, and tetrahedrons. Finite-dimensional basis functions are defined on these dividing elements to discretize the continuous problem, enabling the geometric model to be solved by a computer. In one embodiment of the present invention, the solution domain of the geometric model is divided using hexahedral elements to obtain the fewest possible sub-elements of the solution domain. If the geometric structure is made of multiple materials, hexahedral elements are used to divide the sub-solution domains corresponding to different materials.
[0031] In one embodiment of the present invention, the discretization parameters include: The mesh density parameter of the hexahedral element is used to indicate the mesh accuracy of the geometric mesh sequence at different levels. This mesh density parameter is positively correlated with the size and number of hexahedral elements in each level of the geometric mesh. The mesh density parameter increases progressively with the refinement level (level) of the geometric mesh. It is used not only to indicate the accuracy of each level of the geometric mesh, but also to indicate the approximation accuracy of each level of the geometric mesh to the geometric shape and physical characteristics of the geometric model. Thus, the mesh density parameter at different levels can be used to adapt to the finite element numerical solution requirements of different scales and accuracy requirements.
[0032] The p-order tensor product Lagrange basis function is used to discretize the weak formal equations of the geometric model on hexahedral elements. This p-order tensor product Lagrange basis function is obtained based on a predefined reference element and is a three-dimensional basis function. It can be expressed as the product of three one-dimensional basis functions, as follows:
[0033] Where p is the order of the p-th order tensor product Lagrange basis function, the weak form of the geometric model is embedded into the hexahedral element through this three-dimensional basis function, and the discrete form of the weak form equation can be obtained, which is expressed as follows:
[0034] in, As elements of the global system matrix A, this invention discretizes weak formal equations by embedding them into hexahedral elements, thereby transforming continuous mechanics problems into discrete linear systems.
[0035] The geometric mapping Jacobian matrix is used to transform the physical domain of the geometric model into the reference domain of a predefined reference element; after partitioning the solution domain Ω, a set of hexahedrons T is obtained. h The preset reference unit is =[ 1,1] 3 For any predefined reference unit, K∈T h There exists a geometric mapping ΦK from a preset reference element to a hexahedral element: →K; the Jacobian matrix of the geometric mapping is... This parameter, in addition to realizing the transformation from the physical domain to the reference domain, can also perform the conversion from the physical domain gradient to the reference domain gradient, and the conversion process satisfies... .
[0036] After transforming the physical domain of the geometric model into the reference domain of the preset reference element using the aforementioned geometric mapping Jacobian matrix, the expression for the local matrix-vector product is as follows:
[0037] in, Here is the element stiffness matrix corresponding to the hexahedral element K; Let be the displacement vector to be determined; subscript These are the node numbers of the hexahedral element in the three reference coordinate directions; These are the Gaussian integration points numbered in the three reference coordinate directions; These are the integration weights corresponding to the Gaussian integration points in the three directions, respectively. For geometric mapping Jacobian matrix The absolute value of the determinant; It is the transpose of the inverse of the Jacobian matrix of the geometric mapping; The gradient operator in the reference coordinate system; The p-th order tensor product Lagrange basis function on the preset reference unit; This is the fourth-order material stiffness tensor (constitutive tensor) for isotropic linear elastic materials. Determined by Lamé constant λ and shear modulus μ; This is the finite element approximation function of the displacement field; The coordinates of the Gaussian integration point in the reference domain; the colon ":" indicates the double tensor shrinking operation.
[0038] If the matrix-vector multiplication described above is performed directly, the complexity of the three-dimensional summation is O(p). 6 ).
[0039] See Figure 2 The diagram illustrates the process of discretizing the geometric model. The structure diagram on the left represents the physical domain of the geometric structure, which is also the solution domain. This solution domain is decomposed into multiple hexahedral physical elements through subdomain partitioning. The diagram on the right shows the preset reference elements. =[ 1,1] 3 (Reference Element) Each hexahedral physical element is associated with a reference element through geometric mapping.
[0040] S4: Based on the discretization parameters, construct a matrix-free operator for the geometric mesh based on summation decomposition. The matrix-free operator is used to directly perform matrix-vector product calculations in the geometric model solution process. In one embodiment of the present invention, constructing a matrix-free operator for a geometric mesh based on summation decomposition according to discretization parameters includes: decomposing the three-dimensional summation process of the geometric mesh at the integration points into one-dimensional tensor contraction calculations in three directions according to the discretization parameters. The present invention, through summation decomposition, reduces the computational complexity of matrix-vector multiplication from O(p...)...6 ) decreased to O(p 4 Perform one-dimensional tensor contraction calculations in three directions respectively to construct a matrix-free operator, including: By performing a forward transformation, we obtain the displacement gradient that maps the displacement of the geometric mesh at the element nodes to the displacement at the integration point, as expressed below:
[0041] in, At the integration point after forward transformation The displacement gradient calculated at that location; These are the one-dimensional p-order Lagrange basis functions in the three reference coordinate directions, respectively; For a one-dimensional p-th order Lagrange basis function At the points The derivative value at that point; For unit nodes Displacement value at; Let be the coordinates of the Gaussian integral points in the three directions of the reference domain; the formula decomposes the original three-dimensional tensor summation into a one-dimensional tensor contraction calculation performed sequentially along the first, second, and third reference coordinate directions through three nested summations.
[0042] Based on the displacement gradient, a physical field evaluation is performed at the integration point using pre-defined geometry and material operators, yielding the stress of the geometric mesh at the integration point, expressed as follows:
[0043] in, Including pre-calculated geometric factors , As a material stiffness tensor, this physical field evaluation process is entirely point-wise, making it easy to vectorize.
[0044] By mapping the stress back to the element nodes through a backward transformation, a matrix-free operator is obtained, the expression of which is as follows:
[0045] in, Mapping back to element nodes after backward transformation Equivalent nodal force at the location; For the integration point after evaluation by the physics field The calculated stress value is obtained at the point. This formula is the transpose of the forward transformation, which projects the stress at the integration point back to the element node through three nested summations, completing the element-level matrix-vector product calculation without matrix operators.
[0046] In one embodiment of the present invention, before calculating the one-dimensional tensor contraction in three directions and constructing the matrix-free operator, the solution method further includes: The element node data of the geometric mesh is rearranged from the node-first format AoS to the component-first format SoA. The node-first format is used to indicate that multiple physical components of the same element node are stored contiguously. If a node-first format is used for storage, vectorized instructions (such as AVX-512 or GPU Warp) will experience non-contiguous memory accesses when reading data, greatly reducing bandwidth utilization. Component-first format is used to indicate that the same physical component in different cell nodes is stored contiguously, i.e. By using a component-first storage format, the same VectorRegister can process the physical components of multiple integration points in the same direction simultaneously, thereby fully utilizing the Single Instruction Multiple Data Stream (SIMD) parallel capability of modern CPUs / GPUs to achieve batch computation. Based on forward transformation, physics evaluation, and backward transformation, a target computation kernel is constructed to perform one-dimensional tensor contraction calculations in three directions. Specifically, this kernel performs chained calculations from shape function derivatives to stress tensors directly at the cache level, without writing intermediate variables back to main memory. By reducing the number of read / write operations to main memory, the algorithm's arithmetic intensity is significantly improved, transforming it from memory-bound to computation-bound, thereby greatly increasing the solution speed and effectively avoiding the need to calculate stress in traditional matrix-free implementations. The process is divided into three independent steps: calculating strain, reading material parameters, and calculating stress. This leads to a large amount of intermediate data being frequently written to and read from memory (Memory Traffic), which is limited by the memory bandwidth bottleneck. Based on the one-dimensional tensor contraction in three directions, the calculation order for indicating the one-dimensional tensor contraction calculation in the three directions is determined. The 6×6 (3D) matrix operation involving material coupling is embedded in the innermost loop, thereby maximizing the cache hit rate while ensuring the calculation accuracy.
[0047] See Figure 3 The diagram illustrates the traditional assembly method and linear system, showcasing the core logic and data flow process of the traditional finite element method. As can be seen from the diagram, the traditional method requires explicit assembly of the global system matrix A, with matrix elements derived from... = The calculation yields a linear system Au=f for the geometric model, which is then constructed based on the global system matrix A. The figure clearly shows the intermediate process of matrix assembly, including element matrix generation, global index mapping, and sparse matrix storage. This highlights the dependence of traditional methods on global matrix storage. Even with a sparse storage format, memory usage will still increase dramatically as the scale of problems such as high-order elements or fine meshes expands. See Figure 4 The figure shows the linear system corresponding to the solution method and geometric model of this invention. As can be seen from the figure, this invention does not use an explicit global system matrix A in the solution process. Instead, it uses a matrix-free operator as the core, directly connecting the input displacement field u and the output nodal force vector f. The specific process is as follows: after the input displacement field u is rearranged (AoS rearranged to SoA), it is decomposed into one-dimensional tensor contraction calculations in three directions using a summation decomposition technique. The forward transformation (displacement gradient interpolation), physical field evaluation (stress calculation), and backward transformation (nodal force projection) are completed sequentially, and finally the nodal force vector f is output, realizing the equivalent calculation of the linear system f=Au. This invention skips the matrix assembly step by using a matrix-free operator and directly calculates the matrix-vector product f=Au. The entire process only relies on discretized parameters and geometric mesh data, compared to... Figure 3 The present invention does not require the storage of large-scale matrices, unlike the conventional assembly method shown.
[0048] See Figure 5 As shown, the core logic of data rearrangement and operator fusion before constructing a matrix-free operator is illustrated. The upper part of the figure shows the traditional matrix-free implementation process, where data is stored in node-first AoS format. The calculation process is divided into independent steps such as displacement gradient interpolation, strain calculation, stress solution, and nodal force projection, involving a large amount of intermediate data memory read and write, which is limited by memory bandwidth. The lower part of the figure shows the solution method of this invention. First, the data is rearranged into component-first SoA format to adapt to the continuous memory access requirements of the SIMD instruction set. Then, multiple independent steps are integrated into a single computational kernel through operator fusion, completing the chained calculation from displacement to nodal force at the cache level. Intermediate variables are not written back to main memory, and there is no need for frequent read and write operations on intermediate data, effectively reducing computational costs and improving computational efficiency.
[0049] S5: Based on the matrix-free operator and the geometric mesh sequence, construct a hybrid preconditioner, including: Determine the diagonal inverse matrix based on the matrix-free operator, specifically by using elements of the global system matrix. By defining the displacement field u, we replace it with the square form of the basis functions, i.e., we can calculate it directly:
[0050] in, For the global stiffness matrix and the element nodes The corresponding diagonal element values; GeometricFactors are pre-calculated geometric factors, including Gaussian integral weights. The absolute value of the determinant of the geometric mapping Jacobian matrix and the transpose of the inverse of the geometric mapping Jacobian matrix. The relevant combination quantity; These represent the squares of the one-dimensional basis functions at the corresponding integration points in the three directions; the ellipsis indicates additional product terms related to the squared derivatives of the basis functions and the material stiffness tensor; this formula obtains the diagonal vector of the global stiffness matrix in O(N) time complexity through the same tensor contraction path as the matrix-free operator, without the need to assemble the global matrix.
[0051] By using the same tensor shrinkage path, the diagonal vector of the global system matrix can be obtained accurately in O(N) time, without the need for the matrix assembly step in the traditional solution process, thus reducing computational cost.
[0052] Based on the matrix-free operator and the diagonal inverse matrix, a smoother is constructed. This smoother is used to solve for the remaining geometric meshes in the geometric mesh sequence, excluding the coarsest one. Specifically, the smoother employs a matrix-free Chebyshev iteration. This level consumes the most memory, but the matrix-free method can significantly reduce memory overhead. The iterative formula for this smoother is as follows:
[0053] in, For the first Approximate solution of displacement after the next iteration; For the first Approximate solution of displacement after the next iteration; It is the diagonal inverse matrix; The Chebyshev coefficients are based on the operator's largest eigenvalue λ. max The estimated value is determined, and the largest eigenvalue λ max The result was quickly estimated through several applications of the Power Method on matrix-free operators; This is the right-hand term vector (load vector) of the current grid level. For the first Matrix-free operators corresponding to hierarchical geometric meshes; For the first The residual vector of the next iteration.
[0054] Based on the algebraic multigrid algorithm (AMG), a coarse-grid solver is constructed to solve the coarsest geometric mesh. In the geometric mesh sequence, the coarsest layer is the Level 0 (initial discrete mesh), which has a relatively small number of degrees of freedom. The solution efficiency of this layer significantly affects the overall solution efficiency of the linear algebraic system Au=f corresponding to the geometric model. Therefore, the algebraic multigrid algorithm is used to solve the equations of the coarsest mesh. This can achieve the highest computational performance, among which This is the error correction vector, used to represent the deviation between the current approximate solution and the exact solution. After the solution is completed, this correction vector is superimposed on the current approximate solution to improve the solution accuracy. This is the residual vector, used to quantify the degree to which the current approximate solution satisfies the linear system. This is the sparse stiffness matrix corresponding to the coarsest level of the geometric mesh.
[0055] Based on the V-Cycle multigrid method, an inter-layer transport operator is constructed to transmit data between geometric meshes of adjacent layers. The inter-layer transport operator includes a tensor product-based constraint operator and a tensor product-based interpolation operator, both constructed based on the tensor product structure of hexahedral elements of the geometric mesh, adapting to the discretization characteristics of each layer's geometric mesh. Pre-smoothing is performed at the finest layer (L-level) to calculate the residuals. At this point, the constraint operator in the form of tensor product is used. Project the residuals to the next level (L-1 level): until recursion to L The algebraic multigrid algorithm (AMG) is called to solve the problem, and then an interpolation operator in the form of tensor product is used. Perform correction: This is used to perform post-smoothing on the coarsest level (level 0).
[0056] A hybrid preconditioner is constructed based on a smoother, a coarse mesh solver, and an interlayer transport operator.
[0057] S6: Solve the geometric model based on the hybrid preconditioner to obtain the solution results of the geometric structure. The solution results include the structural data of the geometric structure under the constraints of displacement field or force field.
[0058] In one embodiment of the present invention, the geometric model is solved according to the hybrid preconditioner to obtain the solution result of the geometric structure, including: embedding the hybrid preconditioner into the conjugate gradient method to construct a solver; the solver iteratively processes the algebraic equations of the geometric model to obtain the solution result.
[0059] This invention transforms the original linear system into an equivalent preconditioning system by embedding a hybrid preconditioner into the conjugate gradient method, thereby reducing the condition number and significantly reducing the number of iterations required to achieve a given convergence accuracy, thus accelerating the iteration speed during the solution process.
[0060] The solution results obtained by the solution method of this invention can provide quantitative structural data such as displacement, stress, and strain for the mechanical performance analysis of geometric structures, supporting strength verification and stiffness assessment in engineering design. It can also provide high-precision numerical solutions for predicting the response of geometric structures under complex working conditions such as multi-material composites, nonlinear deformation, and dynamic loads, assisting in engineering safety verification and structural optimization design. Furthermore, through the collaborative solution strategy of matrix-free operators and hybrid preconditioners, the solution of large-scale, high-order finite element problems can be achieved under acceptable computational resource constraints, thereby expanding the applicability of finite element numerical methods in the engineering analysis of large and complex geometric structures.
[0061] In one embodiment of the present invention, the solver's solution framework is a general-purpose solution framework that supports three assembly levels. The user or software can automatically select the following three levels based on hardware resources: Level 1 (Fully Assembled): Suitable for small-scale problems or large memory nodes, all levels assemble the matrix, and the calculation speed is the fastest. Level 2 (Partial Assembly): Suitable for medium-sized problems, assembling only diagonal pieces or specific levels; Level 3 (no matrix assembly): Applicable to large-scale high-order problems, assembly is performed only in the coarsest mesh layer, and the remaining levels have no matrix assembly process.
[0062] This invention directly performs matrix-vector product calculations in the solution process of the geometric model without matrix operators. It can perform matrix-vector product calculations without the construction and storage process of the global system matrix A in the traditional finite element numerical solution process. It does not require memory to store the global system matrix A, nor does it require reassembling the global system matrix A in each iteration, which effectively reduces memory usage, improves solution efficiency, and reduces computational costs.
[0063] The present invention provides a finite element numerical solution system for a geometric structure, which applies the above-described finite element numerical solution method for a geometric structure, including: The model building module is used for: Obtain the structural parameters and boundary constraint parameters of the geometric structure. The boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under displacement or force fields. Based on the structural parameters and boundary constraint parameters, construct the geometric model of the geometric structure; The processing module is used to perform mesh discretization on the geometric model to obtain the discretization parameters of the geometric model and the geometric mesh sequence that constitutes the geometric model. The geometric mesh sequence includes geometric meshes of multiple levels. A matrix-free operator for constructing a geometric mesh based on the discretization parameters and summation decomposition is used. Hybrid preconditioners are constructed based on matrix-free operators and geometric mesh sequences; The solver module is used to solve the geometric model based on the hybrid preconditioner to obtain the solution results of the geometric structure. The solution results include the structural data of the geometric structure under the constraints of displacement field or force field.
[0064] This invention discloses a computer-readable storage medium storing a computer program thereon, which, when executed, implements the aforementioned solution method.
[0065] Optionally, the computer-readable storage medium may include: read-only memory (ROM), random access memory (RAM), solid-state drives (SSDs), or optical discs, etc. The random access memory may include resistive random access memory (ReRAM) and dynamic random access memory (DRAM).
[0066] The computer device of the present invention includes a processor and a memory. The memory stores at least one instruction, at least one program, code set, or instruction set, and the at least one instruction, at least one program, code set, or instruction set is loaded and executed by the processor to implement the above-described solution method.
[0067] The processor can be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc.
[0068] The memory can be used to store the computer program or module. The processor implements various functions of the solution method by running or executing the computer program or module stored in the memory and calling the data stored in the memory. The memory may mainly include a program storage area and a data storage area. The program storage area may store the operating system, at least one application program required for a function, etc.; the data storage area may store data created based on the use of the mobile phone, etc. In addition, the memory may include high-speed random access memory, and may also include non-volatile memory, such as hard disk, RAM, plug-in hard disk, smart media card (SMC), secure digital (SD) card, flash card, at least one disk storage device, flash memory device, or other volatile solid-state storage device.
[0069] The above are merely preferred embodiments of this application and are not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.
Claims
1. A finite element numerical solution method for geometric structures, characterized in that, include: Obtain the structural parameters and boundary constraint parameters of the geometric structure, wherein the boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under a displacement field or a force field; Based on the structural parameters and the boundary constraint parameters, construct the geometric model of the geometric structure; The geometric model is discretized to obtain the discretization parameters of the geometric model and the geometric mesh sequence constituting the geometric model. The geometric mesh sequence includes geometric meshes of multiple levels. Based on the discretization parameters, a matrix-free operator for the geometric mesh is constructed based on summation decomposition. The matrix-free operator is used to directly perform matrix-vector product calculations in the geometric model solution process. Based on the matrix-free operator and the geometric mesh sequence, a hybrid preconditioner is constructed; The geometric model is solved based on the hybrid preconditioner to obtain the solution result of the geometric structure. The solution result includes the structural data of the geometric structure under the constraints of displacement field or force field.
2. The finite element numerical solution method for a geometric structure according to claim 1, characterized in that, The geometric model is discretized into a mesh to obtain the geometric mesh sequence, including: The solution domain of the geometric model is partitioned using hexahedral elements to obtain an initial discrete mesh, wherein the hexahedral elements are the basic geometric elements constituting the initial discrete mesh; The initial discrete mesh is refined multiple times to generate multiple geometric meshes ordered by their coarseness, where the coarseness is determined by the order of the refinement process.
3. The finite element numerical solution method for a geometric structure according to claim 2, characterized in that, The discretization parameters include: The subdivision density parameter of the hexahedral element, which is used to indicate the mesh accuracy of the geometric mesh sequence at different levels; p-order tensor product Lagrange basis function, used to discretize the weak form equations of the geometric model on the hexahedral elements; The geometric mapping Jacobian matrix is used to convert the physical domain of the geometric model into the reference domain of a preset reference cell.
4. The finite element numerical solution method for a geometric structure according to claim 1, characterized in that, Based on the discretization parameters, a matrix-free operator for constructing the geometric mesh is generated using a summation decomposition, including: Based on the discretization parameters, the three-dimensional summation process of the geometric mesh at the integration points is decomposed into three-dimensional tensor contraction calculations in three directions. The one-dimensional tensor contraction calculations in the three directions are performed respectively to construct the matrix-free operator, including: By forward transformation, the displacement gradient at the unit node of the geometric mesh is obtained, which maps the displacement at the integration point to the displacement gradient at the integration point. Based on the displacement gradient, a preset geometry and material operator is applied at the integration point to perform a physical field evaluation, and the stress of the geometric mesh at the integration point is obtained. The stress is mapped back to the element node by a backward transformation to obtain the matrix-free operator.
5. The finite element numerical solution method for a geometric structure according to claim 4, characterized in that, Before constructing the matrix-free operator by calculating the one-dimensional tensor contraction in the three directions respectively, the solution method further includes: The cell node data of the geometric mesh is rearranged from the node-first format to the component-first format. The node-first format is used to indicate that multiple physical components of the same cell node are stored consecutively, and the component-first format is used to indicate that the same physical component of different cell nodes is stored consecutively. Based on forward transformation, physical field evaluation, and backward transformation, a target computation kernel is constructed, which is used to perform one-dimensional tensor contraction calculations in three directions; Based on the one-dimensional tensor contraction in three directions, the computational order used to indicate the calculation of the one-dimensional tensor contraction in the three directions is determined.
6. The finite element numerical solution method for a geometric structure according to claim 1, characterized in that, The hybrid preconditioner is constructed based on the matrix-free operator and the geometric mesh sequence, including: Determine the diagonal inverse matrix based on the matrix-free operator; Based on the matrix-free operator and the diagonal inverse matrix, a smoother is constructed, which is used to solve the remaining geometric grids in the geometric grid sequence except for the coarsest geometric grid; Based on the algebraic multigrid algorithm, a coarse grid solver is constructed to solve the coarsest geometric grid. Based on the multigrid method, an inter-layer transport operator is constructed, which is used to transmit data between geometric grids of adjacent layers. The hybrid preconditioner is constructed based on the smoother, the coarse mesh solver, and the interlayer transport operator.
7. The finite element numerical solution method for a geometric structure according to claim 1, characterized in that, Based on the hybrid preconditioner, the geometric model is solved to obtain the solution results of the geometric structure, including: The hybrid preconditioner is embedded into the conjugate gradient method to construct a solver; The solver iteratively processes the algebraic equations of the geometric model to obtain the solution result.
8. A finite element numerical solution system for a geometric structure, employing the finite element numerical solution method for a geometric structure as described in any one of claims 1-7, characterized in that, include: The model building module is used for: Obtain the structural parameters and boundary constraint parameters of the geometric structure, wherein the boundary constraint parameters are used to indicate the boundary conditions of the geometric structure under a displacement field or a force field; Based on the structural parameters and the boundary constraint parameters, construct the geometric model of the geometric structure; The processing module is used to perform mesh discretization processing on the geometric model to obtain the discretization parameters of the geometric model and the geometric mesh sequence constituting the geometric model, wherein the geometric mesh sequence includes multiple levels of geometric mesh; Based on the discretization parameters, a matrix-free operator for the geometric mesh is constructed using summation decomposition. Based on the matrix-free operator and the geometric mesh sequence, a hybrid preconditioner is constructed; The solver module is used to solve the geometric model according to the hybrid preconditioner to obtain the solution result of the geometric structure. The solution result includes the structural data of the geometric structure under the constraints of displacement field or force field.
9. A computer-readable storage medium, characterized in that, It is a computer-readable storage medium on which a computer program is stored, which, when executed, implements a finite element numerical solution method for a geometric structure as described in any one of claims 1-7.
10. A computer device, characterized in that, The computer device includes a processor and a memory, the memory storing at least one instruction, at least one program, code set, or instruction set, wherein at least one instruction, at least one program, code set, or instruction set is loaded and executed by the processor to implement a finite element numerical solution method for a geometric structure as described in any one of claims 1 to 7.