An efficient Catmull-Clark subdivision volume isogeometric analysis method

Abstract

The present invention discloses an efficient Catmull-Clark subdivision volume isogeometric analysis method, S1, for the CC subdivision volume initial grid input file in .inp format, the geometric information, topological information, load information, constraint information and material information of the CC subdivision volume grid are read in through a data structure; S2, the Catmull-Clark subdivision scheme is applied to subdivide the CC subdivision volume initial grid to a level suitable for analysis and calculate equivalent constraints and equivalent loads; S3, for each hexahedral unit in the subdivided CC subdivision volume grid, its CC subdivision volume unit basis function is calculated according to the grid information and unit mapping; S4, according to the CC subdivision volume unit basis function calculation method, the CC subdivision volume unit stiffness matrix after subdivision is solved; S5, the CC subdivision volume unit stiffness matrix and the linear equation group are assembled to obtain the displacement field U of all degrees of freedom, and the displacement field is visualized by writing into a vtk file. The present invention improves the speed of CC subdivision volume physical simulation under the premise of ensuring accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of an efficient Catmull-Clark subdivision volume isogeometric analysis method, and in particular to an efficient Catmull-Clark subdivision volume isogeometric analysis method. Background Art

[0002] CAD and CAE are key technologies in modern engineering design and analysis. CAD technology can realize the model design of engineering products; CAE technology uses computers to perform physical simulation to evaluate products so that designers can optimize the design of products. How to achieve the integration of industrial design and simulation analysis and perform high-precision physical simulation of complex three-dimensional solid models has been a hot issue in the intersection of multiple disciplines such as computer graphics and geometric computing in recent years. Isogeometric analysis is a new method in recent years that uses NURBS basis functions instead of the Lagrangian basis functions used in traditional finite element analysis to perform physical simulation of CAD models. It uses the same expression form in design and analysis, avoiding the grid generation process that takes about 80% of the entire "design-analysis" process, eliminating the geometric errors of finite element analysis, and has high-order continuity. Its solution accuracy is much higher than traditional finite element analysis.

[0003] Isogeometric analysis has shown its advantages over the finite element method in mechanical problems of plate and shell structures and simple three-dimensional solid models. However, it is very difficult to construct a volume parameterization suitable for isogeometric analysis for three-dimensional solid models with complex shapes and topologies. In contrast, the control grid of CC subdivision volume has great freedom, and can naturally model solid models of arbitrary topology without parameterization, which has high flexibility in the design process. CC subdivision volume can ensure C1 continuity near singular points and C2 continuity in other areas. Its continuity and numerical stability also ensure the stability of calculation when applied to physical simulation. It can be seen that CC subdivision volume is one of the most suitable representation forms for the integration of three-dimensional complex solid design and analysis. However, the related theories and tools of isogeometric analysis of volume subdivision cannot meet the requirements of physical simulation of solid models. Finite element analysis of CC subdivision volume can only analyze subdivision volume of a certain subdivision level, and its vertex and unit positions have certain errors with CC limit subdivision volume. The low-order basis function used in finite element analysis will also lead to further solution accuracy problems. Summary of the invention

[0004] One object of the present invention is to propose an efficient Catmull-Clark subdivision volume and other geometric analysis method, which improves the speed of CC subdivision volume physical simulation while ensuring accuracy.

[0005] An efficient Catmull-Clark subdivision volume isogeometric analysis method according to an embodiment of the present invention comprises the following steps:

[0006] S1. For the CC subdivision volume initial mesh input file in .inp format, the geometric information, topological information, load information, constraint information and material information of the CC subdivision volume mesh are read in through the data structure;

[0007] S2, apply the Catmull-Clark subdivision scheme to subdivide the initial mesh of the CC subdivision volume to a level suitable for analysis and calculate equivalent constraints and equivalent loads;

[0008] S3, for each hexahedral unit in the subdivided CC subdivided volume grid, calculating its CC subdivided volume unit basis function according to the grid information and unit mapping;

[0009] S4, according to the CC subdivision volume unit basis function calculation method, solving the CC subdivision volume unit stiffness matrix after subdivision;

[0010] S5. Assemble the CC subdivision volume element stiffness matrix and linear equations to obtain the displacement field U of all degrees of freedom, and write it into a vtk file to visualize the displacement field.

[0011] Optionally, the geometric information stores the vertex coordinates of each vertex in the grid, and each hexahedral unit stores information of its eight vertices;

[0012] The topological information is that for each hexahedral unit in the grid, six half-faces are stored in the order of XF, YF, ZF, XB, YB, and ZB, and each half-face stores four half-edges and the half-face opposite to it in the specified order. Each half-edge records its starting vertex and ending vertex, and stores the half-edges opposite to it in other half-faces in the current hexahedral unit and the half-edges opposite to it in the half-face opposite to the half-face where it is located, and each vertex stores the half-edge it issues;

[0013] Among them, XF, YF, and ZF are the front, right, and top when viewed from the front, and XB, YB, and ZB are the back, left, and bottom when viewed from the front;

[0014] The load information is to read in the nodes and surfaces where the load is applied as well as the direction and magnitude of the load;

[0015] The constraint information is to read in the nodes, edges and faces to which fixed constraints are imposed;

[0016] The material information is the Young's modulus and Poisson's ratio of the input material.

[0017] Optionally, S2 specifically includes giving a CC subdivision volume initial grid , using the Catmull-Clark subdivision scheme to subdivide the CC subdivision volume initial grid into the CC subdivision volume analysis level grid after k subdivisions , in the subdivision process for the CC subdivision volume initial mesh The cell numbered i in the CC subdivision analysis level grid There are The six half faces of the unit and the order of the half edges in the half faces are the same as the initial mesh of the CC subdivision volume. The direction and order of the cells numbered i in the same, and the CC subdivision volume analysis level grid The units are in a fixed order from arrive serial number.

[0018] Optionally, among the equivalent constraints and equivalent loads, for node loads, the corresponding nodes after subdivision directly inherit their loads; for surface loads, the equivalent node loads are calculated using the equivalent load method on the corresponding four surfaces after subdivision; for node constraints, the corresponding nodes after subdivision directly inherit their constraint conditions; for edge constraints, all nodes on the two edges corresponding to the edge where the constraint is located after subdivision inherit its constraint conditions; for surface constraints, all nodes on the four surfaces corresponding to the surface where the constraint is located after subdivision inherit its constraint conditions.

[0019] Optionally, the S3 specifically includes:

[0020] S31. Analysis of hierarchical grid for CC subdivision The cell i on the CC subdivision volume is calculated by calculating the initial grid of the cell i The corresponding unit , according to the initial grid of the CC subdivision volume in unit j The basis functions of unit i are calculated for the corresponding unit types above;

[0021] S32, when j is the initial grid of CC subdivision A regular CC subdivision volume cell, that is, a cell in which all edges have degree 4 and all vertices have degree 8;

[0022] At this time, in the CC subdivision volume initial grid According to its geometric and topological information, it traverses the u, v, and w directions in order. units, constituting Control vertices are encoded, and CC subdivided into volume unit basis functions is the tensor product of the cubic B-spline basis functions in the u, v, and w directions:

[0023] ;

[0024] in, are the control vertices, N(u,v,w) is the tensor product of the cubic B-spline basis functions;

[0025] S33, when j is the initial grid of CC subdivision The singular CC subdivision unit in the upper layer, that is, the unit has only one edge with degree k, and the degree of its corresponding vertex is also k, and the direction of this edge is w;

[0026] At this time, according to the CC subdivision volume initial grid Each layer of geometric and topological information on the y-axis is traversed sequentially in the u and v directions to obtain (8+k) units, and there are 3 layers in the w direction, so a total of (8+k) units are obtained. units, forming control vertices, the basis function in the w direction is a cubic B-spline basis function, and the basis functions in the u and v directions are solved using Stam's CC limit surface basis function calculation method, denoted as , then the CC volume subdivision unit basis function Depend on and N(w) give:

[0027] ;

[0028] Where N(w) is the cubic B-spline basis function, is the CC limit surface basis function;

[0029] S34, when j is the initial grid of CC subdivision Upper singular CC subdivision volume cells, i.e. cells with more than one edge whose degree is not equal to 4 and whose singular edges contain common vertices;

[0030] At this time, according to the initial grid of the subdivision The geometric and topological information on the vertices around the unit is traversed and the subdivision matrix is ​​used to solve the CC volume subdivision unit basis function :

[0031] ;

[0032] Where V is the eigenvector, is the eigenvalue, A is the subdivision matrix, P is the selection matrix, is the parameter space transformation;

[0033] S35, when j is the initial grid of CC subdivision The CC subdivision unit at the upper boundary is the initial grid of the CC subdivision The outermost unit of the boundary direction is w, and the basis functions in the u and v directions are the tensor product of the B-spline basis functions. The basis function of the CC volume subdivision unit in the w direction is is the boundary basis function C(w):

[0034] ;

[0035] Where N(u,v) is the cubic B-spline basis function tensor product, and C(w) is the boundary basis function;

[0036] S36, when j in If the cell type is not S32-S35, that is, the cell has at least 2 vertices that do not share a common edge and whose degree is not 4, it cannot be subdivided into the initial mesh of the CC volume. Analysis of the CC subdivision level grid To calculate the basis function, it is necessary to subdivide the singular vertices into different units. After subdivision, the CC subdivision volume is analyzed to form a hierarchical mesh. The unit belongs to the unit type S32-S35. At this time, the CC subdivision analysis level grid The cell i on the CC subdivision analysis level grid The geometric and topological information on the CC volume subdivision unit basis function is calculated accordingly .

[0037] Optionally, the S4 includes calculating the CC subdivision analysis level grid When the element stiffness matrix of element i is calculated, the element basis function is subdivided by CC body Calculate the strain displacement matrix B of the unit and the Jacobian matrix of the transformation between physical space and parameter space ;

[0038] For the initial mesh in CC subdivision Calculate the CC volume subdivision unit basis function The initial grid of the CC subdivision volume is obtained The position of the corresponding unit j , according to p, the Jacobian matrix of the transformation between the unit space and the parameter space is calculated ;

[0039] For the initial mesh that cannot be subdivided directly in CC The basis functions are calculated on the CC subdivision analysis level grid Calculate its CC volume subdivision unit basis function And directly calculate the Jacobian matrix of the transformation between unit space and parameter space ;

[0040] The element stiffness matrix is ​​calculated by applying Gaussian integration rule:

[0041] ;

[0042] in, are the coordinates of the Gaussian integration points, , , is the Gaussian integral weight of the Gaussian integral point, B is the strain displacement matrix, D is the material matrix, is the Jacobian matrix of the transformation between physical space and parameter space, is the Jacobian matrix of the transformation between unit space and parameter space;

[0043] For the initial mesh in CC subdivision Direct calculation of basis functions, CC subdivision analysis hierarchical grid In Each unit shares the same set of control vertices and corresponding basis functions.

[0044] Optionally, S5 specifically includes rearranging the constrained nodes and free nodes, solving the unit stiffness matrix of each unit through S4 and assembling the global stiffness matrix K according to the degree of freedom information, applying loads and constraints according to the calculated equivalent load information and equivalent constraint information, assembling to obtain the linear equation group KU=F, solving the displacement field U, and writing it into a .vtk file for visualization.

[0045] The beneficial effects of the present invention are:

[0046] The present invention uses precise geometric expressions to eliminate the errors of finite element analysis methods and approximation methods in the representation of CC subdivision volume unit basis functions, and designs data structures to increase the speed of traversing control vertices, thereby accelerating the speed of solving basis functions. For the isogeometric analysis part based on CC subdivision volume mesh mapping, by constructing a mapping relationship between subdivision grids of different levels, the high-level subdivision volume mesh unit basis functions are calculated on the low-level subdivision volume mesh, which significantly improves the unit stiffness matrix assembly speed, and improves the speed of CC subdivision volume physical simulation while ensuring accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0048] Figure 1 A flowchart of an efficient Catmull-Clark subdivision volume isogeometric analysis method proposed by the present invention;

[0049] Figure 2 : Shows the internal topological information of the unit in the data structure designed by the present invention.

[0050] Figure 3 : It shows the mapping relationship of subdivision units stipulated in the present invention.

[0051] Figure 4: The initial model of Example 1 of the present invention and its constraint and load information are shown.

[0052] Figure 5 : The subdivided model and its constraints and load information of Example 1 of the present invention are shown.

[0053] Figure 6 :The present invention demonstrates the spatial transformation when using mapping relationship to calculate the unit stiffness matrix.

[0054] Figure 7 : Shows the isogeometric analysis results of Example 1 of the present invention. DETAILED DESCRIPTION

[0055] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, which only illustrate the basic structure of the present invention in a schematic manner, and therefore only show the components related to the present invention.

[0056] refer to Figure 1 , an efficient Catmull-Clark subdivision volume isogeometric analysis method, comprising the following steps:

[0057] S1. For the CC subdivision volume initial mesh input file in .inp format, the geometric information, topological information, load information, constraint information and material information of the CC subdivision volume mesh are read in through the data structure;

[0058] S2, apply the Catmull-Clark subdivision scheme to subdivide the initial mesh of the CC subdivision volume to a level suitable for analysis and calculate equivalent constraints and equivalent loads;

[0059] S3, for each hexahedral unit in the subdivided CC subdivided volume grid, calculating its CC subdivided volume unit basis function according to the grid information and unit mapping;

[0060] S4, according to the CC subdivision volume unit basis function calculation method, solving the CC subdivision volume unit stiffness matrix after subdivision;

[0061] S5. Assemble the CC subdivision volume element stiffness matrix and linear equations to obtain the displacement field U of all degrees of freedom, and write it into a vtk file to visualize the displacement field.

[0062] In this embodiment, the geometric information stores the vertex coordinates of each vertex in the grid, and each hexahedral unit stores the information of its eight vertices;

[0063] The topological information is that for each hexahedral unit in the grid, six half-faces are stored in the order of XF, YF, ZF, XB, YB, and ZB. Each half-face stores four half-edges and the half-face opposite to it in the specified order. Each half-edge records its starting vertex and ending vertex, and stores the half-edges opposite to it in other half-faces in the current hexahedral unit and the half-edges opposite to it in the half-face opposite to its half-face. Each vertex stores the half-edge it issues.

[0064] Among them, XF, YF, and ZF are the front, right, and top when viewed from the front, and XB, YB, and ZB are the back, left, and bottom when viewed from the front;

[0065] The load information is to read in the nodes and surfaces where the load is applied as well as the direction and magnitude of the load;

[0066] The constraint information is to read in the nodes, edges and faces to which fixed constraints are imposed;

[0067] The material information is the Young's modulus and Poisson's ratio of the input material.

[0068] In this embodiment, S2 specifically includes giving a CC subdivision volume initial grid , using the Catmull-Clark subdivision scheme to subdivide the CC subdivision volume initial grid into the CC subdivision volume analysis level grid after k subdivisions , in the subdivision process for the CC subdivision volume initial mesh The cell numbered i in the CC subdivision analysis level grid There are The six half faces of the unit and the order of the half edges in the half faces are the same as the initial mesh of the CC subdivision volume. The direction and order of the cells numbered i in the same, and the CC subdivision volume analysis level grid The units are in a fixed order from arrive serial number.

[0069] Optionally, among the equivalent constraints and equivalent loads, for node loads, the corresponding nodes after subdivision directly inherit their loads; for surface loads, the equivalent node loads are calculated using the equivalent load method on the corresponding four surfaces after subdivision; for node constraints, the corresponding nodes after subdivision directly inherit their constraint conditions; for edge constraints, all nodes on the two edges corresponding to the edge where the constraint is located after subdivision inherit its constraint conditions; for surface constraints, all nodes on the four surfaces corresponding to the surface where the constraint is located after subdivision inherit its constraint conditions.

[0070] In this implementation, in equivalent constraints and equivalent loads, for node loads, the corresponding nodes after subdivision directly inherit their loads; for surface loads, the equivalent node loads are calculated using the equivalent load method on the corresponding four surfaces after subdivision; for node constraints, the corresponding nodes after subdivision directly inherit their constraint conditions; for edge constraints, all nodes on the two edges corresponding to the edge where the constraint is located after subdivision inherit its constraint conditions; for surface constraints, all nodes on the four surfaces corresponding to the surface where the constraint is located after subdivision inherit its constraint conditions.

[0071] In this implementation, S3 specifically includes:

[0072] S31. Analysis of hierarchical grid for CC subdivision The cell i on the CC subdivision volume is calculated by calculating the initial grid of the cell i The corresponding unit , according to the initial grid of the CC subdivision volume in unit j The basis functions of unit i are calculated for the corresponding unit types above;

[0073] S32, when j is the initial grid of CC subdivision A regular CC subdivision volume cell, that is, a cell in which all edges have degree 4 and all vertices have degree 8;

[0074] At this time, in the CC subdivision volume initial grid According to its geometric and topological information, it traverses the u, v, and w directions in order. units, constituting Control vertices are encoded, and CC subdivided into volume unit basis functions is the tensor product of the cubic B-spline basis functions in the u, v, and w directions:

[0075] ;

[0076] in, are the control vertices, N(u,v,w) is the tensor product of the cubic B-spline basis functions;

[0077] S33, when j is the initial grid of CC subdivision The singular CC subdivision unit in the upper layer, that is, the unit has only one edge with degree k, and the degree of its corresponding vertex is also k, and the direction of this edge is w;

[0078] At this time, according to the CC subdivision volume initial grid Each layer of geometric and topological information on the y-axis is traversed sequentially in the u and v directions to obtain (8+k) units, and there are 3 layers in the w direction, so a total of (8+k) units are obtained. units, forming control vertices, the basis function in the w direction is a cubic B-spline basis function, and the basis functions in the u and v directions are solved using Stam's CC limit surface basis function calculation method, denoted as , then the CC volume subdivision unit basis function Depend on and N(w) give:

[0079] ;

[0080] Where N(w) is the cubic B-spline basis function, is the CC limit surface basis function;

[0081] S34, when j is the initial grid of CC subdivision Upper singular CC subdivision volume cells, i.e. cells with more than one edge whose degree is not equal to 4 and whose singular edges contain common vertices;

[0082] At this time, according to the initial grid of the subdivision The geometric and topological information on the vertices around the unit is traversed and the subdivision matrix is ​​used to solve the CC volume subdivision unit basis function :

[0083] ;

[0084] Where V is the eigenvector, is the eigenvalue, A is the subdivision matrix, P is the selection matrix, is the parameter space transformation;

[0085] S35, when j is the initial grid of CC subdivision The CC subdivision unit at the upper boundary is the initial grid of the CC subdivision The outermost unit of the boundary direction is w, and the basis functions in the u and v directions are the tensor product of the B-spline basis functions. The basis function of the CC volume subdivision unit in the w direction is is the boundary basis function C(w):

[0086] ;

[0087] Where N(u,v) is the cubic B-spline basis function tensor product, and C(w) is the boundary basis function;

[0088] S36, when j in If the cell type is not S32-S35, that is, the cell has at least 2 vertices that do not share a common edge and whose degree is not 4, it cannot be subdivided into the initial mesh of the CC volume. Analysis of the CC subdivision level grid To calculate the basis function, it is necessary to subdivide the singular vertices into different units. After subdivision, the CC subdivision volume is analyzed to form a hierarchical mesh. The unit belongs to the unit type S32-S35. At this time, the CC subdivision analysis level grid The cell i on the CC subdivision analysis level grid The geometric and topological information on the CC volume subdivision unit basis function is calculated accordingly .

[0089] In this embodiment, S4 includes calculating the CC subdivision analysis level grid When the element stiffness matrix of element i is calculated, the element basis function is subdivided by CC body Calculate the strain displacement matrix B of the unit and the Jacobian matrix of the transformation between physical space and parameter space ;

[0090] For the initial mesh in CC subdivision Calculate the CC volume subdivision unit basis function The initial grid of the CC subdivision volume is obtained The position of the corresponding unit j , according to p, the Jacobian matrix of the transformation between the unit space and the parameter space is calculated ;

[0091] For the initial mesh that cannot be subdivided directly in CC The basis functions are calculated on the CC subdivision analysis level grid Calculate its CC volume subdivision unit basis function And directly calculate the Jacobian matrix of the transformation between unit space and parameter space ;

[0092] The element stiffness matrix is ​​calculated by applying Gaussian integration rule:

[0093] ;

[0094] in, are the coordinates of the Gaussian integration points, , , is the Gaussian integral weight of the Gaussian integral point, B is the strain displacement matrix, D is the material matrix, is the Jacobian matrix of the transformation between physical space and parameter space, is the Jacobian matrix of the transformation between the unit space and the parameter space;

[0095] For the initial mesh in CC subdivision Direct calculation of basis functions, CC subdivision analysis hierarchical grid In Each unit shares the same set of control vertices and corresponding basis functions.

[0096] In this embodiment, S5 specifically includes rearranging the constrained nodes and free nodes, solving the unit stiffness matrix of each unit through S4 and assembling the global stiffness matrix K according to the degree of freedom information, applying loads and constraints according to the calculated equivalent load information and equivalent constraint information, assembling the linear equation group KU=F, solving the displacement field U, and writing it into the .vtk file for visualization.

[0097] Embodiment 1:

[0098] first step:

[0099] For the CC subdivision volume initial mesh input file in ".inp" format, the geometry, topology information, load, constraint, material and other physical simulation information of the CC subdivision volume mesh are read in through the self-designed data structure, so as to realize the fast fixed order traversal of the control vertices of the CC subdivision volume unit basis function, and facilitate the calculation of the mapping matrix from the parameter space to the unit space. For the specific order, see Figure 2 .

[0100] CC subdivision volume initial mesh Contains 1612 nodes, 1196 cells, reads in the CC subdivision volume initial grid Material, load, and constraint information: Initial mesh of CC subdivision volume Apply 4 node loads (concentrated forces) in the z direction with a value of 0.1 to the parts shown in the figure; The parts shown in the figure are constrained to 9 faces. CC subdivision volume initial mesh The corresponding model and its loads and constraints are as follows Figure 4 shown.

[0101] Step 2: Apply the Catmull-Clark subdivision scheme to subdivide the mesh to a level suitable for analysis and calculate equivalent constraints and equivalent loads, as follows:

[0102] Subdividing the initial mesh in CC Apply the Catmull-Clark subdivision scheme once to obtain the CC subdivision volume analysis level grid , CC subdivision analysis level grid Contains 11097 vertices and 9568 cells. CC subdivision volume initial mesh The cell numbered i in the figure is subdivided into CC subdivision analysis level grids The cells numbered 8i to 8i+7 in the CC subdivision analysis level grid The order of half faces and half edges of the 8 units and the initial mesh of the CC subdivision volume The units i in the equation are exactly the same, which makes it easy to calculate the Jacobian matrix of the transformation between the unit space and the parameter space, such as Figure 3 .

[0103] For CC subdivision initial mesh The node loads in the calculation are calculated on the CC subdivision analysis level grid The equivalent load in is the node load with a value of 0.1 in the z direction on the corresponding node after 4 subdivisions; for the initial grid of the CC subdivision volume The surface fixed constraint in the calculation is calculated at the CC subdivision analysis level mesh The equivalent constraints in are the node fix constraints on the corresponding surface after 49 subdivisions. CC subdivision volume analysis level grid The corresponding model and its calculated equivalent loads and constraints are as follows: Figure 5 shown.

[0104] Step 3: For different types of CC subdivision volume units after subdivision, different rules are used to calculate their subdivision basis functions according to the unit mapping and CC subdivision volume basis function calculation, as follows:

[0105] For CC subdivision initial mesh It does not belong to the four cases where basis function units can be calculated, such as Figure 4 As shown in the CC subdivision analysis level grid The corresponding units in are all units that can directly calculate the basis function. Calculate its basis function on the initial grid of CC subdivision volume The basis function unit i can be calculated in the calculation of CC subdivision analysis level grid When the corresponding 8 element stiffness matrices are obtained, the CC subdivision analysis level grid can be obtained in one traversal The control vertices of the previous 8 units are obtained to obtain the common subdivision basis function representation of the 8 units.

[0106] Step 4: According to the CC subdivision body parameterization method and mapping relationship, solve the subdivided CC subdivision body unit stiffness matrix, as follows:

[0107] The calculation formula of the element stiffness matrix is ​​as follows:

[0108] ;

[0109] in is the mapping matrix from physical space to parameter space, It is the mapping matrix from parameter space to unit space, which is calculated based on the different positions of the subdivided units. By subdividing the basis function definition:

[0110] ;

[0111] The constitutive matrix D shared by all elements of the grid is calculated based on Young's modulus and Poisson's ratio, and the element stiffness matrix can be obtained through Gaussian integration.

[0112] Step 5: Assemble the global stiffness matrix and linear equations and solve them as follows:

[0113] According to the degree of freedom number of the above element stiffness matrix, the global stiffness matrix K is assembled. And according to the linear elastic problem equation system KU=F, for the CC subdivision body analysis level grid The equivalent node constraints on U are fixed to 0. For the CC subdivision analysis level grid The equivalent node load on F is set to the load value. Use SuperLU to solve the linear equations and get the final displacement U of all degrees of freedom. Write it into the "vtk" file for visualization, such as Figure 7 shown.

[0114] The present invention uses precise geometric expressions to eliminate the errors of finite element analysis methods and approximation methods in the representation of CC subdivision volume unit basis functions, and designs data structures to increase the speed of traversing control vertices, thereby accelerating the speed of solving basis functions. For the isogeometric analysis part based on CC subdivision volume mesh mapping, by constructing a mapping relationship between subdivision grids of different levels, the high-level subdivision volume mesh unit basis functions are calculated on the low-level subdivision volume mesh, which significantly improves the unit stiffness matrix assembly speed, and improves the speed of CC subdivision volume physical simulation while ensuring accuracy.

[0115] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. An efficient Catmull-Clark subdivision volume isogeometric analysis method, characterized in that: The steps include: S1. For the CC subdivision volume initial mesh input file in .inp format, the geometric information, topological information, load information, constraint information and material information of the CC subdivision volume mesh are read in through the data structure; S2, apply the Catmull-Clark subdivision scheme to subdivide the initial mesh of the CC subdivision volume to a level suitable for analysis and calculate equivalent constraints and equivalent loads; S3, for each hexahedral unit in the subdivided CC subdivided volume grid, calculating its CC subdivided volume unit basis function according to the grid information and unit mapping; S4, according to the CC subdivision volume unit basis function calculation method, solving the CC subdivision volume unit stiffness matrix after subdivision; S5, assemble the CC subdivision volume unit stiffness matrix and linear equations to obtain the displacement field U of all degrees of freedom, and write it into a vtk file to visualize the displacement field; S4 includes calculating the CC subdivision analysis level grid When the element stiffness matrix of element i is calculated, the element basis function is subdivided by CC body Calculate the strain displacement matrix B of the unit and the Jacobian matrix of the transformation between physical space and parameter space ; For the initial mesh in CC subdivision Calculate the CC volume subdivision unit basis function The initial grid of the CC subdivision volume is obtained The position of the corresponding unit j , according to p, the Jacobian matrix of the transformation between the unit space and the parameter space is calculated ; For the initial mesh that cannot be subdivided directly in CC The basis functions are calculated on the CC subdivision analysis level grid Calculate its CC volume subdivision unit basis function And directly calculate the Jacobian matrix of the transformation between unit space and parameter space ; The element stiffness matrix is ​​calculated by applying Gaussian integration rule: ; in, are the coordinates of the Gaussian integration point, , , is the Gaussian integral weight of the Gaussian integral point, B is the strain displacement matrix, D is the material matrix, is the Jacobian matrix of the transformation between physical space and parameter space, is the Jacobian matrix of the transformation between unit space and parameter space; For the initial mesh in CC subdivision Direct calculation of basis functions, CC subdivision analysis hierarchical grid In Each unit shares the same set of control vertices and corresponding basis functions.

2. The efficient Catmull-Clark subdivision volume isogeometric analysis method according to claim 1, characterized in that: The geometric information stores the vertex coordinates of each vertex in the grid, and each hexahedral unit stores the information of its eight vertices; The topological information is that for each hexahedral unit in the grid, six half-faces are stored in the order of XF, YF, ZF, XB, YB, and ZB, and each half-face stores four half-edges and the half-face opposite to it in the specified order. Each half-edge records its starting vertex and ending vertex, and stores the half-edges opposite to it in other half-faces in the current hexahedral unit and the half-edges opposite to it in the half-face opposite to the half-face where it is located, and each vertex stores the half-edge it issues; Among them, XF, YF, and ZF are the front, right, and top when viewed from the front, and XB, YB, and ZB are the back, left, and bottom when viewed from the front; The load information is to read in the nodes and surfaces where the load is applied as well as the direction and magnitude of the load; The constraint information is to read in the nodes, edges and faces to which fixed constraints are imposed; The material information is the Young's modulus and Poisson's ratio of the input material.

3. The efficient Catmull-Clark subdivision volume isogeometric analysis method according to claim 1, characterized in that: S2 specifically includes giving a CC subdivision volume initial grid , using the Catmull-Clark subdivision scheme to subdivide the CC subdivision volume initial grid into the CC subdivision volume analysis level grid after k subdivisions , in the subdivision process for the CC subdivision volume initial mesh The cell numbered i in the CC subdivision analysis level grid There are The six half faces of the unit and the order of the half edges in the half faces are the same as the initial mesh of the CC subdivision volume. The direction and order of the cells numbered i in the same, and the CC subdivision volume analysis level grid The units are in a fixed order from arrive serial number.

4. The efficient Catmull-Clark subdivision volume isogeometric analysis method according to claim 1, characterized in that: Among the equivalent constraints and equivalent loads, for node loads, the corresponding nodes after subdivision directly inherit their loads; for surface loads, the equivalent node loads are calculated using the equivalent load method on the corresponding four surfaces after subdivision; for node constraints, the corresponding nodes after subdivision directly inherit their constraint conditions; for edge constraints, all nodes on the two edges corresponding to the edge where the constraint is located after subdivision inherit its constraint conditions; for surface constraints, all nodes on the four surfaces corresponding to the surface where the constraint is located after subdivision inherit its constraint conditions.

5. The efficient Catmull-Clark subdivision volume isogeometric analysis method according to claim 3, characterized in that: The S3 specifically includes: S31. Analysis of hierarchical grid for CC subdivision The cell i on the CC subdivision volume is calculated by calculating the initial grid of the cell i The corresponding unit , according to the initial grid of the CC subdivision volume in unit j The basis functions of unit i are calculated for the corresponding unit types above; S32, when j is the initial grid of CC subdivision A regular CC subdivision volume cell, that is, a cell in which all edges have degree 4 and all vertices have degree 8; At this time, in the CC subdivision volume initial grid According to its geometric and topological information, it traverses the u, v, and w directions in order. units, constituting Control vertices are encoded, and CC subdivided into unit basis functions is the tensor product of the cubic B-spline basis functions in the u, v, and w directions: ； in, are control vertices, N(u,v,w) is the tensor product of cubic B-spline basis functions, u, v, w are the three directions of parameter space; S33, when j is the initial grid of CC subdivision The singular CC subdivision unit in the upper layer, that is, the unit has only one edge with degree k, and the degree of its corresponding vertex is also k, and the direction of this edge is w; At this time, according to the CC subdivision volume initial grid Each layer of geometric and topological information on the y-axis is traversed sequentially in the u and v directions to obtain (8+k) units, and there are 3 layers in the w direction, so a total of (8+k) units are obtained. units, forming control vertices, the basis function in the w direction is a cubic B-spline basis function, and the basis functions in the u and v directions are solved using Stam's CC limit surface basis function calculation method, denoted as , then the CC volume subdivision unit basis function Depend on and N(w) give: ; Where N(w) is the cubic B-spline basis function, is the CC limit surface basis function; S34, when j is the initial grid of CC subdivision Upper singular CC subdivision volume cells, i.e. cells with more than one edge whose degree is not equal to 4 and whose singular edges contain common vertices; At this time, according to the initial grid of the subdivision The geometric and topological information on the vertices around the unit is traversed and the subdivision matrix is ​​used to solve the CC volume subdivision unit basis function : ; Where V is the eigenvector, is the eigenvalue, A is the subdivision matrix, P is the selection matrix, is the parameter space transformation; S35, when j is the initial grid of CC subdivision The CC subdivision unit at the upper boundary is the initial grid of the CC subdivision The outermost unit of the boundary direction is w, and the basis functions in the u and v directions are the tensor product of the B-spline basis functions. The basis function of the CC volume subdivision unit in the w direction is is the boundary basis function C(w): ; Where N(u,v) is the cubic B-spline basis function tensor product, and C(w) is the boundary basis function; S36, when j in If the cell type is not S32-S35, that is, the cell has at least 2 vertices that do not share a common edge and whose degree is not 4, it cannot be subdivided into the initial mesh of the CC volume. Analysis of the CC subdivision level grid To calculate the basis function, it is necessary to subdivide the singular vertices into different units. After subdivision, the CC subdivision volume is analyzed to form a hierarchical mesh. The unit belongs to the unit type S32-S35. At this time, the CC subdivision analysis level grid The cell i on the CC subdivision analysis level grid The geometric and topological information on the CC volume subdivision unit basis function is calculated accordingly .

6. The efficient Catmull-Clark subdivision volume isogeometric analysis method according to claim 1, characterized in that: The S5 specifically includes rearranging the constrained nodes and free nodes, solving the unit stiffness matrix of each unit through S4 and assembling the global stiffness matrix K according to the degree of freedom information, applying loads and constraints according to the calculated equivalent load information and equivalent constraint information, assembling to obtain the linear equation group KU=F, solving the displacement field U, and writing it into the .vtk file for visualization.

Images