Oral cavity finite element simulation calculation method based on collision detection

Abstract

The invention discloses an oral cavity finite element simulation calculation method based on collision detection. The method comprises the following steps: 1) acquiring oral scanning data of teeth, generating surface grid data of the invisible appliance according to the oral scanning data, and then generating a tetrahedral grid from the surface grid of the invisible appliance for subsequent simulation calculation; 2) on the basis of collision detection, calculating a pushing vector of teeth to the invisible appliance, and taking the pushing vector as a displacement boundary condition of subsequent calculation; 3) constructing a three-dimensional finite element model taking voxels as units, calculating an element stiffness matrix, assembling an overall stiffness matrix, loading displacement boundary conditions, forming a finite element solution equation set, and obtaining a model node displacement result; and 4) calculating stress and strain of the finite element model according to a displacement result. Aiming at the orthodontic mechanical analysis of the invisible orthodontic appliance, the invention provides a modeling method for the interaction between the teeth and the invisible orthodontic appliance based on a collision detection technology, completes the subsequent simulation solution of a finite element model, and gives consideration to the solution precision and speed.

Description

Technical Field

[0001] This invention belongs to the field of computer-aided engineering, specifically relating to a collision detection-based finite element simulation calculation method for the oral cavity, which establishes a three-dimensional finite element model based on oral cavity scanning data and realizes the analysis and calculation process. Background Technology

[0002] The performance analysis of clear aligners primarily focuses on their mechanical properties. Currently, the most commonly used in vitro research method in the field of orthodontics is three-dimensional finite element analysis (FEM). In modern engineering and scientific research, the finite element method (FEM) is widely used to simulate the mechanical properties of various structures and materials. This method, by decomposing complex structures into smaller, manageable units (or finite elements), can simulate and analyze the response of materials under various forces while ensuring high simulation accuracy.

[0003] Currently, most finite element analyses of clear aligners rely on complex commercial software. Their boundary condition application methods typically involve applying displacement to the portion of the clear aligner corresponding to the treated tooth, calculating the stress distribution within the aligner, or directly assuming the aligner provides ideal orthodontic force to the target tooth and then analyzing the stress distribution in the jawbone model. However, in reality, the interaction between the clear aligner and the teeth is a holistic process. The constraint of the teeth on the clear aligner should not be limited to the treated tooth; untreated teeth also provide support for the aligner. Furthermore, the clear aligner exerts force not only on the target tooth but also on other teeth; the influence on other teeth should be considered. Therefore, when determining the boundary conditions for clear aligners, a contact analysis of the entire clear aligner and the teeth should be performed.

[0004] Contact problems between deformable bodies exist in practical engineering. The application of invisible aligners to teeth for orthodontic treatment is a typical example of contact mechanics problems. These problems typically involve triple nonlinearity in geometry, materials, and boundaries, making them among the most challenging nonlinear problems. Numerical methods are the most effective means of solving such problems. Extensive analysis shows that the contact algorithm is a crucial factor affecting the accuracy of numerical calculations. Furthermore, determining the contact interface and contact state during the contact process is extremely time-consuming; contact calculations typically account for more than half of the total solution time. Therefore, developing efficient and high-precision contact algorithms is urgently needed for practical engineering applications. Summary of the Invention

[0005] The present invention provides a finite element simulation calculation method for oral cavity based on collision detection, which solves the problems mentioned in the background art.

[0006] This invention considers the interaction area between teeth and braces holistically, using the braces collision depth as a displacement boundary condition in the simulation calculation to model the interaction process between the braces and teeth in real-world differential regions. The model is then meshed and subjected to finite element analysis (FEM). Dedicated storage structures are used for memory optimization of computational resources, and specialized data structures are employed to optimize the computational speed. Finally, based on the displacement results obtained from the calculation, stress and strain are calculated for reference by orthodontists when designing treatment plans.

[0007] This invention provides a finite element simulation calculation method for the oral cavity based on collision detection, comprising the following steps:

[0008] Step 1: Scan the dental model image data, apply the doctor's treatment plan to the dental model, and obtain the desired post-treatment dental model. Import the desired post-treatment dental model into Geomagic to obtain the surface mesh model of the invisible aligner.

[0009] Dental model image data, i.e., the surface mesh model of the teeth, is obtained through oral scanning.

[0010] Step 2: Process the desired post-treatment tooth model and the clear aligner model using a solver to obtain the stress distribution of the clear aligner. The specific process is as follows:

[0011] 2.1) Use the Delaunay algorithm to generate a tetrahedral mesh model, and set the corresponding physical properties according to the actual invisible aligner material used;

[0012] 2.2) Extract the surface nodes of the invisible aligner mesh, calculate the distance from the node inside the tooth model to the tooth model (i.e., the intrusion distance), and assign the intrusion vector as a displacement boundary condition to the node. All data is stored in a dedicated data structure for subsequent simulation calculations.

[0013] 2.3) Based on the previous data, the stiffness matrix of the tetrahedral elements is calculated using multithreading, and then assembled into the global stiffness matrix. Displacement boundary conditions are applied to construct the final finite element solution equations. The equations are calculated using the LU decomposition method to obtain the solution vector of the equations, which is the displacement vector of the finite element model.

[0014] 2.4) Calculate the stress and strain of the finite element model based on the displacement vector and derive the results.

[0015] Step 1 specifically involves:

[0016] First, apply the orthodontic treatment plan provided by the doctor to the dental model obtained by oral scanning, then import it into Geomagic, and offset it outward (towards the direction of the tooth surface normal vector) by 2mm (this offset value corresponds to the actual thickness of the manufactured clear aligner). Then use the shelling function to obtain the clear aligner model.

[0017] In step 2.2), the specific process of extracting the surface nodes of the invisible orthodontic appliance mesh is as follows:

[0018] After the tetrahedral mesh is generated, each geometric element will have a corresponding data structure for storage, including point elements, edge elements, triangular facet elements, and tetrahedral elements. Triangular facet elements are essentially elements composed of nodes on the surface of the mesh model. By traversing all triangular facet elements, the corresponding node ID is obtained and stored in a `std::set` structure. This fully utilizes the unique storage characteristic of `std::set` elements to achieve automatic deduplication.

[0019] The intrusion calculations for each node and the tooth model are independent of each other, and this step can be accelerated using multithreading. Multithreading traverses all surface nodes, first using ray casting to determine if a node is inside the model. An odd number of intersections between the ray and the model indicates an interior node, while an even number indicates an exterior node. For nodes inside, the nearest point to the tooth surface is found, and this vector is assigned to the node as a boundary condition for subsequent finite element calculations. This step simulates the process of teeth pushing and deforming braces when actually wearing them.

[0020] Calculating the nearest point on the tooth surface involves using a kd-tree and an r-tree for faster lookup. The kd-tree is used to divide the tooth into nodes that make up the tooth, and the r-tree is used to divide it into triangular facets. First, the kd-tree is searched to find the tooth node closest to the surface node. Then, a region is drawn with these two points as radii, and the r-tree is searched to obtain all facets within this region. The nearest point of the surface node among all facets is calculated; the final nearest point is the closest point from the surface node to the tooth model.

[0021] Step 2.3) specifically refers to:

[0022] 1. The integral formula for calculating the element stiffness matrix includes material matrix and gradient matrix information; the material matrix is ​​obtained based on the input material properties, and the gradient matrix is ​​obtained from finite element theory; the Gaussian quadrature rule is used in the calculation, and multi-threaded parallel computation is employed.

[0023] First, the maximum number of concurrent connections supported by the hardware is obtained. Then, each thread is assigned a task to calculate the element stiffness matrix. Each thread calculates the material matrix based on the input material properties.

[0024] First, the maximum number of concurrent connections supported by the hardware is obtained. Then, each thread is assigned a task to calculate the element stiffness matrix. Each thread calculates the material matrix based on the input material properties and then obtains the gradient matrix based on the finite element theory. For each element, the Gaussian quadrature method is used to calculate the element stiffness matrix based on the material matrix and the gradient matrix, thus obtaining a 12*12 element stiffness matrix, which is stored using triples.

[0025] The tetrahedral elements are mapped back to standard tetrahedrons by calculating the Jacobian matrix. The element stiffness matrix is ​​calculated using the Gaussian quadrature method. Specifically, the four vertices of the standard tetrahedron are used as four sampling points. The calculation results of the four sampling points are accumulated to obtain the integral result. Finally, the integral result is multiplied by the mapped Jacobian matrix to obtain the final 12*12 element stiffness matrix, which is stored using triples.

[0026] Finally, all element stiffness matrices are assembled to obtain the final global stiffness matrix. The global stiffness matrix uses the classic COO sparse matrix storage format, which automatically accumulates the values ​​at corresponding positions when generating from the set of triples, achieving the purpose of fast assembly.

[0027] 2. For displacement constraints, we must ensure that the corresponding nodal results after solving the equations are the values ​​we preset. We use the method of multiplying the diagonal elements by a large number, multiplying the corresponding positions of the overall stiffness matrix and the corresponding positions of the right-hand side terms by a maximum value. This ensures that the solution result at the current position is not affected by other parts of the equation system and remains constant as our preset value.

[0028] 3. The lu decomposition method for solving linear equations is as follows: For solving a large-scale linear system of equations Ax = b, directly calculating the inverse matrix of A is prohibitively expensive. However, A can be decomposed using Gaussian transformation into a unit lower triangular matrix L and a unit upper triangular matrix U. The computational cost of calculating the inverse matrices of L and U is very small. Therefore, solving the original system of equations becomes: 1. Substituting backwards to solve Ly = b, 2. Substituting backwards to solve Ux = y.

[0029] Step 2.4) specifically refers to:

[0030] First, transform the tetrahedral element into a standard tetrahedron and calculate the corresponding Jacobian matrix. Then, calculate the strain matrix of the four vertices in the tetrahedral element, multiply it by the displacement vector on the left to obtain the vertex strain tensor; then multiply the strain tensor by the material matrix on the left to obtain the vertex stress tensor, and calculate the vertex principal stress from the stress tensor.

[0031] The beneficial effects of this invention are:

[0032] 1) This invention uses Geomagic software to obtain the invisible orthodontic appliance model, and then divides it into tetrahedral meshes that better fit the irregular model for finite element simulation calculation, making the results closer to reality.

[0033] 2) This invention sets the boundary conditions for solving based on the collision depth, which more realistically simulates the interaction between the orthodontic appliance and the teeth during the orthodontic process. At the same time, it uses a specific data structure to accelerate the solution process, improving efficiency while ensuring high accuracy of the results. Attached Figure Description

[0034] Figure 1 The flowchart illustrates the finite element simulation calculation method for the oral cavity based on collision detection, as provided in this embodiment of the invention.

[0035] Figure 2 For the surface mesh model of teeth

[0036] Figure 3 Surface mesh model of invisible aligners

[0037] Figure 4 A schematic diagram of the principle of LU decomposition.

[0038] Figure 5 A visualization diagram of the model solution results provided for the embodiment. Detailed Implementation

[0039] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not limit the scope of protection of the invention.

[0040] like Figure 1 This invention provides a finite element simulation calculation method for oral cavity based on collision detection, including the following steps:

[0041] 1. Based on the oral cavity scan data, certain repairs, restorations, and fillings are performed to obtain the following results: Figure 2 The image shows a surface mesh model of the tooth. The dentist's treatment plan is applied to the tooth model to obtain the desired post-treatment tooth model. The desired post-treatment tooth model is imported into Geomagic, the outer shell is offset outwards by 1.3mm, and then subtracted from the original model to obtain the result shown below. Figure 3 The surface mesh model of the invisible orthodontic device is shown.

[0042] 2. Read in the tooth and clear aligner model obtained in step 1, divide it into tetrahedral meshes, construct a three-dimensional finite element model, calculate the element stiffness matrix, and assemble the overall stiffness matrix.

[0043] Specifically, a tetrahedral element finite element model is first constructed: two model files and material physical property data are sequentially read in. Physical properties include Young's modulus and Poisson's ratio, stored in the Material class. Both model files are triangular meshes. The braces model needs to be further subdivided into tetrahedral meshes for finite element calculations, while the tooth model serves as a boundary, provided for collision depth calculations in the braces model, and does not require further subdivision. A tetrahedral mesh model is generated based on the Delaunay algorithm. The generated model contains three parts of data: node data is stored in a structure of size NodeNumber*3, containing the node positions; tetrahedral element data is stored in a structure of size TetNumber*4, containing the node IDs corresponding to the tetrahedrons; and triangular elements are stored in a structure of size TriNumber*3, containing the node IDs corresponding to the triangular elements forming the mesh surface.

[0044] The element stiffness matrix is ​​then calculated using multi-threaded parallel computation. The maximum concurrency supported by the current device is obtained as the number of threads, and the task of calculating the stiffness matrix is ​​then grouped and assigned to each thread. Since the tetrahedrons have varying shapes, we utilize the concept of isoparametric elements, transforming the tetrahedrons to a local coordinate system for calculating the element stiffness matrix. The specific approach is as follows: The integral formula for the element stiffness matrix... The calculation of matrix B involves the partial derivatives of the shape functions, derived from the geometric equations of the tetrahedron. We calculate the Jacobian matrix of the isoparametric element corresponding to the current element, and then, based on the expression of the inverse of the Jacobian matrix and the partial derivatives of the shape function N with respect to the local coordinate system, we obtain the final partial derivatives of the shape functions and assemble them into matrix B. Based on the material's physical properties, we calculate the material matrix D. Using the Gaussian quadrature rule, each element is calculated using four sampling points, i.e., the four vertices of the tetrahedral element, yielding the integral result, which is a 12*12 element stiffness matrix.

[0045] Next, the global stiffness matrix is ​​assembled: Since the global stiffness matrix is ​​a large-scale sparse matrix, we store it using the COO format to save space. The calculated element stiffness matrix data is stored using triples, and the global stiffness matrix is ​​constructed from the set of triples, with terms in the same position added together.

[0046] 3. Perform collision calculations on the braces mesh and the tooth model, and set the boundary conditions for the finite element model. Specifically, calculate the intrusion depth of nodes on the braces surface into the teeth. The calculation is divided into two steps: first, search for nodes that collide, and then calculate the nearest point from each collided node to the tooth mesh surface.

[0047] Determining whether a node has collided essentially means determining whether the node is inside the tooth mesh. This is done using a raycasting method: radiate a ray outward from the current node, calculate the intersection points of the ray and the tooth mesh. If the number of intersection points is odd, the node is inside the mesh; otherwise, it is outside.

[0048] For a collision node, find its nearest point on the tooth surface. This is accelerated using kd-trees and r-trees. Obtain the tooth mesh nodes and construct a kd-tree, where each tree node corresponds to a mesh node; obtain the tooth surface triangles and construct an r-tree, where each tree node corresponds to a triangle's bounding box. First, search the kd-tree to find the tooth mesh node closest to the collision node. Then, using the distance between these two points as the radius and the tooth mesh node as the center, draw a rectangular region and search the r-tree to extract the triangular faces within this rectangular region. These triangular faces are the faces containing potential nearest points. The collision node sequentially searches for its nearest point on these triangular faces, and the final closest point is the collision node's nearest point on the tooth surface. Calculate a displacement vector for the collision node and the nearest point and assign it to the collision node.

[0049] After the above steps, each collision node has a corresponding displacement vector, which we need to set as a boundary condition in the final linear equation system.

[0050] There are several methods for setting displacement boundary conditions, such as direct substitution and setting diagonal elements to 1. We use the method of multiplying diagonal elements by a large number. This method does not modify the order or structure of the equation system. Each condition only requires changing the corresponding term in the global stiffness matrix and the corresponding term in the right-hand side. It is applicable to both zero and non-zero values. The specific operation is as follows: When there is a given nodal displacement, multiply the main diagonal element of that row of the equation by a large number α (on the order of e10), and change the value in the corresponding right-hand side to the product of α and the original main diagonal element and the given displacement value. After making the above correction for all given displacement values, solving the corrected equation will yield the displacement values ​​of all nodes.

[0051] 4. For example Figure 4 As shown, the LU decomposition method is used to solve the system of equations, and the final solution displacement is obtained.

[0052] Solving the linear equation system Ax = b typically involves finding the inverse of A. However, for finite element analysis, the coefficient matrix is ​​usually a large-scale sparse matrix, making the cost of finding its inverse prohibitive. LU decomposition refers to the ability of matrix A to be decomposed into a product of LUs, where L is a unit lower triangular matrix and U is a unit upper triangular matrix, such as... Figure 4 As shown. The calculation formula can be derived as follows:

[0053] u kj =akj -(l kj u 1j +···+l k，k-1 u k-1，j )

[0054] l ik =(a ik -l i1 u 1k -···-l i，k-1 u k-1，k ) / u kk

[0055] After decomposition, the original equation is transformed into LUx = b. Solving this equation involves two steps: forward substitution to solve for Ly = b, and backward substitution to solve for Ux = y. The final solution is the displacement of the model, such as... Figure 5 As shown.

[0056] Then calculate the stress: calculate the strain matrix B of the four nodes of the element, multiply it by the displacement field to obtain the nodal strain tensor, and then multiply the strain tensor by the material matrix D to obtain the nodal stress tensor.

Claims

1. A finite element simulation calculation method for oral cavity based on collision detection, characterized in that, Includes the following steps: Step 1: Obtain a tooth model through oral scanning, apply the orthodontic treatment plan to the tooth model to obtain the desired post-treatment tooth model, and import the desired post-treatment tooth model into Geomagic to obtain the surface mesh model of the invisible aligner. Step 2: Process the desired post-treatment tooth model and the clear aligner model using a solver to obtain the stress distribution of the clear aligner. The specific process is as follows: 2.1) Use the Delaunay algorithm to generate a tetrahedral mesh model, and set the corresponding physical properties according to the actual invisible aligner material used; 2.2) Extract the surface nodes of the invisible aligner mesh, calculate the distance from the nodes inside the tooth model to the tooth model (i.e., the intrusion depth) using multi-threaded calculation, and use the intrusion depth as the displacement boundary condition. 2.3) Based on steps 2.1) and 2.2), the stiffness matrix of the tetrahedral element is calculated using multithreading, then assembled into the global stiffness matrix, and displacement boundary conditions are applied to construct the final finite element solution equation set; the equation set is calculated using the lu decomposition method to obtain the solution vector of the equation set, which is the displacement vector of the finite element model. 2.4) Calculate the stress and strain of the finite element model based on the displacement vector, and derive the results; The force exerted on the teeth by the invisible aligner was calculated based on the stress and strain of the finite element model.

2. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 1, the process of importing the tooth model into Geomagic for shelling to obtain the surface mesh model of the invisible aligner is as follows: First, apply the orthodontic treatment plan provided by the doctor to the dental model obtained by oral scanning, then import it into Geomagic, offset it outward by 2mm along the normal vector of the tooth surface, and then use the shelling function to obtain the invisible aligner model.

3. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 2.2), the specific process of extracting the surface nodes of the invisible orthodontic appliance mesh is as follows: Traverse all triangular facets in the tetrahedral mesh, obtain the corresponding node IDs, and record and store them.

4. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 2.2), the multi-threaded calculation of the intrusion depth specifically involves: Traverse all surface nodes, first use the ray method to determine whether a node is inside the model: if the number of intersections between the ray and the model is odd, it means the node is inside; if the number of intersections is even, it means the node is outside. For nodes inside the tooth, find the nearest point to the tooth surface. The vector from the node to the nearest point is the intrusion depth. Use this vector as the boundary condition related to the current internal node in the subsequent finite element calculation. Specifically, calculating the nearest point on the tooth surface involves using kdtree and rtree to accelerate the query; using kdtree to divide the nodes that make up the tooth, and using rtree to divide the triangular facets that make up the tooth. First, search the kdtree to find the tooth node that is closest to the node on the surface of the clear aligner. Then, draw a region with the node on the surface of the clear aligner and its corresponding nearest tooth node as the radius. Then, find all the faces within the drawn region by searching the rtree. Calculate the nearest point of the surface node of the invisible aligner among all the facets obtained by searching the rtree. The final nearest point is the nearest point from the surface node to the tooth model.

5. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 2.2), all data is stored in a dedicated data structure for subsequent simulation calculations. The dedicated data structure is as follows: Using std::set as the storage structure for surface nodes, we can fully utilize the uniqueness of std::set in storing elements to achieve automatic deduplication. std::map is used to record node IDs and corresponding invasion depths, thus saving storage space; The Material class is designed to store the physical properties of the invisible aligner, which can be customized by the user.

6. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 2.3), the stiffness matrix is ​​calculated as follows: First, the maximum number of concurrent connections supported by the hardware is obtained. Then, each thread is assigned a task to calculate the element stiffness matrix. Each thread calculates the material matrix based on the input material properties and then obtains the gradient matrix based on the finite element theory. For each element, the Gaussian quadrature method is used to calculate the element stiffness matrix based on the material matrix and the gradient matrix, thus obtaining a 12*12 element stiffness matrix, which is stored using triples. Finally, the stiffness matrices of all elements are assembled to obtain the final global stiffness matrix.

7. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 2.3), the application of displacement boundary conditions specifically involves: The method of multiplying diagonal elements by a large number is adopted, which multiplies the corresponding positions of the overall stiffness matrix and the corresponding positions of the right-hand side terms by a maximum value.

8. The finite element simulation calculation method for oral cavity based on collision detection according to claim 1, characterized in that, In step 4, the algorithm for calculating stress and strain is as follows: The tetrahedral elements are mapped back to standard tetrahedrons by calculating the Jacobian matrix. The strain matrix of the four vertices in the tetrahedral element is calculated by using the material matrix and gradient matrix. Multiplying the strain matrix by the displacement vector on the left yields the vertex strain tensor. Multiplying the strain tensor by the material matrix on the left yields the vertex stress tensor.

Images