Patent Publication Number: US-2019178767-A1

Title: Method To Carry Out Accurate Finite Element Analysis Over A Tangled Mesh

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is a continuation of U.S. application Ser. No. 13/828,673, filed Mar. 14, 2013. 
    
    
     FIELD OF THE INVENTION 
     This invention relates generally to finite element analysis, and in particular, to provide a method for carrying accurate finite element analysis over a tangled mesh. 
     BACKGROUND AND SUMMARY OF THE INVENTION 
     The finite element analysis (FEA) is a numerical analysis technique that allows a user to solve complex problems by breaking it down into a finite number of simple problems. FEA is often used by engineers and designers in new product design and in the refinement of existing products. Using FEA, a user is able to verify a proposed design of a product will perform in accordance with desired specifications prior to manufacturing or construction of such product. 
     In FEA, the geometry of the proposed design is identified. The geometry is broken up into a discrete representation, known as a mesh or grid. The mesh is made up of a plurality of finite elements defined by simple polynomial shape functions, wherein the vertices of the shape define nodes. Boundary conditions (e.g. stress, constraints and/or loads) are applied to the mesh and the displacement of the elements is determined by the nodal displacements. Once the nodal displacements are known, element stresses and strains can be calculated. The governing equations are assembled into matrix form and are solved numerically. 
     In classic FEA, the underlying mesh is required to satisfy certain topological and geometrical properties. One such critical property is that none of the elements must overlap, i.e., the mesh must not be tangled. Because tangled meshes generate physically invalid solutions, it is imperative that the meshes be untangled prior to solving the governing equations. However, as hereinafter described, the tangling of a finite element mesh is unavoidable in modern FEA. 
     The tangling of finite element meshes can occur during various operations. For example, in an effort to generate all-hex meshes, mesh generators can inadvertently introduce tangling. Similarly, mesh optimizers can introduce tangling. Further, in large-scale deformation excessive node movement relative to the local size of the finite elements can result in tangling. Finally, tangling can occur when an existing mesh is morphed to reflect a design change. The untangling of a mesh is perhaps as difficult as mesh generation and optimization, and has received considerable attention from the finite element community. As such, it would be highly desirable to provide an alternative to ‘untangling’ of a finite element mesh which would allow a user to generate physically valid solutions to the equations governing a proposed design. 
     Therefore, it is a primary object and feature of the present invention to provide a method for carrying accurate finite element analysis over a tangled mesh. 
     It is a further object and feature of the present invention to provide a method for carrying out accurate finite element analysis over a tangled mesh which is simple and inexpensive to implement. 
     It is a still further object and feature of the present invention to provide a method for carrying out accurate finite element analysis over a tangled mesh which allows a user to generate physically valid solutions to the equations governing a proposed design. 
     In accordance with the present invention, a method is provided for carrying accurate finite element analysis. The method includes the step of meshing a domain subjected to a field. The mesh is defined by a plurality of finite elements. A stiffness contribution due to the plurality of finite elements is computed and a stiffness contribution due to the finite elements overlapping corresponding finite elements is computed. The stiffness contribution due to the plurality of finite elements and the stiffness contribution due to the overlapping finite elements are combined. 
     Each of the finite elements is defined by a plurality of nodes and each finite element has an orientation. In addition, each node has a nodal shape function. The nodal shape functions are defined according to the expression: 
     
       
         
           
             
               
                 φ 
                 i 
               
                
               
                 ( 
                 • 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   j 
                   ∈ 
                   
                     C 
                      
                     
                       ( 
                       i 
                       ) 
                     
                   
                 
               
                
               
                 
                   Θ 
                   j 
                 
                  
                 
                   
                     N 
                     
                       i 
                       , 
                       j 
                     
                   
                    
                   
                     ( 
                     • 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein ϕ i (⋅) are the nodal shape functions; j is a finite element; i is a node; C i  is a set of finite elements connected to node i; Θ j  is the orientation of the finite element j; and N i,j (⋅) are element shape functions. 
     The stiffness contribution due to the plurality of finite elements may be calculated according to the expression: 
     
       
         
           
             
               K 
               standard 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∫ 
                   
                     E 
                     j 
                   
                 
                  
                 
                   
                     
                       ( 
                       
                         ∇ 
                         
                           N 
                           i 
                         
                       
                       ) 
                     
                     · 
                     
                       ( 
                       
                         ∇ 
                         
                           N 
                           j 
                         
                       
                       ) 
                     
                   
                    
                   d 
                    
                   
                       
                   
                    
                   Ω 
                 
               
             
           
         
       
     
     wherein: K standard  is a stiffness matrix due to the plurality of finite elements; j is a finite element; E j  is a region covered by each of the finite elements; ∇N i  is a spatial gradient of a function N i ; ∇N j  is a spatial gradient of a function N j ; and dΩ is a boundary of a region covered by each finite element. 
     The stiffness contribution due to the overlapping finite elements may calculated according to the expression: 
     
       
         
           
             
               K 
               overlapping 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∑ 
                   
                     k 
                     ≠ 
                     j 
                   
                 
                  
                 
                   
                     ∫ 
                     
                       
                         E 
                         j 
                       
                       ⋂ 
                       
                         E 
                         k 
                       
                     
                   
                    
                   
                     
                       Θ 
                       j 
                     
                      
                     
                       Θ 
                       k 
                     
                      
                     
                       
                         ∇ 
                         
                           N 
                           j 
                         
                       
                       · 
                       
                         ∇ 
                         
                           N 
                           k 
                         
                       
                     
                      
                     d 
                      
                     
                         
                     
                      
                     Ω 
                   
                 
               
             
           
         
       
     
     wherein: K overlapping  is a stiffness matrix of the overlapping elements; j is a finite element; k is a finite element overlapping finite element j; E j  ∩ E k  is an overlapping region between finite elements, j and k; Θ j  is the orientation of the first finite element j; Θ k  is the orientation of the second finite element k; ∇N j  is a spatial gradient of a function N j ; ∇N k  is a spatial gradient of a function N k ; and dΩ is the boundary of a region covered by the elements, j and k. 
     The step of combining the stiffness of the finite elements and the stiffness of the overlapping finite elements may include the additional step of summing the stiffness of the plurality of finite elements and the stiffness of the overlapping finite elements. 
     In accordance with a further aspect of the present invention, a method is provided for carrying out finite element analysis. The method includes the steps of meshing a domain under a field with a plurality of finite elements and detecting if each finite element overlaps a finite element. A stiffness contribution due to the plurality of finite elements is calculated and a stiffness contribution due to the overlapping finite elements are calculated. The stiffness contribution due to the plurality of finite elements and the stiffness contribution due to the tangled finite elements are combined. 
     Each of the finite elements is defined by a plurality of nodes and each finite element has an orientation. Each node has a nodal shape function and the nodal shape functions are defined according to the expression: 
     
       
         
           
             
               
                 φ 
                 i 
               
                
               
                 ( 
                 • 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   j 
                   ∈ 
                   
                     C 
                      
                     
                       ( 
                       i 
                       ) 
                     
                   
                 
               
                
               
                 
                   Θ 
                   j 
                 
                  
                 
                   
                     N 
                     
                       i 
                       , 
                       j 
                     
                   
                    
                   
                     ( 
                     • 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein ϕ i (⋅) are the nodal shape functions; j is a finite element; i is a node; C♭i? is a set of finite elements connected to node i; Θ j  is the orientation of finite element j; and N i,j (⋅) are element shape functions. 
     The stiffness contribution due to the plurality of finite elements may be calculated according to the expression: 
     
       
         
           
             
               K 
               standard 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∫ 
                   
                     E 
                     j 
                   
                 
                  
                 
                   
                     
                       ( 
                       
                         ▽N 
                         i 
                       
                       ) 
                     
                     · 
                     
                       ( 
                       
                         ▽N 
                         j 
                       
                       ) 
                     
                   
                    
                   d 
                    
                   
                       
                   
                    
                   Ω 
                 
               
             
           
         
       
     
     wherein: K standard  is a stiffness matrix of the plurality of finite elements; j is a finite element; E j  is a region covered by each of the finite elements; ∇N i  is a spatial gradient of a function N i ; ∇N j  is a spatial gradient of a function N j ; and dΩ is a boundary of a region covered by each finite element. 
     The stiffness contribution of the overlapping finite elements may be calculated according to the expression: 
     
       
         
           
             
               K 
               overlapping 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∑ 
                   
                     k 
                     ≠ 
                     j 
                   
                 
                  
                 
                   
                     ∫ 
                     
                       
                         E 
                         j 
                       
                       ⋂ 
                       
                         E 
                         k 
                       
                     
                   
                    
                   
                     
                       Θ 
                       j 
                     
                      
                     
                       Θ 
                       k 
                     
                      
                     
                       
                         
                           ▽ 
                            
                           N 
                         
                         j 
                       
                       · 
                       
                         
                           ▽ 
                            
                           N 
                         
                         k 
                       
                     
                      
                     d 
                      
                     
                         
                     
                      
                     Ω 
                   
                 
               
             
           
         
       
     
     wherein: K overlapping  is a stiffness matrix of the overlapping elements; j is a finite element; k is a finite element overlapping finite element j; E j  ∩ E k  is an overlapping region between first and second finite elements, j and k; Θ j  is the orientation of the first finite element j; Θ k  is the orientation of the second finite element k; ∇N j  is a spatial gradient of a function N j ; ∇N k  is a spatial gradient of a function N k ; and dΩ is the boundary of a region covered by the elements, j and k. 
     The step of combining the stiffness of the plurality of finite elements and the stiffness of the tangled finite elements includes the step of summing the stiffness of the plurality of finite elements and the stiffness of the tangled finite elements. In addition, the step of detecting if each finite element is tangled includes the step of determining if a finite element overlaps an adjacent finite element and if the finite element overlaps the adjacent finite element, determining that the finite element is tangled. 
     In accordance with a still further aspect of the present invention, a method is provided for carrying out finite element analysis. The method includes the steps of meshing a domain under a field with a plurality of finite elements and determining a stiffness contribution due to the plurality of finite elements. A stiffness contribution due to a subset of the plurality of finite elements is determined and combined with the stiffness contribution due to the plurality of finite elements. 
     It is contemplated for the finite element to be triangle. Each triangle is defined by a first point, a second point and a third point in order. The first and second points define a first segment, the second and third points define a second segment and the third and first points define a third segment. The first, second and third segments define an interior of the triangle. Alternatively, each finite element may be a tetrahedron. Each tetrahedron is defined by a first point, a second point, a third point and a fourth point in order. The first, second and third points define a first plane; the second, third, and fourth points define a second plane; the third, fourth, and first points define a third plane; and the fourth, first and second points define a fourth plane. The first, second, third and fourth planes define an interior of the tetrahedron. The orientation of each finite element is determined by computing a determinant of the Jacobian. 
     Each of the finite elements is defined by a plurality of nodes and each finite element has an orientation. In addition, each node has a nodal shape function. The nodal shape functions are defined according to the expression: 
     
       
         
           
             
               
                 φ 
                 i 
               
                
               
                 ( 
                 • 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   j 
                   ∈ 
                   
                     C 
                      
                     
                       ( 
                       i 
                       ) 
                     
                   
                 
               
                
               
                 
                   Θ 
                   j 
                 
                  
                 
                   
                     N 
                     
                       i 
                       , 
                       j 
                     
                   
                    
                   
                     ( 
                     • 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein ϕ i (⋅) are the nodal shape functions; j is a finite element; i is a node; C♭i? is a set of finite elements connected to node i; Θ j  is the orientation of the finite element j; and N i,j (⋅) are element shape functions. 
     The stiffness contribution due to the plurality of finite elements may be calculated according to the expression: 
     
       
         
           
             
               K 
               standard 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∫ 
                   
                     E 
                     j 
                   
                 
                  
                 
                   
                     
                       ( 
                       
                         ▽N 
                         i 
                       
                       ) 
                     
                     · 
                     
                       ( 
                       
                         ▽N 
                         j 
                       
                       ) 
                     
                   
                    
                   d 
                    
                   
                       
                   
                    
                   Ω 
                 
               
             
           
         
       
     
     wherein: K standard  is a stiffness matrix of the plurality of finite elements; j is a finite element; E j  is a region covered by each of the finite elements; ∇N j  is a spatial gradient of a function N i ; ∇N j  is a spatial gradient of a function N j ; and dΩ is a boundary of a region covered by each finite element. 
     The stiffness contribution due to the subset of finite elements may calculated according to the expression: 
     
       
         
           
             
               K 
               overlapping 
             
             = 
             
               
                 ∑ 
                 j 
               
                
               
                 
                   ∑ 
                   
                     k 
                     ≠ 
                     j 
                   
                 
                  
                 
                   
                     ∫ 
                     
                       
                         E 
                         j 
                       
                       ⋂ 
                       
                         E 
                         k 
                       
                     
                   
                    
                   
                     
                       Θ 
                       j 
                     
                      
                     
                       Θ 
                       k 
                     
                      
                     
                       
                         
                           ▽ 
                            
                           N 
                         
                         j 
                       
                       · 
                       
                         
                           ▽ 
                            
                           N 
                         
                         k 
                       
                     
                      
                     d 
                      
                     
                         
                     
                      
                     Ω 
                   
                 
               
             
           
         
       
     
     wherein: K overlapping  is a stiffness matrix due to the subset of finite elements; j is a finite element; k is a finite element overlapping finite element j; E j  ∩ E k  is an overlapping region between finite elements, j and k; Θ j  is the orientation of the finite element j; Θ k  is the orientation of the finite element k; ∇N j  is a spatial gradient of a function N j ; ∇N k  is a spatial gradient of a function N k ; and dΩ is the boundary of a region covered by the elements, j and k. 
     The step of combining the stiffness of the finite elements and the stiffness of the subset of finite elements may include the additional step of summing the stiffness of the plurality of finite elements and the stiffness of the subset of finite elements. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The drawings furnished herewith illustrate a preferred construction of the present invention in which the above advantages and features are clearly disclosed as well as others which will be readily understood from the following description of the illustrated embodiment. 
       In the drawings: 
         FIG. 1  is a schematic view of an untangled finite element mesh; 
         FIG. 2  is a schematic view of a tetrahedral finite element; 
         FIG. 3  is a schematic view of a tangled finite element mesh; and 
         FIG. 4  is a schematic view of an invalid, tangled finite element mesh. 
     
    
    
     DETAILED DESCRIPTION OF THE DRAWINGS 
     In classic FEA, one attempts to verify that a proposed design of a product will perform in accordance with desired specifications by approximating a solution to an underlying field equation, such as the Poisson equation: 
     
       
         
           
             
               
                 
                   
                     ∇ 
                     
                       · 
                       
                         ∇ 
                         u 
                       
                       · 
                     
                   
                   = 
                   
                     
                       - 
                       f 
                     
                      
                     
                         
                     
                      
                     in 
                      
                     
                         
                     
                      
                     Ω 
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     1 
                     ) 
                   
                 
               
             
             
               
                 
                   
                     u 
                      
                     
                        
                       
                         ∂ 
                         Ω 
                       
                     
                   
                   = 
                   0 
                 
               
               
                 
                     
                 
               
             
           
         
       
     
     wherein: u is the unknown field, such as temperature; ƒ is a body heat; Ω is a region; and ∂Ω is a boundary of the region of the product. Other field problems such as structural or magnetic may be posed in a similar form. The first step towards solving Equation (1) is to create a finite element mesh over the domain Ω. For simplicity, only simplicial meshes, i.e., 2-D triangular,  FIGS. 1 , and 3-D tetrahedral elements are considered hereafter. However, other types of meshes may be utilized without deviating from the scope of the present invention. Heretofore, in classic FEA, the mesh needed to satisfy certain topological and geometrical properties. For example, in 2-D triangular meshes, the following conditions must be satisfied: 
     1. Every element must be defined by exactly three nodes. 
     2. Every element must be positively oriented, i.e., the determinant of the Jacobian must be positive. 
     3. No two elements must overlap. 
     4. Every node must be connected to at least one element. 
     Meshes which satisfy the above-noted condition are referred to as being valid. 
     In FEA, approximate functions are constructed by seeking solutions of the form: 
     
       
         
           
             
               
                 
                   
                     u 
                      
                     
                       ( 
                       • 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       i 
                     
                      
                     
                       
                         
                           u 
                           ^ 
                         
                         i 
                       
                        
                       
                         
                           φ 
                           i 
                         
                          
                         
                           ( 
                           • 
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     2 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: û i  is the degree of freedom associated with node i, and ϕ i (⋅) are nodal shape functions or “hat functions” that take a value of 1 at that node, and go to zero beyond the set of elements connected to node i,  FIG. 2 . It is noted that in Equation (2), ϕ i (⋅) is an abbreviated notation for ϕ i (x, y) in 2D, and ϕ i (x, y, z) in 3D. 
     In addition, the hat function is usually decomposed into disjoint element shape functions N i,j (⋅) defined over each element j connected to node i. Unfortunately, as hereinafter described, when a mesh is tangled or in other words when the above conditions are not satisfied, the above process leads to erroneous results. 
     For a 2-D simplicial element,  FIG. 1 , such as a triangle, the phrase ‘orientation of an element’ may be defined by three points, P 0 , P 1 , and P 2  (in that order) such that if the interior of the element lies to the left of the oriented segment P 0 P 1 , (or equivalently P 1 P 2  or P 2 P 0 ), the triangle is positively oriented. Otherwise, the triangle is negatively oriented. For a 3-D simplicial element,  FIG. 2 , such as a tetrahedron, the phrase ‘orientation of an element’ may be defined by four points, P 0 , P 1 , P 2  or P 3  (in that order) such that if the interior of the element lies to the left of the oriented plane P 0 P 1 P 2  (or equivalently P 1 P 2 P 3 , P 2 P 3 P 0 or P 3 P 0 P 1 ), the tetrahedron is positively oriented. Otherwise, the tetrahedron is negatively oriented. It can be appreciated that the above definitions of ‘orientation of an element’ for simplicial elements is equivalent to other definitions of the ‘orientation of an element’ involving the determinant of the Jacobian. It can be understood that in a valid mesh, as heretofore described, all elements are positively oriented. 
     On the other hand, a mesh is tangled if: (a) one or more elements have a negative orientation; and (b) the mesh can be untangled, i.e., the mesh may be converted into a valid mesh with all positively oriented elements through some hypothetical movement of nodes. It is noted that the definition of a tangled mesh only requires the existence of node movement. This, in turn, eliminates meshes that are topologically corrupt. Further, it is assumed that the tangling of the mesh is internal, i.e., that the overall geometry defined by the tangled mesh is exactly the same as the original geometry Ω defined via the valid mesh. Thus, the mesh shown in  FIG. 4  wherein a node goes beyond the boundary of the original geometry is not considered here. 
     Referring to  FIG. 3 , in a tangled mesh, it is noted there must exist a pair of neighboring elements  10  and  12  with opposite orientation, and such pairs must overlap, i.e., there are points  14  within Ω that belong to both elements. The boundary  16  shared by two oppositely oriented elements is of significant importance. More specifically, the common boundary of neighboring and oppositely oriented elements is referred to as a fold. Since only simplicial elements are being considered, the fold must necessarily occur at the common boundary. 
     If a classic FEA approximation, as in Equation (2), is sought over a tangled mesh, the nodal shape function ϕ i (⋅) can no longer be decomposed into disjoint element shape functions due to the guaranteed presence of overlapping elements. As a result, the required continuity of the functional space is lost. While it is possible to ignore this violation and carry out a finite element assembly, followed by a solution to obtain the nodal values over the tangled mesh, erroneous results appear later in numerical results. 
     In order to handle tangled meshes, the methodology of the present invention contemplates a modification to classic FEA. More specifically, for any node i in a tangled mesh, C(i) may be defined as the set of elements connected to it. As such, the nodal shape functions maybe defined as: 
     
       
         
           
             
               
                 
                   
                     
                       φ 
                       i 
                     
                      
                     
                       ( 
                       • 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         j 
                         ∈ 
                         
                           C 
                            
                           
                             ( 
                             i 
                             ) 
                           
                         
                       
                     
                      
                     
                       
                         Θ 
                         j 
                       
                        
                       
                         
                           N 
                           
                             i 
                             , 
                             j 
                           
                         
                          
                         
                           ( 
                           • 
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     3 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein ϕ i (⋅) are the nodal shape functions; j is a finite element; i is a node; C♭i? is a set of finite elements connected to node i; Θ j  is the orientation of the finite element j; and N i,j (⋅) are element shape functions. 
     It is noted that if all elements connected to node i are positively or negatively oriented, then Equation (3) is the classic hat function, heretofore described. In the event that a portion of the elements in C♭i? are negatively oriented, while others are positively oriented, it can be proven that: 
     1. ϕ i (⋅)=0 over ∂Ω i ; and 
     2. ϕ i (⋅) is continuous in Ω i . 
     First, consider the statement: ϕ i (⋅)=0 over ∂Ω i . Since only simplicial elements are being considered, as heretofore described, points on ∂Ω i must also lie on the boundary of one or more elements in C♭i?. If a point lies on the boundary of a single element (example: point p in  FIG. 3 ), then, by definition, the element shape function of that element goes to zero at that point. Further, per Equation (3), the nodal shape function ϕ i (⋅) depends only on that element at that point, and therefore, also goes to zero at that point. Now, suppose a point on ∂Ω i  lies on the boundary of more than one element in C(i) (example: point q in  FIG. 3 ), then the point must lie on a fold. Since the element shape functions are defined to be continuous across element boundaries, and since the two elements have opposite orientations, then, per Equation (3), their net contribution is zero, i.e., ϕ i (⋅) goes to zero at point q. 
     Now consider the statement: ϕ i (⋅) is continuous in Ω i . Here again, there are only two cases to consider. First, referring to  FIG. 3 , consider the continuity at a point r that lies in the interior of one or more elements. Since the element shape functions N i,j (⋅) are defined to be continuous in the interior of elements and since ϕ i (⋅) is linear combination of such functions, the continuity of the latter is guaranteed. Second, consider the continuity at a point s that lies in the interior of Ω i , but on the boundary of the elements. The element shape functions N i,j (⋅) at the boundary of neighboring elements are guaranteed to continuous, once again by construction. Therefore, since ϕ i (⋅) is a linear sum of such functions, it is also continuous. Finally, by generalization, any point that falls into both of the above categories, i.e., lies in the interior of certain elements and lies on the boundary of other elements, is also continuous. 
     In order to implement the methodology of the present invention, it is contemplated not to use the nodal shape functions ϕ i (⋅) directly. Instead, element shape functions N i,j (⋅) may be used for reasons of computational ease and efficiency. More specifically, the field (e.g. stress, thermal or magnetic fields) over a mesh may be expressed at any point as the sum over all elements that contain that point: 
     
       
         
           
             
               
                 
                   
                     u 
                      
                     
                       ( 
                       • 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       j 
                     
                      
                     
                       
                         u 
                         j 
                       
                        
                       
                         ( 
                         • 
                         ) 
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     4 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: u(⋅) is the field (e.g. stress, thermal or magnetic fields) over a mesh; and u j (⋅) the field over an element j. 
     The field over an element j (u j (⋅)) is simply the oriented sum of contributions from all nodes i that define that element and may be calculated according to the expression: 
     
       
         
           
             
               
                 
                   
                     
                       u 
                       j 
                     
                      
                     
                       ( 
                       • 
                       ) 
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                          
                         
                           p 
                           ∈ 
                           
                             E 
                             j 
                           
                         
                       
                     
                      
                     
                       
                         Θ 
                         j 
                       
                        
                       
                         
                           ∑ 
                           
                             i 
                              
                             
                               j 
                               ∈ 
                               
                                 C 
                                  
                                 
                                   [ 
                                   i 
                                   ] 
                                 
                               
                             
                           
                         
                          
                         
                           
                             
                               u 
                               ^ 
                             
                             i 
                           
                            
                           
                             
                               N 
                               
                                 i 
                                 , 
                                 j 
                               
                             
                              
                             
                               ( 
                               • 
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     5 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: u j (⋅) the field over an element j; Θ j  is the orientation of the finite element j; E j  is the region of element j; i is a node; C♭i? is a set of finite elements connected to node i; û i  is the degree of freedom associated with node i; and N i,j (⋅) are element shape functions. 
     The stiffness contribution due to each element, whether positively or negatively orientated, is calculated and assembled according to the standard expression: 
     
       
         
           
             
               
                 
                   
                     K 
                     standard 
                   
                   = 
                   
                     
                       ∑ 
                       j 
                     
                      
                     
                       
                         ∫ 
                         
                           E 
                           j 
                         
                       
                        
                       
                         
                           
                             ( 
                             
                               ▽N 
                               j 
                             
                             ) 
                           
                           · 
                           
                             ( 
                             
                               ▽N 
                               j 
                             
                             ) 
                           
                         
                          
                         d 
                          
                         
                             
                         
                          
                         Ω 
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     6 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: K standard  is the stiffness matrix of the entire set of elements; j is a finite element; E j  is a region covered by each element; ∇N j  is the spatial gradient of the function N i ; ∇N j  is the spatial gradient of the function N j ; and dΩ is an infinitesimal region. 
     The stiffness contribution due to each pair of overlapping elements is calculated and assembled according to the expression: 
     
       
         
           
             
               
                 
                   
                     K 
                     overlapping 
                   
                   = 
                   
                     
                       ∑ 
                       j 
                     
                      
                     
                       
                         ∑ 
                         
                           k 
                           ≠ 
                           j 
                         
                       
                        
                       
                         
                           ∫ 
                           
                             
                               E 
                               j 
                             
                             ⋂ 
                             
                               E 
                               k 
                             
                           
                         
                          
                         
                           
                             Θ 
                             j 
                           
                            
                           
                             Θ 
                             k 
                           
                            
                           
                             
                               ∇ 
                               
                                 N 
                                 j 
                               
                             
                             · 
                             
                               ∇ 
                               
                                 N 
                                 k 
                               
                             
                           
                            
                           d 
                            
                           
                               
                           
                            
                           Ω 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     7 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: K overlapping  is the stiffness matrix due to each pair of overlapping elements; j is a finite element; k is a finite element overlapping finite element j; E j  ∩ E k  is an overlapping region between the elements, j and k; Θ j  is the orientation of the finite element j; Θ k  is the orientation of the finite element k; ∇N j  is the spatial gradient of the function N j ; ∇N k  is the spatial gradient of the function N k ; and dΩ is an infinitesimal region. It may be appreciated that in an untangled mesh, no two elements overlap, and therefore, the K overlapping  matrix is zero. 
     The combined stiffness matrix is calculated according to the expression: 
         K=K   standard   +K   overlapping    Equation (8)
 
     wherein: K is the combined stiffness; K standard  is the stiffness of each element; and K overlapping  is the stiffness due to each pair of overlapping elements. 
     The external force over a tangled mesh is computed according to the expression: 
     
       
         
           
             
               
                 
                   
                     f 
                     oriented 
                   
                   = 
                   
                     
                       ∑ 
                       j 
                     
                      
                     
                       
                         ∫ 
                         
                           E 
                           j 
                         
                       
                        
                       
                         
                           Θ 
                           j 
                         
                          
                         
                           N 
                           j 
                         
                          
                         fd 
                          
                         
                             
                         
                          
                         Ω 
                       
                     
                   
                 
               
               
                 
                   Equation 
                    
                   
                       
                   
                    
                   
                     ( 
                     9 
                     ) 
                   
                 
               
             
           
         
       
     
     wherein: ƒ oriented  is the external force over a tangled mesh; j is a finite element; E j  is a region covered by each element; Θ j  is the orientation of the finite element j; ƒ is the body force; and dΩ is an infinitesimal region. It can be appreciated that if the orientation is all positive, then Equation (9) results in the standard body force. 
     It can be appreciated that the terms of above equations may be rewritten in the matrix form as follows: 
         Kû=ƒ   oriented    Equation (10)
 
     wherein: K is the combined stiffness, as heretofore described; û is the collection of all nodal degrees of freedom; and ƒ oriented  is external force over a tangled mesh. 
     The methodology, as heretofore described, is an extension to the underlying mathematics of the classic finite element formulation. The methodology of the present invention allows for FEA to be used in conjunction with tangled meshes that were previously considered unacceptable. In addition, the methodology can be easily incorporated into classic FEA with minor modifications. 
     Various modes of carrying out the invention are contemplated as being within the scope of the following claims particularly pointing out and distinctly claiming the subject matter, which is regarded as the invention.