A Method for Generating Complex Multi-Constraint Structural Grids Based on Frame Fields

Abstract

The present invention discloses a frame field-based complex multi-constraint structured grid generation method for models with complex geometric feature constraints and dimensional constraints. Existing quadrilateral grid generation methods are difficult to achieve a method with high robustness while satisfying complex constraints and ensuring high-quality grid results. The frame field-guided grid method adopted in the present invention has a high-quality guarantee on the one hand. In addition, the frame field only plays a certain role in the generation of the initial streamline, and subsequent methods are independent of the frame field, avoiding the situation where the solution of the frame field is not applicable to grid generation. In this method, under various complex constraint conditions tested, the grid can be generated with high-quality results while ensuring geometric constraints and satisfying dimensional constraints.

Description

Technical Field

[0001] The present invention belongs to the field of computer-aided design, and particularly relates to a method for generating quasi-structured quadrilateral meshes for CAD models. Background Art

[0002] Regarding the generation of structured quadrilateral meshes, generally there are two main fields of research. One is engineering simulation analysis. Quadrilateral meshes have better geometric support in finite element and finite volume methods compared to triangular meshes. The other is the field of computer graphics. Besides numerical simulation, quadrilateral meshes are also very useful in surface modeling and texturing.

[0003] In recent decades, the engineering analysis field has been dedicated to the generation of quadrilateral meshes, resulting in two major categories of industrial technologies. One is to manually divide quadrilateral regions using some semi-automatic assistance methods, and then generate standardized quadrilateral mesh filling regions according to user needs. Such a method is very time-consuming and is mainly used for demanding numerical simulations such as CFD. The other is some fully automatic methods such as the advancing front method based on a frame field to generate right-angled triangular meshes, and then perform operations such as merging and subdivision to generate full quadrilateral meshes. However, such methods have a large number of irregular vertices and lack user control.

[0004] On the other hand, the field of computer graphics has thoroughly developed and explored the quadrilateral mesh partitioning technology based on a frame field. The core is to automatically generate quadrilateral topological partitions. However, such topological partitions are not necessarily guaranteed to be extractable, and such technology lacks sufficient robustness and stability when applied to engineering analysis.

[0005] The quasi-quadrilateral mesh partitioning proposed in this paper divides topological partitions based on a frame field, but does not require them to be necessarily quadrilateral topological partitions. In such topological partitions, the template method is used to generate meshes for each partition, and the robustness is guaranteed to a certain extent. After testing, it can basically meet the user's specified size field requirements, and the generated meshes have few singular points and high quality. Summary of the Invention

[0006] The purpose of the present invention is to overcome the deficiencies of the prior art and propose a method for generating a complex multi-constraint quadrilateral mesh structure based on finite element frame field calculation for CAD models. According to the complex geometric feature constraints and size constraints, the frame field in the mesh is topologically extracted, reasonable topological simplification is performed, and the template method is used to generate meshes for the topology.

[0007] The present invention is specifically as follows:

[0008] Step 1: First, convert the input background triangular mesh and the size constraints in the corresponding size constraints and geometric feature constraints into corresponding geometric constraints;

[0009] 1.1 Use the marching triangle method to extract the isocontours for each triangular face;

[0010] 1.2 Connect the extracted isocontours and perform a certain degree of smoothing;

[0011] 1.3 Filter the isocontours with a small number;

[0012] 1.4 Re - mesh at the isocontours and store the isocontours as geometric feature constraints;

[0013] Step 2: Use the finite element method to calculate the frame field of the background mesh and extract the singular points therein;

[0014] Cross representation:

[0015] Use the standard unit vector to represent u = (cos(4θ), sin(4θ)) = (u1, u2), where θ is the angle of one of the four branches and the local reference range. Contrary to the general angle representation θ, the vector representation u is invariant under quadrilateral symmetry and is suitable for linear finite - element interpolation. Since the linear combination of unit vectors is usually not a unit vector, the interpolation within a triangle is not a strict cross - representation, but the cross can be easily recovered by projecting onto the unit circle, that is, normalizing the vector representation.

[0016] Finite - element discretization:

[0017] To discretize the vector - field representation of the cross - field, define a cross at each edge of the triangulation and use Crouzeix - Raviart interpolation, as shown in. The unknowns of the cross - field problem are the representation vector - field components u1, u2 at each edge eij. Since it is desired that the cross - field aligns with the surface boundary and the characteristic lines, the Dirichlet boundary condition is u = (cos(0), sin(0)) = (1, 0) at each boundary and characteristic line.

[0018] Smooth the frame field using the heat - conduction equation

[0019] To have a smooth cross - field within the domain, a natural approach is to minimize the Dirichlet energy:

[0020]

[0021] Since This objective function is non - linear. By directly solving the Laplacian of the vector representation (i.e., )It is clearly inappropriate to minimize the Dirichlet energy because these values may collapse to 0 far from the boundary. The solution to this problem is to solve a Ginzburg-Landau nonlinear equation to penalize values that leave the cross manifold, but this method is too time-consuming for complex problems. An effective alternative is the MBO method, in which heat diffusion (Equation 1) and projection (Equation 2) steps are alternated. This method uses decreasing diffusion time steps, selects them based on grid size considerations, and groups them again to allow reuse of the matrix factorization computed by a direct linear solver.

[0022] Diffusion:

[0023] Projection:

[0024] Detecting singular points:

[0025] For each point x on the surface, cross-field singularities are detected by calculating the angular difference along a small closed circle γ centered at x. This defines a singularity index:

[0026]

[0027] By definition, when the index an of a cross-field singularity is -1, it corresponds to an irregular vertex with index k = -1 in a quadrilateral mesh, i.e., valence five.

[0028] In the discrete setting, extract the oriented edge single loop (ei1,..., ein) around vertex i and calculate the sum of the angular differences. In practice, there are three cases: the sum is zero and the vertex is regular, or the sum is one or negative one and the vertex is a singularity of the frame field. Due to the Crouzeix-Raviart discretization of the frame field (one angle per edge), cross-field singularities may be located at the vertices, edges, or triangles of the mesh. When a singularity is located on an edge, the indices of its two adjacent vertices are not zero, and when a singularity is located on a triangle, the indices of its three adjacent vertices are not zero.

[0029] Step 3. Implement the Runge-Kutta method at the singular points and the calculated partial boundary points to extend the streamlines;

[0030] 3.1. First, calculate the first discrete point of the starting point in a one-neighborhood and a two-neighborhood at the above points:

[0031] ① First, calculate the first discrete point in the relevant neighborhood of the boundary singular point

[0032] ② Calculate whether there is an internal singular point at the boundary singular point. If so, execute ③, otherwise jump to ④;

[0033] ③If the singular point is within the surface, extend this discrete point to the edge / point of the triangle; otherwise, directly use it as a discrete point, and then delete this singular point;

[0034] ④Screen all the discrete points obtained above that are related to the boundary singular points, and remove the discrete points with similar directions

[0035] ⑤Calculate the first discrete point in the relevant domain of the remaining internal singular points

[0036] ⑥If there are other singular points in the second domain of a singular point, if so, execute ⑦; otherwise, jump to ⑧

[0037] ⑦Select a singular point, and use the other singular point as a discrete point of this singular point. If this discrete point is within the surface, extend it to the edge / point of the triangle; otherwise, directly use it as a discrete point, and then delete the other singular point;

[0038] ⑧Screen all the discrete points obtained above that are related to the internal singular points, and remove the discrete points with similar directions

[0039] 3.2 Extend the streamline of the obtained singular points and the corresponding first discrete points according to the Runge-Kutta method. The termination conditions are as follows:

[0040] ①Until the first domain of the singular point, terminate with the singular point

[0041] ②Terminate with the nearest triangular grid point at the termination of the constraint edge / boundary edge;

[0042] ③The same streamline passes through the same background grid triangle in the same direction and terminates.

[0043] Step 4. Simplify the streamline and construct the topological relationship:

[0044] 4.1 For the streamline structure obtained in Step 3, detect its streamline situation and perform different simplifications for different situations until no simplifiable situation appears. All simplification situations are roughly divided into three cases:

[0045] ①Two streamlines have the same start and end points, denoted as A and B respectively, and the lengths are quite different. Modify the longer streamline to A->A, and keep the shorter streamline unchanged

[0046] ②Two streamlines have the same start and end points, and the lengths are not much different. Remove one of the streamlines

[0047] ③One of the start and end points of two streamlines is the same, denoted as A, and the additional points are B and C respectively. Without loss of generality, assume A->B is longer than A->C. Modify the longer streamline to B->C;

[0048] 4.2 Construct the topological partition relationship according to the final streamline structure.

[0049] Step 5. Topological simplification of chord folding: Based on the dual structure of the quadrilateral mesh and keeping the number of singular points and the basic topological relationships in the mesh unchanged, a topological simplification operation and optimization framework for quadrilateral meshes is designed, with high algorithm execution efficiency and good simplification effect, and it can be applied to the simplification of various quadrilateral meshes.

[0050] The definitions of chord, horizontal edge, and vertical edge are as follows:

[0051] Chord: In the topological division of the quadrilateral region, the quadrilateral region formed by the set of opposite sides corresponding to each side.

[0052] Horizontal edge: All the opposite sides from the boundary to the opposite side of another boundary within the chord.

[0053] Vertical edge: The remaining edges within the chord except for the horizontal edges (the edges connecting the vertices of the horizontal edges).

[0054] Chord folding operation:

[0055] Delete an entire chord from the dual representation of the model, which means deleting all the quadrilateral patches passed by this chord. The specific deletion operation is completed by folding all the horizontal edges passed by the chord, that is, merging the vertices of the chord horizontal edges pairwise, and merging the regional lines on the corresponding vertical edges pairwise. The dual chord folding in the quadrilateral mesh does not affect the connectivity of the mesh, and it can ensure that the folded mesh is still a pure quadrilateral mesh.

[0056] Folding operation of the quadrilateral region:

[0057] The folding of the quadrilateral region is a topological region deletion operation achieved by merging the two opposite vertices of the multi - block regions formed by combining some quadrilateral regions within the chord. Compared with the folding of the dual chord, this is a local mesh simplification operation. The quadrilateral region can be imagined as two triangles, and the edge between the vertices to be merged is the common edge, and the folding of the quadrilateral can be regarded as the folding of the edge of the triangle. However, this common edge is linearly interpolated and synthesized by the two vertical edges of the chord. Since the vertical edge is not a simple straight line, but is composed of multiple regional lines and the number of points on the edge is inconsistent, the interpolation process also needs to fit the same number of points for each regional line. The corresponding calculation formula is as follows:

[0058] q i =(1 - t i )×a i +t i ×b i ,q i+1 =(1 - t i+1 )×a i+1 +t i+1 ×b i+1

[0059]

[0060]

[0061] where \(i\) represents the serial number of the regional line, \(a\) i a i+1 , \(b\) i b i+1 are the regional lines on the longitudinal side within the chord, \(a\) i b i are the transverse sides within the chord. The regional line is composed of many broken lines, \(j\) represents the serial number of the broken line on the regional line, \(t\) i is x i,j a point on the longitudinal side regional line

[0062] The folding operation influence range of the quadrilateral region is only the region connected to the folded quadrilateral, and the whole operation satisfies the restriction of not introducing non - quadrilateral regions.

[0063] Step 5 includes the following sub - steps:

[0064] ①: In the topological partition of the quadrilateral region, divide the chord according to the opposite sides (search for the opposite side from the boundary side until another boundary side or start from this side and go around a circle back to this side)

[0065] ②: Classify the chords by type:

[0066] (1) There are no singular points or singular points exist on the same longitudinal side

[0067] (2) There are two singular points on the same transverse side

[0068] (3) Singular points exist on both longitudinal sides and are not on the same transverse side

[0069] ③: Check the three types of chords:

[0070] (1) If the longitudinal side to be merged is a boundary side or a constraint side, discard this chord.

[0071] (2) The first - type chords do not need to be processed, and the second - and third - type chords are transferred to step four.

[0072] ④: Sort the selected chords according to the energy formula of the priority metric, and select the chord with the maximum energy and greater than 0 for topological region simplification. If there is no chord with energy greater than 0, end the simplification.

[0073] The priority metric includes:

[0074] Use the energy formula to select and delete those long and narrow and irregular quadrilateral regions, and simplify a quadrilateral region topological dissection structure into an optimal topological structure with the fewest singular points; the energy formula is as follows:

[0075]

[0076] In the formula, K controls the chord energy ratio;

[0077] ⑤: Check the type of the chord to perform region simplification:

[0078] (1) For the second type of chord, directly perform the folding operation on the quadrilateral region;

[0079] (2) The third type of chord will be divided into non-zip-patch regions and zip-patch regions according to the positions of the singular points (zip-patch is the region between adjacent opposite-corner singular points)

[0080] For the non-zip-patch region, delete the inner opposite-side longitudinal edge and points of the chord closest to its singular point, and then merge it with the quadrilateral region outside the longitudinal edge. For the zip-patch region, linearly interpolate and fit a new longitudinal edge according to the original transverse and longitudinal edges, delete the original two longitudinal edges, and then divide the transverse edge into two according to the new points and merge it with the outside of the chord.

[0081] (3) Return to ①.

[0082] Step 6: For the above-obtained region division, calculate the number of segments of each region edge based on the background size length, and perform certain optimization on this number of segments:

[0083] Since generating a mesh using the template method depends on the number of segments set for each edge, region line quantization is very important. First, perform mesh dissection on the model region lines, and then generate templates for the model region surfaces. Limited by curve discretization, this works well for generating triangular meshes (or tetrahedral meshes), but for quadrilateral meshes (or hexahedral meshes), the quadrilateral topological constraints (chords) are global and pass through the model curves.

[0084] Quantifying the region lines is very simple. Use the integer rounded value of the ideal number of edges calculated by integrating the global size field, and apply equality to the opposite sides, except for the region lines of the quadrilateral regions. This non-optimal quantization will result in many necessary dipoles in the final quadrilateral meshes of different region surfaces.

[0085] Step 6 includes the following sub-steps:

[0086] Sub-step 1, using the quantization calculation method, calculate the integer rounded value of the ideal number of edges for all region lines in combination with the size field; the quantization calculation method for the number of segments on the region lines:

[0087] Sub-step 2: Using the topological constraints of quadrilateral region quantization, impose equality on the number of segments of some quadrilateral regions;

[0088] Sub-step 3: Perform secondary quantization on all region lines with the obtained final ideal number of edges.

[0089] Quantization calculation of the number of segments on the region line

[0090] Consider a region line C, which is parameterized by t varying from a to b. The floating-point number of edges on the curve according to the dimension map is: To obtain an integer number of edges n e , use a simple rounding:

[0091] The parameter t associated with the i-th interior point of the region line i , i ∈ [1, n e - 1] is such that:

[0092]

[0093] This integral equation can be solved by numerical integration, by adding values along the region line until the sum equals Perform linear interpolation between samples. To mesh the region line C, calculate the parameter t using (Equation 11) i and evaluate the region line parameterization: x i = f(t i ) to calculate the vertex positions x i , i ∈ [1, n e - 1]. Using this method, the points are well placed on the curve according to the dimension map, smoothly transitioning from smaller feature regions to coarser regions.

[0094] Topological constraints of quadrilateral regions

[0095] The integration of the dimension field does not consider the particularity of the quadrilateral mesh topology, so the quadrilaterals are organized into topological chords (double loops of adjacent elements). In this work, the number of edges n on the opposite sides of the quadrilateral region face is adjusted e . Consider a region face whose boundary is composed of four region lines (C1, C2, C3, C4), and force the number of points on the opposite sides to be equal, i.e., n e1 = n e3 and n e2 = n e4 , unless the integral values are very different. When two adjacent quadrilateral region faces share a region line C c , the n on the two faces ecThe values must be the same. This means that the equality constraints are being propagated. To solve the propagation problem, topological strings associated with the faces of the quadrilateral regions are constructed. The strings consist of topologically parallel region lines, all of which receive the same fixed number of points, and these points are calculated by averaging the ideal values previously computed on each curve. Through this simple propagation, the string propagates only when two adjacent quadrilateral faces share a curve, which is one of their four edges.

[0096] Step 7: Use the template method to generate a mesh for each partition. Steps:

[0097] 1) Input the corner coordinates of the polygon region and the subdivision numbers of the corresponding sides in clockwise / counterclockwise order.

[0098] 2) Select a suitable topological template according to the input subdivision numbers, and determine the topological structure of the patch by calculating the corresponding parameters.

[0099] 3) Calculate the coordinates of the points on the boundary, and use Laplacian smoothing to determine the geometric position of the patch.

[0100]

[0101] subject to v i = ω i , i ∈ C

[0102] where ε is the set of edges of the mesh, C is the set of boundary points, and ω i ∈ R 2 is the coordinate of the i-th fixed boundary point input.

[0103] Simplify the input:

[0104] Since considering all possible input cases is a challenge, the problem can be considered to be reduced to an equivalent and easier-to-solve sub-problem. This operation greatly reduces the cases that need to be considered, thereby ensuring the generality of the algorithm.

[0105] For the input (l0,..., l N-1 ), for some k, l k-1 and l k+1 are both greater than 1. Define d = min(l k-1 , l k+1 ) - 1. Then the input can be simplified to (l0,..., l N-1 ), where

[0106]

[0107] The above simplification operation can be repeated until no further simplification is possible, thus obtaining the maximally simplified input.

[0108] The template selection is formulated as an ILP:

[0109] Each topological template defines a linear relationship between the inputs (l0,..., l N-1 ) and the parameters, where the parameters include the padding amount p i of the i-th edge, the number of edge flows x, and y.

[0110] The above general form is denoted as:

[0111] Ax = b

[0112] where A is an N×M matrix, where M is the number of parameters of the template, x represents the M-dimensional vector of template parameters (to be solved), b is an N-dimensional vector jointly determined by the inputs (l0,..., l N-1 ), and it is required that be non-negative integers, and it can be expressed as an integer linear programming (ILP):

[0113] argmax c T x

[0114] subject to Ax = b

[0115] x ≥ 0, x ∈ Z M

[0116] where c represents the M-dimensional vector of the objective to be maximized. As long as the ILP has a solution, it means that the template is feasible for the input. Since M, the number of parameters of a certain template, is very small (at most 10), it can be solved quickly. Since N ≤ M, there are usually multiple solutions to the above ILP, which means that for the same topological template, there are multiple parameter combinations that meet the input requirements. Setting the objective to maximize Under this objective, when the corresponding variable represents the boundary padding amount, c j = 1; otherwise, c j = 0, and the number of singular points included in the topological generation at this time is the least.

[0117] Given the region division result of the problem domain and the number of subdivisions of the corresponding edges of each sub-region, the quadrilateral mesh generation of each sub-region is obtained using the above algorithm, and then they are stitched together to form the final mesh effect.

[0118] Step 8, Mesh smoothing:

[0119] ① Regular vertices: Given a fixed boundary, solve the Winslow non-linear elliptic PDE Δ xu = 0, where u(x) is a certain coordinate calculation space and x is the coordinate in the physical space. The advantage of this is that the two coordinate components are coupled, and as a result, the quadrilateral shape is very good, and some orthogonality can be enforced even under large distortions. By applying finite difference (FDM) discretization to the Winslow equation, the kernel of the regular vertices in the locally smoothed quadrilateral mesh can be derived. Assume that (x1,..., x8) are the vertices of the template orderly surrounding the regular vertex x, and its new position is given by:

[0120]

[0121] where α0 = (x3 - x7)·(x3 - x7), α1 = (x1 - x5)·(x1 - x5), and β = (x1 - x5)·(x3 - x7)

[0122] ② Irregular vertices: Use the angle-based smoothing method. The idea is to move the vertices on the single-ring bisector. As the Winslow kernel above, it works well when the mesh is not too restricted, but it may create flipped cells in complex configurations.

[0123] The substantial effect of the present invention lies in that the robustness is guaranteed, the user's specified size field requirements are met, and the generated mesh has few singular points and high quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0124] Figure 1 is the flow chart of the present invention;

[0125] Figure 2-1 is the triangular background mesh diagram of a model;

[0126] Figure 2-2 is Figure 2-1 the constraint diagram of the geometric feature lines of the model;

[0127] Figure 2-3 is Figure 2-1 the background size constraint diagram of the model;

[0128] Figure 3 is the schematic diagram of the frame field calculation result and the singular point distribution;

[0129] Figure 4 is the schematic diagram of the initial streamline result;

[0130] Figure 5 is the schematic diagram of the topological structure after simplifying the streamline;

[0131] Figure 6 is the schematic diagram of the simplified topological structure;

[0132] Figure 7 is the schematic diagram of the initial background quadrilateral mesh generated according to the final topology;

[0133] Figure 8 Schematic diagram of the final quadrilateral mesh after grid geometry optimization. Specific implementation manner

[0134] The present invention will be further described below in conjunction with the accompanying drawings.

[0135] As Figure 1 shown, for the quadrilateral mesh generation method for a complex constraint model, the specific steps are as follows:

[0136] Step 1: First, convert the size constraints in the input background triangular mesh and the corresponding size constraints and geometric feature constraints into corresponding geometric constraints;

[0137] 1.1 Use the marching triangle method to extract the isocontours for each triangular face;

[0138] 1.2 Connect the extracted isocontours and perform a certain degree of smoothing;

[0139] 1.3 Filter the isocontours with a small number;

[0140] 1.4 Re-mesh at the isocontours and store the isocontours as geometric feature constraints;

[0141] It is worth mentioning that if the geometric feature constraints can include the places where the mesh size changes greatly in this step, there is no need to add such isocontours specifically. In this example, no isocontours are added. Figure 2-1 Schematic diagram of the background triangular mesh, Figure 2-2 Schematic diagram of the geometric feature line constraint, Figure 2-3 Background size field constraint.

[0142] Step 2: Use the finite element method to calculate the frame field of the background mesh and extract the singular points therein;

[0143] Cross representation:

[0144] Use the standard unit vector to represent u = (cos(4θ), sin(4θ)) = (u1, u2), where θ is the angle of one of the four branches and the local reference range. Contrary to the general angle representation θ, the vector representation u is invariant to the quadrilateral symmetry and is suitable for linear finite element interpolation. Since the linear combination of unit vectors is usually not a unit vector, the interpolation within a triangle is not a strict cross representation, but the cross can be easily restored by projecting onto the unit circle, that is, normalizing the vector representation. Since it is desired that the cross field aligns with the surface boundary and the feature lines, the Dirichlet boundary condition is u = (cos(0), sin(0)) = (1, 0) at each boundary and feature line.

[0145] Smoothing with the heat equation

[0146] To have a smooth cross-field within the domain, the natural approach is to minimize the Dirichlet energy:

[0147]

[0148] In this work, decreasing diffusion time steps are used, chosen according to grid size considerations and regrouped to allow reuse of the matrix factorization computed by the direct linear solver.

[0149] Diffusion:

[0150] Projection:

[0151] Detecting singular points:

[0152] For each point x on the surface, cross-field singularities are detected by computing the angular difference along a small closed circle γ centered at x. This defines a singularity index:

[0153]

[0154] By definition, when the index an of a cross-field singularity is -1, it corresponds to an irregular vertex with index k = -1 in a quadrilateral mesh, i.e., valence five.

[0155] In the discrete setting, an oriented edge single-loop (ei1,..., ein) around vertex i is extracted and the sum of the angular differences is computed. In practice, there are three cases: the sum is zero and the vertex is regular, or the sum is one or negative one and the vertex is a singularity of the frame field. Due to the Crouzeix-Raviart discretization of the frame field (one angle per edge), cross-field singularities may lie on vertices, edges, or triangles of the mesh. When a singularity lies on an edge, the indices of its two adjacent vertices are both non-zero, and when a singularity lies on a triangle, the indices of its three adjacent vertices are all non-zero.

[0156] As Figure 3 , for the frame field and singularity information calculated from the background mesh in Figure 2-1 and the geometric feature constraints in 2-2.

[0157] Step 3: Implement the Runge-Kutta method at the singular points and the calculated partial boundary points to extend the streamlines;

[0158] 3.1. First, calculate the first discrete point of the starting point in a one-domain and a two-domain at the above points:

[0159] ① First, calculate the first discrete point in the relevant domain of the boundary singularity

[0160] ② Calculate whether there are internal singular points at the boundary singular points. If so, execute ③; otherwise, jump to ④;

[0161] ③ If the singular point is within the surface, extend this discrete point to the edge / point of the triangle; otherwise, directly use it as a discrete point, and then delete this singular point;

[0162] ④ Screen all the discrete points obtained above that are related to the boundary singular points, and remove the discrete points with close directions

[0163] ⑤ Calculate the first discrete point in the relevant domain of the remaining internal singular points

[0164] ⑥ If there are other singular points in the two domains of a singular point, if so, execute ⑦; otherwise, jump to ⑧,

[0165] ⑦ Select a singular point, and use the other singular point as a discrete point of this singular point. If this discrete point is within the surface, extend it to the edge / point of the triangle; otherwise, directly use it as a discrete point, and then delete the other singular point;

[0166] ⑧ Screen all the discrete points obtained above that are related to the internal singular points, and remove the discrete points with close directions

[0167] 3.2. Extend the streamline of the obtained singular points and the corresponding first discrete points according to the Runge-Kutta method. The termination conditions are as follows:

[0168] ① Until the one-domain of the singular point, terminate with the singular point

[0169] ② Terminate with the nearest triangular grid point at the termination of the constraint edge / boundary edge;

[0170] ③ The same streamline passes through the same background grid triangle in the same direction, and terminate;

[0171] Figure 4 is a schematic diagram of the initial streamline.

[0172] Step 4. Simplify the streamline structure and construct a topological partition

[0173] 4.1 For the streamline structure obtained in Step 3, detect its streamline situation, and perform different simplifications for different situations until no simplifiable situation appears. All simplification situations are roughly divided into three cases:

[0174] ① The starting and ending points of two streamlines are the same, denoted as A and B respectively, and the lengths differ greatly. Modify the longer streamline to A->A, and keep the shorter streamline unchanged

[0175] ② The starting and ending points of two streamlines are the same, and the lengths are not much different. Remove one of the streamlines

[0176] ③ If one of the starting and ending points of two streamlines is the same, denoted as A, and the additional points are B and C respectively. Without loss of generality, assume A->B is longer than A->C, and change the long streamline to B->C;

[0177] 4.2 Find the intersection points between streamlines and construct the topological relationship.

[0178] Figure 5 It is a schematic diagram of the initial topological partition.

[0179] Step 5, Topological simplification of chord folding. Based on the dual structure of quadrilateral meshes and keeping the number of singular points and the basic topological relationship in the mesh unchanged, a topological simplification operation and optimization framework for quadrilateral meshes is designed. The algorithm has high execution efficiency and good simplification effect, and can be applied to the simplification of various quadrilateral meshes.

[0180] The definitions of chord, horizontal edge, and vertical edge are as follows:

[0181] Chord: In the topological partition of a quadrilateral region, the quadrilateral region formed by the set of opposite sides corresponding to each side.

[0182] Horizontal edge: All the opposite sides from the boundary to the opposite side of another boundary within the chord.

[0183] Vertical edge: The remaining edges within the chord except the horizontal edges (the edges connecting the vertices of the horizontal edges).

[0184] Chord folding operation:

[0185] Delete an entire chord from the dual representation of the model, which means deleting all the quadrilateral patches passed by this chord. The specific deletion operation is completed by folding all the horizontal edges passed by the chord, that is, merging the vertices of the chord horizontal edges pairwise, and merging the regional lines on the corresponding vertical edges pairwise. The dual chord folding in the quadrilateral mesh does not affect the connectivity of the mesh, and it can ensure that the folded mesh is still a pure quadrilateral mesh.

[0186] Folding operation of quadrilateral region:

[0187] The folding of the quadrilateral region is a topological region deletion operation achieved by merging the two opposite vertices of the multi-block regions formed by combining some quadrilateral regions within the chord. Compared with the folding of the dual chord, this is a local mesh simplification operation. The quadrilateral region can be imagined as two triangles, and the edge between the vertices to be merged is the common edge, and the folding of the quadrilateral can be regarded as the folding of the edge of the triangle. However, this common edge is linearly interpolated and synthesized by the two vertical edges of the chord. Since the vertical edge is not a simple straight line, but is composed of multiple regional lines and the number of points on the edge is inconsistent, the interpolation process also needs to fit the same number of points for each regional line. The corresponding calculation formula is as follows:

[0188] qi =(1 - t i )×a i +t i ×b i , q i+1 =(1 - t i+1 )×a i+1 +t i+1 ×b i+1

[0189]

[0190]

[0191] Where i represents the serial number of the regional line, a i a i+1 , b i b i+1 are the regional lines on the longitudinal side within the chord, a i b i are the transverse sides within the chord. The regional line is composed of many broken lines. j represents the serial number of the broken line on the regional line, t i is x i,j the point on the longitudinal side regional line

[0192] The influence range of the folding operation of the quadrilateral region is only the region connected to the folded quadrilateral, and the whole operation satisfies the restriction of not introducing non - quadrilateral regions.

[0193] Priority metric:

[0194] To improve the connectivity of the model and ensure its geometric fidelity, it is particularly important to use an algorithm that can intelligently select elements to be deleted. The priority of the folding operation is sorted by calculating the energy formula for the chords in the original topological region, and as many long and irregular quadrilateral regions as possible are deleted to simplify a quadrilateral region topological dissection structure into an optimal topological structure with the fewest singular points. The energy formula is as follows:

[0195]

[0196] Where K controls the chord energy ratio, which is set to 4 in this algorithm example.

[0197] The specific steps of the algorithm are as follows:

[0198] ①: In the topological division of the quadrilateral region, divide the chord according to the opposite sides (search for the opposite side from the boundary side until the other boundary side or start from this side and go around a circle back to this side)

[0199] ②: Classify the chords by type:

[0200] (1) There are no singular points or the singular points lie on the same vertical edge

[0201] (2) There are two singular points on the same horizontal edge

[0202] (3) Singular points exist on both vertical edges and they are not on the same horizontal edge

[0203] ③: Check the three types of chords:

[0204] (1) If the merged vertical edge is a boundary edge or a constraint edge, discard the chord.

[0205] (2) Chords of the first type do not need to be processed. Chords of the second and third types are transferred to step four.

[0206] ④: Sort the selected chords according to the energy formula of the priority metric, and select the chord with the maximum energy greater than 0 for topological region simplification. If there is no chord with energy greater than 0, end the simplification.

[0207] ⑤: Check the type of the chord for region simplification:

[0208] (1) For chords of the second type, directly perform the folding operation on the quadrilateral region

[0209] (2) Chords of the third type will be divided into non - zip - patch regions and zip - patch regions according to the positions of the singular points (zip - patch is the region between adjacent opposite - corner singular points)

[0210] In the non - zip - patch region, delete the inner opposite - side vertical edge and points of the chord closest to its singular point, and then merge it with the quadrilateral region outside the vertical edge. In the zip - patch region, linearly interpolate and fit a new vertical edge according to the original horizontal and vertical edges, delete the original two vertical edges, and then divide the horizontal edge into two according to the new points and merge it with the outside of the chord.

[0211] (3) Return to ①

[0212] Figure 6 It is the schematic diagram of the simplified topological partition.

[0213] Step 6: For the above - obtained region division, calculate the number of segments of each region's edge based on the background size length, and perform a certain optimization on this number of segments:

[0214] Since generating a mesh using the template method depends on the segment number settings for each edge, regional line quantization is important. First, mesh the model regional lines, and then generate templates for the model regional faces. Limited by curve discretization, this works well for generating triangular (or tetrahedral) meshes, but for quadrilateral (or hexahedral) meshes, the quadrilateral topological constraints (strings) are global and pass through the model curves. For example, imagine a simple rectangular region. According to the continuous size field, the ideal number of edges for its four sides is [9.4, 4.1, 10.6, 3.9]. Simple quantization by integer rounding would separately select [9, 4, 11, 4] edges. Although there are quadrilateral meshes with such quantization, they are not structured because they must include a pair of +1 / –1 irregular vertices to transition from 9 edges on one side to 11 edges on the other side.

[0215] Quantization of the regional lines is very simple. Use the integer-rounding value of the ideal number of edges calculated by integrating the global size field, and impose equality on opposite sides, except for the regional lines of quadrilateral regions. This non-optimal quantization results in many necessary dipoles (a pair of +1 / –1 irregular vertices) in the final quadrilateral meshes of different regional faces.

[0216] Quantization calculation of the number of segments on regional lines

[0217] Consider a regional line C, which is parameterized by t varying from a to b. The floating-point number of edges on the curve according to the size map is: where dl = ||x'(t)||dt. To obtain an integer number of edges n e , use a simple rounding:

[0218] The parameter t associated with the i-th interior point of the regional line i , i ∈ [1, n e - 1] is such that:

[0219]

[0220] This integral equation can be solved by numerical integration, by adding values along the regional line until the sum equals and performing linear interpolation between samples. To mesh the regional line C, calculate the parameter t using (Equation 11) i and evaluate the regional line parameterization: x i = f(t i ) to calculate the vertex positions x i , i ∈ [1, n e - 1]. Using this method, the points are well placed on the curve according to the size map, smoothly transitioning from smaller feature regions to coarser regions.

[0221] Topological Constraints of Quadrilateral Regions

[0222] The integrated size field does not consider the particularity of the topology of quadrilateral meshes, so that quadrilaterals are organized into topological chords (double loops of adjacent elements). In this work, the number of sides n on the opposite sides of the quadrilateral region face is adjusted e . Consider a region face whose boundary consists of four region lines (C1, C2, C3, C4), and force the number of points on the opposite sides to be equal, i.e., n e1 = n e3 and n e2 = n e4 , unless the integral values are very different. When two adjacent quadrilateral region faces share a region line C c , the n ec values on the two faces must be the same. This means that the equality constraint is propagating. To solve the propagation problem, topological chords related to the quadrilateral region face are constructed. The chords consist of topologically parallel region lines, and they all receive the same fixed number of points, which are calculated by averaging the ideal values previously calculated on each curve. Through this simple propagation, the chord only propagates when two adjacent quadrilateral faces share a curve, which is one of their four sides

[0223] Step 7. Generate meshes for each partition using the template method

[0224] Steps

[0225] 1) Input the corner coordinates of the polygon region and the subdivision numbers of the corresponding sides in clockwise / counterclockwise order

[0226] 2) Select a suitable topological template according to the input subdivision numbers, and determine the topological structure of the patch by calculating the corresponding parameters

[0227] 3) Calculate the coordinates of the points on the boundary, and use Laplacian smoothing to determine the geometric position of the patch

[0228]

[0229] subject to v i = ω i , i ∈ C

[0230] where ε is the edge set of the mesh, C is the boundary point set, and ω i ∈ R 2 is the input coordinate of the i-th fixed boundary point

[0231] Simplify the input

[0232] Since considering all possible input cases is a challenge, the problem can be considered simplified to an equivalent and easier-to-solve sub-problem. This operation greatly reduces the cases to be considered, thus ensuring the generality of the algorithm.

[0233] For the input (l0,..., l N-1 ), for some k, l k-1 and l k+1 are both greater than 1. Define d = min(l k-1 , l k+1 ) - 1. Then the input can be simplified to (l0,..., l N-1 ), where

[0234]

[0235] The above simplification operation can be repeated until no further simplification is possible, thus obtaining the maximally simplified input.

[0236] The template selection is formulated as an ILP

[0237] Each topological template defines a linear relationship between the input (l0,..., l N-1 ) and the parameters, which include the padding amount p i of the i-th edge, the number of edge flows x and y.

[0238] The above general form is denoted as:

[0239] Ax = b

[0240] where A is an N×M matrix, M is the number of parameters of the template, x represents the M-dimensional vector of template parameters (to be solved), b is an N-dimensional vector jointly determined by the input (l0,..., l N-1 ) and the template, and it is required that is a non-negative integer. It can be expressed as an integer linear programming (ILP):

[0241] argmax c T x

[0242] subject to Ax = b

[0243] x ≥ 0, x ∈ Z M

[0244] Among them, c represents an M-dimensional vector of the target to be maximized. As long as the ILP has a solution, it means that the template is feasible for the input. Since M, that is, the number of parameters of a certain template, is very small (at most 10), it can be solved quickly. Since N ≤ M, there are usually multiple solutions to the above ILP, which means that for the same topological template, there are multiple parameter combinations that meet the input requirements. The goal is set to maximize Under this goal, when the corresponding variable represents the boundary filling amount, c j = 1; otherwise, c j = 0, and the number of singular points included in the topological generation at this time is the least.

[0245] Given the region division result of the problem domain and the number of subdivisions of the corresponding sides of each sub-region, the quadrilateral mesh generation of each sub-region is obtained using the above algorithm, and then they are stitched together to form the final mesh effect.

[0246] Figure 7 It is a schematic diagram of the result of generating a mesh by the template method according to the partition points calculated in step 7.

[0247] Step 8, Mesh smoothing

[0248] ① Regular vertices: Given a fixed boundary, solve the Winslow nonlinear elliptic PDE Δ x u = 0, where u(x) is a certain coordinate calculation space and x is the coordinate in the physical space. The advantage of doing this is that the two coordinate components are coupled, and the resulting quadrilateral shape is good, and some orthogonality can be enforced even under large distortions. By applying finite difference (FDM) discretization to the Winslow equation, the kernel of regular vertices in the locally smoothed quadrilateral mesh can be derived. Assume that (x1,..., x8) are the vertices of the template ordered around the regular vertex x, and its new position is given by:

[0249]

[0250] where α0 = (x3 - x7)·(x3 - x7), α1 = (x1 - x5)·(x1 - x5), β = (x1 - x5)·(x3 - x7)

[0251] ② Irregular vertices: Use an angle-based smoothing method. The idea is to move the vertices on the single-ring bisector. As the above Winslow kernel, it works well when the mesh is not very restricted, but may create flipped cells in complex configurations.

[0252] Figure 8 It is a schematic diagram of the final mesh for geometric optimization of the mesh result in step 7.

[0253] In summary, by combining the extraction of the simple topology of the frame field with the generation of the mesh from the template, the present invention obtains a quadrilateral mesh that meets the requirements of high precision, high quality, and effectiveness for complex geometric feature constraints and dimensional constraints, and thus can be applied to the meshing research and design of high-level CAD and CAE models.

Claims

1. A method for generating complex multi-constrained structured grids based on a frame field, characterized in that Including the following steps: Step 1: Convert the input background triangular mesh and the corresponding dimensional constraints in the dimensional constraints and geometric feature constraints into corresponding geometric constraints; 1.1 Use the marching triangle method to extract the isoclines for each triangular face; 1.2 Connect the extracted isoclines and perform a certain degree of smoothing; 1.3 Filter the isoclines with a small number; 1.4 Re-mesh at the isoclines and store the isoclines as geometric feature constraints; Step 2: Use the finite element method to calculate the frame field of the background mesh and extract the singular points therein; Step 3: Implement the Runge-Kutta method at the singular points and some of the calculated boundary points to extend the streamlines; 3.1 Calculate the first discrete point of the starting point in the one-neighborhood and two-neighborhood at the singular points and some of the calculated boundary points: ① Calculate the first discrete point in the relevant neighborhood of the boundary singular point; ② Calculate whether there are internal singular points at the boundary singular point. If so, execute ③, otherwise jump to ④; ③ If the singular point is inside the face, extend this discrete point to the side / point of the triangle, otherwise directly use it as a discrete point, and then delete this singular point; ④ Screen all the discrete points related to the boundary singular point obtained above and remove the discrete points with close directions; ⑤ Calculate the first discrete point in the relevant neighborhood of the remaining internal singular points; ⑥ If there are other singular points in the two-neighborhood of the current singular point, execute Step ⑦; Otherwise, execute Step ⑧; ⑦ Select a singular point and use the other singular point as a discrete point of this singular point. If this discrete point is inside the face, extend it to the side / point of the triangle, otherwise directly use it as a discrete point, and then delete the other singular point; ⑧ Screen all the discrete points related to the internal singular point obtained above and remove the discrete points with close directions; 3.2 Extend the streamlines of the obtained singular points and the corresponding first discrete points according to the Runge-Kutta method. The termination conditions are as follows: ① Until the one-neighborhood of the singular point, terminate with the singular point; ② Terminate with the nearest triangular mesh point at the termination of the constrained edge / boundary edge; ③ The same streamline passes through the same background mesh triangle in the same direction, terminate; Step 4: Simplify the streamline structure and construct the initial topology; 4.1 For the streamline structure obtained in Step 3, detect its streamline situation and simplify it until there is no simplifiable situation. The simplification situations include: ① Two streamlines have the same starting and ending points, denoted as A and B respectively, and the lengths differ greatly. Change the longer streamline to A->A, and the shorter streamline remains unchanged; ② Two streamlines have the same starting and ending points and the lengths are not much different. Remove one of the streamlines; ③ Two streamlines have one same starting and ending point denoted as A, and the additional points are B and C respectively. Without loss of generality, assume A->B is longer than A->C. Change the longer streamline to B->C; 4.2 Find the intersection points between the streamlines and construct the topological relationship; Step 5: Topological simplification by chordal collapse. Based on the quadrilateral mesh dual structure and keeping the number of singular points and the basic topological relationship in the mesh unchanged, design a quadrilateral mesh topological simplification operation and optimization framework, which is applicable to the simplification of various quadrilateral meshes; The definitions of chord, horizontal side, and vertical side are as follows: Chord: In the topological partitioning of a quadrilateral region, a quadrilateral region formed by the set of opposite sides corresponding to each side; Horizontal side: All opposite sides within the chord starting from one boundary to the opposite boundary; Vertical side: The remaining sides within the chord excluding the horizontal sides; Folding operation of the chord: Merge the vertices of the horizontal sides of the chord in pairs, and merge the region lines on the corresponding vertical sides in pairs; Folding operation of the quadrilateral region: The folding of the quadrilateral region is a topological region deletion operation achieved by merging the two opposite vertices of the multiple regions formed by combining some of the quadrilateral regions within the chord; imagining the quadrilateral region as two triangles, the side between the vertices to be merged is the common side, and the common side is linearly interpolated and synthesized by the two vertical sides of the chord. The corresponding calculation formula is as follows: q i =(1 - t i )×a i +t i ×b i ,q i+1 =(1 - t i+1 )×a i+1 +t i+1 ×b i+1 where i represents the serial number of the regional line, a i a i+1 , b i b i+1 are the regional lines on the longitudinal side within the chord, a i b i are the transverse sides within the chord. The regional line is composed of many broken lines. j represents the serial number of the broken line on the regional line, t i is x i,j is the point on the longitudinal side regional line; Step 5 includes the following sub-steps: ①: Divide the chord according to the opposite sides in the topological partitioning of the quadrilateral region; ②: Classify the chords by type: (1) There are no singular points or singular points exist on the same vertical side; (2) There are two singular points on the same horizontal side; (3) Singular points exist on both vertical sides and are not on the same horizontal side; ③: Check the three types of chords in ②: (1) When it is judged that the vertical side to be merged is a boundary side or a constraint side, discard the chord; (2) The first type of chord does not need to be processed, and the second and third types of chords are transferred to ④; ④: Sort the selected chords according to the energy formula of the priority metric, and select the chord with the maximum energy and greater than 0 for topological region simplification; if there is no chord with energy greater than 0, end the simplification; The priority metric includes: Use the energy formula to select and delete those long and irregular quadrilateral regions, and simplify the topological dissection structure of a quadrilateral region into an optimal topological structure with the fewest singular points; the energy formula is as follows: Where K controls the chord energy ratio; ⑤: Check the type of the chord for region simplification: (1) For the second type of chord, perform the folding operation of the quadrilateral region; (2) The third type of chord will be divided into non-zip-patch regions and zip-patch regions according to the positions of the singular points; For the non-zip-patch region, delete the in-chord opposite side vertical side and points of the singular point closest to it, and then merge with the quadrilateral region outside the vertical side; For the zip-patch region, linearly interpolate and fit a new vertical side according to the original horizontal side and vertical side, delete the original two vertical sides, and then divide the horizontal side into two parts according to the new points and merge with the outside of the chord; (3) Return to ①; Step 6: For the region partitioning obtained above, calculate the number of segments of each region edge based on the background size length and optimize this number of segments: Step 6 includes the following sub-steps: Sub-step 1, use the quantization calculation method to calculate the integer rounded value of the ideal number of edges for all region lines in combination with the size field; the quantization calculation method for the number of segments on the region line: Let there be a region line C, which is parameterized by t varying from a to b; the number of floating-point edges on the curve according to the dimension drawing is: where dl = ||x'(t)||dt; to obtain an integer number of edges n e , a simple rounding is used: Sub-step 2, use the topological constraints of the quadrilateral region quantization to impose equality on the number of segments of some quadrilateral regions; The topological constraints of the quadrilateral region: Adjust the number of sides n on opposite sides of the quadrilateral region surface e ; Calculate these points by averaging the ideal values previously calculated on each curve; Through this simple propagation, the chord propagates only when two adjacent quadrilateral faces share a curve, which is one of their four sides; Sub-step 3, perform secondary quantization on all region lines with the final ideal number of edges obtained; The parameter t associated with the i-th interior point of the region line i , where i ∈ [1, n e -1] is as follows: By adding values along the regional line until the sum equals Perform linear interpolation between samples; to mesh the regional line C, calculate the parameter t using Equation 11 i And evaluate the regional line parameterization: x i = f(t i ) to calculate the vertex positions x i , i ∈ [1, n e - 1]; using this method, the points are well placed on the curve according to the dimension drawing, smoothly transitioning from the smaller feature area to the coarser area; Step 7: Use the template method to generate grids for each partition; Step: 1) Input the corner coordinates of the polygon area and the subdivision number of the corresponding sides in a clockwise / counterclockwise manner; 2) Select a suitable topological template according to the input subdivision number, and determine the topological structure of the patch by calculating the corresponding parameters; 3) Calculate the coordinates of the points on the boundary, and use Laplacian smoothing to determine the geometric position of the patch; the formula for Laplacian smoothing is as follows: subject to v i = ω i , i ∈ C where ε is the edge set of the grid, C is the boundary point set, and ω i ∈R 2 is the coordinate of the i-th fixed boundary point as the input; Step 8, Mesh Smoothing ① Regular vertex: Given a fixed boundary, solve the Winslow non-linear elliptic PDE △ x u = 0, where u(x) is a certain coordinate calculation space and x is the coordinate in the physical space; by applying finite difference discretization to the Winslow equation, the kernel of the regular vertex in the locally smoothed quadrilateral mesh is derived; let (x1,..., x8) be the vertices of the template that orderly surround the regular vertex x, and its new position is given by: where α0 = (x3 - x7)·(x3 - x7), α1 = (x1 - x5)·(x1 - x5), β = (x1 - x5)·(x3 - x7) ② Irregular vertices: Use an angle-based smoothing method, the idea of which is to move the vertices on the single-ring bisector; as the Winslow kernel above.

Images