Adaptive Creation of OrthogonalAnisotropic Triangular Meshes Using TargetMatricesbyJose David Zuniga VazquezB.Sc., Instituto Tecnologico de Hermosillo, 2009A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2015c© Jose David Zuniga Vazquez 2015AbstractWe present a new procedure to adaptively produce anisotropic metric-orthogonal meshes in this thesis. Theapproach is based on mesh optimization techniques: point insertion designed to improve mesh alignment aswell as conforming to a metric, swapping in the metric space, point movement defined by target elementsfrom a metric, and point deletion based on quality and metric length. These techniques are intended toproduce quasi-structured meshes which have the advantages of flexibility for complex geometries of un-structured meshes and the directional accuracy of structured meshes. The result is reliable alignment foranisotropic meshes reducing our previous work’s reliance on smoothing for good alignment. The methodol-ogy is implemented in 2D and extended to 3D. Examples of analytical metrics and error estimation metricson a numerical simulation of a flow are shown.iiPrefaceAll of the work presented henceforth was an intellectual product of a working relationship between JoseDavid Zuniga Vazquez and Dr. Carl Ollivier-Gooch. All the data analysis, concept formation, implemen-tation of the methods, and the manuscript preparations were done by Jose David Zuniga Vazquez withinvaluable guidance from Carl Ollivier-Gooch throughout the process.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Objective and Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Mesh Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Mesh Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Metric-Based Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4 Target Matrix Optimization Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.1 2D Mesh Optimization Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 2D General Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3 Tri-to-Quad Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4 3D Mesh Optimization Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.5 3D General Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.1 2D Analytic Metric Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.2 2D Flow Around an Airfoil Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.3 Tri-to-Quad Conversion Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66ivTable of Contents4.4 3D Analytic Metric Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 835.1 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 835.2 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86vList of Tables3.1 Mesh size during 2D adaptation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 Mesh size during 3D adaptation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.1 Mesh size in adaptation for a quarter circle curved anisotropic feature in 2D . . . . . . . . . 554.2 Mesh size in adaptation for perpendicular anisotropic features in 2D . . . . . . . . . . . . . 574.3 Comparison of lift and drag for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦ . . . . . . . . . 644.4 Mesh size in adaptation for a quarter circle curved anisotropic feature in 3D . . . . . . . . . 704.5 Mesh size in adaptation for perpendicular anisotropic features in 3D . . . . . . . . . . . . . 76viList of Figures1.1 Smoothed and non smoothed for a subsonic viscous flow over a NACA 0012, Ma= 0.5, Re=5000,α = 0◦. Figures courtesy of M. Sharbatdar [67] . . . . . . . . . . . . . . . . . . . . . 41.2 Close up view of top and bottom surface. Figures courtesy of M. Sharbatdar [67] . . . . . . 41.3 Close up view of trailing edge. Figures courtesy of M. Sharbatdar [67] . . . . . . . . . . . 52.1 Classification of common meshes in 2D by elements . . . . . . . . . . . . . . . . . . . . . 62.2 Classification of common meshes in 3D by elements . . . . . . . . . . . . . . . . . . . . . 72.3 Mixed element meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.4 Structured mesh and unstructured mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.5 Comparison between unstructured, quasi structured, and structured meshes . . . . . . . . . 92.6 Comparison between isotropic and anisotropic meshes in 2D . . . . . . . . . . . . . . . . . 102.7 Classification of skinny tetrahedra collapsing to an edge . . . . . . . . . . . . . . . . . . . . 112.8 Classification of skinny tetrahedra collapsing to a plane . . . . . . . . . . . . . . . . . . . . 112.9 Configurations suitable for swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.10 Configurations not suitable for swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.11 Configurations for a 2-3 / 3-2 swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.12 Configurations for a 5-6 / 6-5 swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.13 Configurations for a 4-4 swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.14 Configurations not suitable for swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.15 Insertion of a vertex in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.16 Insertion of a vertex in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.17 Edge insertion: 4 to 8 tetrahedra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.18 Deleting a vertex in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.19 Deleting a vertex from interior cells in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.20 Deleting a vertex from a boundary cells in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 172.21 Smoothing a vertex in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.22 Smoothing a vertex in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.23 Deformation of tensors Fp and Fq . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.24 Distance and angle measured for a certain point . . . . . . . . . . . . . . . . . . . . . . . . 202.25 Weighted Jacobian matrix definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.26 Target matrix in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.27 Target matrix in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24viiList of Figures2.28 Possible 2D configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.29 Possible 3D configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.30 Quality with different metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.1 2D swapping in metric space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Insertion for a long boundary edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3 Insertion for a long interior edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.4 Deletion of a vertex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.5 Representation of quarter circle metric function behavior . . . . . . . . . . . . . . . . . . . 353.6 2D adaptation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.7 1st example of combining neighboring triangular cells into a quad cell . . . . . . . . . . . . 383.8 2nd example of combining neighboring triangular cells into a quad cell . . . . . . . . . . . . 393.9 3D swapping configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.10 3D maxmin sine in physical space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.11 3D maxmin sine in metric space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.12 3D volume-length in physical space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.13 3D volume-length in metric space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.14 3D insertion on an interior edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.15 3D insertion on a boundary edge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.16 Process of searching for conflicts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.17 3D adaptation process: initial mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.18 3D adaptation process: 1st quality insertion (2nd step) . . . . . . . . . . . . . . . . . . . . 493.19 3D adaptation process: 2nd quality insertion (5th step) . . . . . . . . . . . . . . . . . . . . 503.20 3D adaptation process: long edge insertion (6th step) . . . . . . . . . . . . . . . . . . . . . 513.21 3D adaptation process: swapping after coarsening (8th step) . . . . . . . . . . . . . . . . . 524.1 Initial mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.2 1st adaptation for a quarter circle curved anisotropic feature in 2D with a close up . . . . . . 544.3 2nd adaptation for a quarter circle curved anisotropic feature in 2D with a close up . . . . . 544.4 1st adaptation, Shape Orientation quality and histogram for a quarter circle curved anisotropicfeature in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.5 2nd adaptation, Shape Orientation quality and histogram for a quarter circle curved anisotropicfeature in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.6 1st adaptation for perpendicular anisotropic features in 2D with a close up . . . . . . . . . . 564.7 2nd adaptation for perpendicular anisotropic features in 2D with a close up . . . . . . . . . 564.8 1st adaptation, Shape Orientation quality and histogram for perpendicular anisotropic fea-tures in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.9 2nd adaptation, Shape Orientation quality and histogram for perpendicular anisotropic fea-tures in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.10 Initial mesh around a NACA 0012 airfoil with 2543 vertices and 4941 cells . . . . . . . . . 58viiiList of Figures4.11 Mesh and solutions for a subsonic viscous flow over a NACA 0012, Ma= 0.5, Re= 5000,α =0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.12 Closeup of airfoil NACA 0012 after 3rd adaptation . . . . . . . . . . . . . . . . . . . . . . 614.14 Comparison between smoothed and non-smoothed mesh after 3rd adaptation . . . . . . . . 624.13 ShapeOrient quality after 3rd adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.15 Drag convergence for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦ . . . . . . . . . . . . . . . 644.16 Distribution behavior per adaptation for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦ . . . . . 654.17 Distribution behavior for top and bottom surfaces of 3rd Adaptation for NACA 0012, Ma =0.5, Re = 5000,α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.18 Tri-to-quad for curved anisotropic feature . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.19 Tri-to-quad for perpendicular anisotropic features . . . . . . . . . . . . . . . . . . . . . . . 674.20 Tri-to-quad for subsonic viscous flow over a NACA 0012, Ma = 0.5, Re = 5000,α = 0◦ . . 684.21 Initial mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.22 Close up of 1st adaptation quarter circle curved anisotropic feature in 3D . . . . . . . . . . . 714.23 Close up of 2nd adaptation quarter circle curved anisotropic feature in 3D . . . . . . . . . . 714.24 1st adaptation for a quarter circle curved anisotropic feature in 3D . . . . . . . . . . . . . . 724.25 2nd adaptation for a quarter circle curved anisotropic feature in 3D . . . . . . . . . . . . . . 734.26 ShapeOrient quality for 1st adaptation of quarter circle curved anisotropic feature in 3D . . . 744.27 ShapeOrient quality for 2nd adaptation of quarter circle curved anisotropic feature in 3D . . 754.28 2nd adaptation for perpendicular anisotropic features in 3D . . . . . . . . . . . . . . . . . . 774.29 Close up of 2nd adaptation perpendicular anisotropic features in 3D . . . . . . . . . . . . . 784.30 4th Adaptation for perpendicular anisotropic features in 3D . . . . . . . . . . . . . . . . . . 794.31 Close up of 4th adaptation perpendicular anisotropic features in 3D . . . . . . . . . . . . . 804.32 ShapeOrient quality for 2nd adaptation of perpendicular anisotropic features in 3D . . . . . 814.33 ShapeOrient quality for 4th adaptation of perpendicular anisotropic features in 3D . . . . . . 82ixNomenclature∆ Local aspect ratio parameterCDν Viscosity drag coefficientCDp Pressure drag coefficientC f Skin friction factorCL Lift coefficientCP Pressure coefficientl LengthlM Metric lengthO Local orientation parameterR Unitary eigenvector matrixS Local shape parameterU Local shape and size parameterΓ Radious metric componentΛ Diagonal matrix of eigenvaluesλ EigenvalueΦ Local angle parameterΘ Angle metric componentϒ Quality metric criteria~V EigenvectorΞ Local size parameterA Jacobian matrixC CellxNomenclatureE Edge of a cellF Frobenius matrix normM Metric tensorMa Mach numberQ Cell qualityT Weighted Jacobian matrixV Vertex of a cellW Target matrixx,y,z Cartesian coordinatesxiAcknowledgementI would like to express my gratitude to my supervisor Carl Ollivier-Gooch for all his support throughout myresearch. I’m thankful for his aspiring guidance, invaluably constructive criticism and experienced advices.I also want to thank Mahkame Sharbatdar and the rest of my lab-mates for their advices and support duringthe different struggles that I had related to research progress.xiiDedicationI want to thank my family for their support, my parents Norma and Juan for their constant support andmotivation from a far away country, my brother Miguel for his guidance and advice at the beginning of thisjourney, and my brother Daniel for reminding me not to give up.I would like to make a special dedication for my wife / love of my life Gaby, who left her home, herfriends, and family to come with me. She’s the main reason why I have accomplished my goals.xiiiChapter 1IntroductionThe numerical solution of the conservation equations for fluids is known as computational fluid dynamics(CFD). The goal of the CFD field is to analyze the physical characteristics and behavior of fluid flow in adomain. This is related to the action and interaction of phenomena such as dissipation, diffusion, convection,turbulence, shock waves, and boundary layers. These phenomena are governed by the compressible Navier-Stokes equations; further detail is given in [39]. CFD is commonly used in today’s industry in a wide varietyof disciplines such as medical research, aerospace, fluid dynamics, meteorology, and more. The advantagesof CFD have already been proven including reduced reliance on physical experimentation, improved designreliability, and reduced of cost of product and process development, among others.To solve problems in CFD, we must define flow physics and problem geometry. The physics definitionmust include models, initial conditions, and boundary conditions. Once the problem geometry is defined,often using computer aided design (CAD) software, this domain must be discretized or meshed. Thereare two fundamentally different approaches, structured meshes characterized by regular connectivity andunstructured meshes characterized by irregular connectivity, which are described in Sec.2.1. Once a meshis defined, the equations are discretized over the mesh and the resulting system of equations is solved. Aftera satisfactory iterative convergence is reached, the solution is analyzed in a postprocessor (more details aredescribed by D. L. Davidson [13]).There are many areas of on going research within CFD; the area of interest for this research is un-structured mesh adaptation. Many CFD simulations have common where small regions of the domain re-quire higher mesh density. Common features include shock waves, boundary layers, and ignition frontsamong others. To improve the efficiency and accuracy of the numerical solution, a higher density ofpoints would be expected in these regions. Many different mesh adaptation schemes have been developed[1, 12, 57, 60, 68, 69, 72, 73, 75] for this purpose.Essential to mesh adaptation is the ability to control the size, shape and orientation of mesh elementsthroughout the domain. If the scheme adapts the size of the mesh elements and maintains the original shape(in case of triangular elements an equilateral shape) the adaptation scheme is an isotropic mesh adaptationscheme. Isotropic mesh adaptation schemes use more mesh elements than needed in regions with largesolution errors where the solution changes more significantly in one direction than the other. To handlethis problem, anisotropic mesh adaptation schemes can adjust the size, shape and orientation of the meshelements according to the behavior of the physical solution (which usually means large aspect ratio meshelements). More details are described by W. Huang [26].11.1. Motivation1.1 MotivationIn the context of numerical solutions, researchers are continuously investigating numerical strategies toresolve problems while maintaining good accuracy at a low computational cost, with a goal of reducingthe burden of mesh generation. For a wide range of engineering problems, a mesh that is adaptive in size,shape and orientation can result in accurate solutions for problems with complex geometries and rapidlychanging solutions. For this reason many researchers opt for using mesh adaptation schemes. Unstructuredmesh adaptation has largely proved its efficiency for improving the accuracy of the numerical solution andto capture the behavior of physical phenomena in complex domains better than structured meshes [46, 75].When the governing equations are discretized across the mesh, numerical fluxes are usually computedfor each face of the control volume by an approximate Riemann solver. These fluxes are usually calculated ina reference frame that is aligned with the mesh. For well suited structured meshes it is possible in principleto align faces to be either tangent or normal to a flow feature, but such meshes are difficult to generate.Unstructured meshes are easier to generate, but physical features may misrepresent since the control volumefaces may be at arbitrary angles to the flow features; this lack of alignment degrades accuracy. Minimumerror should occur when the normal of the control volume interface coincides with the normal of the flowfeature with a small spacing. A lack of alignment of the mesh elements may cause accuracy problems thatstill need to be handled; more details are shown in [24].An approach for handling this problem, which this thesis is focused on, is by applying mesh modificationtechniques to align the mesh with the physical features. We adaptively modify a mesh based on a metric,obtained analytically or defined by the corresponding solution, so the mesh becomes aligned with solutionfeatures. Most CFD related problems do not have an analytical solution, so it is necessary to determine anappropriate error estimate. This approach to mesh adaptation requires computation of a suitable error, whichwill drive the adaptation scheme to modify the mesh.One of the most difficult aspects of mesh adaptation is finding a good criterion for adapting the mesh;approaches for this include featured-based indicators, adjoint methods, and truncation errors, among others[66]. Feature-based adaptation methods often use solution gradients, solution curvature, or solution features[5, 15, 66]. This method often results in some feature being over resolved while others are not resolvedenough, producing an increase in the solution error. Adjoint-based adaptation can provide targeted meshadaptation due to an estimate of local contribution of each cell to the error in some output quantity, but itscomplexity and code intrusivenesses is a major drawback [10, 44, 50, 59]. Truncation error-based adaptationaims to equally distribute the truncation error (difference between the partial differential equation and itsdiscrete approximation) over the entire domain to reduce the total discretization error [11]. Metric-basedadaptation transforms an error estimate into a metric tensor field, for which the mesh elements can be adaptedto match the metric in a later process [7, 41, 56, 67].Efficient adaptation of a high quality anisotropic mesh with low computational cost is a non-trivialproblem. The advantages of metric based adaptation schemes are well known [1, 8, 25, 38, 40, 42, 56, 68].By transforming an error estimate into a metric tensor field, the mesh elements can be adapted to matchthe metric and thus improve accuracy. For this purpose different adaptation schemes have been developed.An approach which uses an advancing-front type point placement has been used for mesh adaptation in21.2. Objective and Outlinedifferent schemes [9, 62]. S.H. Lo [38] proposed a simple anisotropic refinement procedure for three-dimensional tetrahedral mesh based on successive bisection of edges. Refinement is done by dividing thelines of the mesh based on the lengths calculated by a metric tensor. A simple six tetrahedral elementmesh was adapted to several different metric fields, but the local optimization and node relocation did nothave a effective criteria for optimization based on element swap and node shifting. Carlo L. Bottasso [8]describes a Gauss–Seidel algorithm for optimizing a three-dimensional unstructured grid to conform to agiven metric. He analyzed different possible choices for the objective function, based on the maximum valueof an element residual measuring the distance of any simplex in the grid to the local target. Weizhang Huang[25, 26] developed a general formula for metric tensor to be used in any spatial dimension, based on errorestimates for polynomial preserving interpolation on simplicial elements. The numerical results presentedshow that the defined metric tensor is able to produce an anisotropic mesh with overall good quality whenused with an existing meshing strategy. A more recent scheme was developed by D. Marcum and F. Alauzet[42], which uses advancing-front type point placement that aligns the elements with the underlying metricfield. Even though they have accomplished a considerable alignment, a metric field background mesh witha certain refinement is essential, and their scheme essentially generates a new mesh at every adaptationcycle. For robustness reasons, most strategies rely on local mesh modification operators such as swappingfaces, deleting points, adding points, or smoothing [4, 8, 21, 37, 57, 68]. Pagnutti and Ollivier-Gooch[57] proposed a method for anisotropic Delaunay mesh refinement to higher-order numerical solutions,this was accomplished by assigning metrics to vertices that approximate the error in a region. Numericalresults for subsonic viscous and transonic inviscid flow around an airfoil showed improvement during theadaptation process. Sharbatdar and Ollivier-Gooch [68] extend the approach by adding node-shifting basedwith Target Matrix Optimization Paradigm [32, 33, 34] to aligned the mesh into a quasi structured mesh.Numerical results from subsonic viscous and transonic inviscid flows showed an improvement comparedwith the previous version.1.2 Objective and OutlineIn general, metric based anisotropic schemes can produce highly anisotropic adaptation, but a better elementalignment based the eigen-analysis of the provided metric field is still desired. The Delaunay anisotropicadaptation scheme developed by Pagnutti and Ollivier-Gooch [57] captured a certain anisotropy in a sub-sonic viscous flow and a transonic inviscid flow, but it lacked of any cell alignment with the flow features.Sharbatdar and Ollivier-Gooch [68] extended this scheme by adding a anisotropic adding node-shiftingbased with Target Matrix Optimization Paradigm, in order to achieve a quasi structured mesh. A compar-ison between this schemes is shown for the subsonic viscous flow around an airfoil in Fig. 1.1, a close upview of the top and bottom surface is shown in Fig. 1.2 and a close up view of the trailing edge is shownin Fig. 1.3. This case shows physical behavior that commonly occurs in aerodynamics, an anisotropic fea-ture parallel to the flow within the boundary layer and wake. A certain level of alignment is obtained afteradding the smoothing phase to the algorithm, and a comparison of numerical results presented in [68] showan improvement in efficiency and accuracy after adding the smoothing phase to the algorithm, but a better31.2. Objective and OutlineSmoothedNon-smoothedFigure 1.1: Smoothed and non smoothed for a subsonic viscous flow over a NACA 0012, Ma = 0.5, Re =5000,α = 0◦. Figures courtesy of M. Sharbatdar [67]SmoothedNon-smoothed(a) Close view of top surface with smoothingSmoothedNon-smoothed(b) Close view bottom surface without smoothingFigure 1.2: Close up view of top and bottom surface. Figures courtesy of M. Sharbatdar [67]alignment is still desired.The goal of my research is to create a scheme that adaptively produces anisotropic metric-orthogonalmeshes accordingly to a metric field. Sharbatdar’s introduction of metric-based smoothing [67] improvedPagnutti’s adaptation scheme [56], but Pagnutti’s Delaunay-based point insertion does not contribute tomesh alignment. By creating a new algorithm for inserting a point in a mesh which enhances cell alignment,the efficiency and accuracy of the solution should be improved. This is done by matching the anisotropicelements in the mesh with the anisotropic features of the flow properties, i.e. boundary layers in viscousflows . The effectiveness of the adaptation scheme in producing aligned meshes will not only rely onsmoothing, thus a higher cell alignment should be obtained.The organization of this thesis is as follows: In Chapter 2, we recall some background general terms usedin the process of this research. Flow domains need to be split into smaller sub domains in order to analyzefluid flows so the governing equations can be discretized and solved inside each of these sub domains, whichis called a mesh. In this chapter gives a brief definition of different mesh types based on:• element classification;41.2. Objective and OutlineSmoothedNon-smoothed(a) Close view of trailing edge with smoothingSmoothedNon-smoothed(b) Close view of trailing edge without smoothingFigure 1.3: Close up view of trailing edge. Figures courtesy of M. Sharbatdar [67]• connectivity, such as structured mesh, quasi structured mesh and unstructured mesh; and• aspect ratio of the elements, such as isotropic and anisotropic.To improve the efficiency and accuracy of a numerical solution, a mesh can be adapted in a specific way toresemble the features of the subject being analyzed. One common approach of adapting a mesh is by localmesh operations, such as swapping, coarsening, smoothing and insertion of points. A general definition ofthese operations is given in this chapter. Much of the research in this thesis uses quality measured based onthe Target-Matrix Optimization Paradigm (TMOP), which is described in Sec. 2.4.In Chapter 3, the methodology used in this research is described. The procedure to adaptively produceanisotropic metric-orthogonal meshes based on local mesh operations designed to conform a metric, suchas swapping, insertion, smoothing and coarsening is described. The criteria for each operation is describedas well as the general algorithm. An advantage of this algorithm is the ability of facilitating the merge oftriangular cells into quad cells due to the alignment obtained. A simple process for merging triangular cellsinto quad cells is explained. The process is explained for 2D and then extended to 3D.In Chapter 4, the results obtained from the proposed algorithm will be shown. Results are shown forexamples with analytical metrics which are characterized by a smooth metric field in 2D / 3D and for a 2Dflow around an airfoil with anisotropic features. This to prove the algorithm’s adaptation efficiency and toshow the effectiveness of the adaptation scheme in producing aligned meshes. Mesh quality is computedwith quality metrics based on TMOP and numerical comparisons are shown with other adaptation schemes.In Chapter 5 the thesis concludes with a summary of the research, describing contributions, discussionof the results and recommending some future work.5Chapter 2BackgroundThis chapter gives an overview for the basic concepts required for the proposed adaptation algorithm.Section 2.1 briefly discusses different types of mesh classifications based on type of elements, connectivity,and aspect ratio. Section 2.2 describes the essential local mesh modification operators, in addition to thefundamental approach of metric adaptation used in this research. Section 2.4 introduces the Target-MatrixOptimization Paradigm to define quality, which is also used for mesh optimization via node movement.2.1 Mesh ClassificationThere exist different ways of classifying a mesh depending on its characteristics. An element-based classifi-cation is categorized according to their dimensionality and choice of elements. Meshes are also categorizedaccording to the elements’ connectivity and topology, which include structured, unstructured, and quasi-structured meshes. Another way of classifying a mesh depends on the size, shape and orientation of meshelements throughout the domain. In the case of triangular elements, an equilateral shape would imply anisotropic mesh while large aspect ratio mesh elements would imply an anisotropic mesh (where aspect ratiois the ratio of the longest edge to shortest edge).2.1.1 Element-Based ClassificationMeshes can be classified according to the shapes of their elements. The most common meshes for 2D consistof triangular or quadrilateral element meshes. An example of these meshes is shown in Fig. 2.1. Triangularelements are the simplest polygons having three vertices and three faces for which often a 60◦ angle isdesired. Quadrilateral elements are four-sided polygons. Their sides do not need to be parallel; however,convexity is often required. In quadrilateral elements an angle of 90◦ often is desired.(a) Triangular element and mesh (b) Quadrilateral element and meshFigure 2.1: Classification of common meshes in 2D by elements62.1. Mesh Classification(a) Tetrahedral element and mesh (b) Hexahedral element and meshFigure 2.2: Classification of common meshes in 3D by elements(a) Triangular and quadrilateral mesh (b) Triangular mesh with quadrilateral boundariesFigure 2.3: Mixed element meshesFor 3D, common mesh elements include tetrahedral and hexahedral elements. An example is shown inFig. 2.2. Tetrahedral elements are the simplest of all the polyhedra, having four vertices and four triangularsurfaces. Hexahedral elements are polyhedra with six quadrilateral faces; their faces do not need to be par-allel or even planar. All-hexahedral meshes can present difficulties for complicated domains for topologicalreasons. Other type of meshes can be created by combining different elements; examples of mixed mesheswith triangular and quadrilateral elements are shown in Fig. 2.3.2.1.2 Structured and Unstructured MeshesIn the context of CFD, meshes can be categorized based on their connectivity as structured and unstructured.In a structured mesh every cell and vertex has the same connectivity and topology. A cell (i, j) is locatedbetween cells (i− 1, j),(i+ 1, j),(i, j− 1), and (i, j+ 1) as shown in Fig. 2.4a. An unstructured mesh isidentified by irregular connectivity as shown in Fig.2.4b. Connectivity data must explicitly store neighboring72.1. Mesh Classificationcells in computer memory for a given vertex as described in [22]. Structured meshes are generally composedof quadrilaterals in 2D and hexahedra in 3D, while unstructured meshes are commonly composed of oneor more types of elements: triangles or quadrilaterals in 2D and tetrahedra, prisms or hexahedra in 3D.Examples of these types of meshes were shown in Sec. 2.1.1. [𝑖 + 1, 𝑗] [𝑖, 𝑗] [𝑖 − 1, 𝑗] [𝑖, 𝑗 + 1] [𝑖, 𝑗 − 1] (a) Structured mesh 1 6 5 7 10 11 2 8 3 9 4 (b) Unstructured meshFigure 2.4: Structured mesh and unstructured meshSome points of comparison between structured and unstructured meshes are:• Accuracy:– For simpler geometries, structured meshes can be more accurate than unstructured meshes.– For complex geometries, more accurate numerical solutions can be obtained with unstructuredmeshes due to its easier adaptive capability.– Meshing is easier for complex geometries with unstructured meshes than with structured meshes.• Memory:– Structured meshes use less memory for data structures due to their predicable topology– Unstructured meshes need more storage due to their irregular connectivity.– For unstructured meshes it is common to capture the entire mesh with fewer elements.• CPU time:– The neighboring states must be known for flux computations for the finite volume method. Forstructured meshes this is straight forward, hence less CPU time is needed to compute the solu-tion.82.1. Mesh Classification– Complex problems have been proved to be a considerable challenge for structured meshes, caus-ing often an over refinement of the mesh or time-consuming human interventions.– Unstructured meshes require a polynomial reconstruction of the solution over each control vol-ume, which requires more CPU time. Even though it may be more costly for the solver per meshelement, the ability of unstructured meshes to capture complex geometries with considerablyfewer elements may result in a significant CPU time advantage [47].– From a performance point of view, for unstructured meshes there are added difficulties in un-derstanding both the processing time on a single processor and the scaling characteristics whenusing large-scale parallel systems [43].2.1.3 Quasi Structured MeshesA quasi structured mesh can be defined as a triangular mesh which resembles a quadrilateral mesh withdiagonals added in each quad cell, thus the topology is not predictable and it must be treated as an unstruc-tured mesh. A comparison between a unstructured mesh, a quasi structured mesh and a structured meshis shown in Fig. 2.5. A quasi structured mesh takes advantage of the positive properties of both structuredmesh elements and unstructured mesh elements, for which a high grid quality can be achieved throughoutthe domain [74].(a) Unstructured mesh(b) Quasi structured mesh(c) Structured meshFigure 2.5: Comparison between unstructured, quasi structured, and structured meshes2.1.4 Isotropic and Anisotropic MeshesA mesh can be categorized as isotropic or anisotropic depending on the aspect ratio of the elements in themesh, regardless of the element topology. For the purpose of the terms in this section, we will refer to92.1. Mesh Classification(a) Isotropic mesh (b) Anisotropic meshFigure 2.6: Comparison between isotropic and anisotropic meshes in 2Dtriangular meshes in 2D and tetrahedral meshes in 3D. The ratio of the longest edge to the shortest edgedefines aspect ratio; for isotropic meshes an aspect ratio of near one is desired, meaning almost equilateraltriangles or equilateral tetrahedra; meanwhile higher aspect ratio denotes an anisotropic mesh. A comparisonin 2D is shown in Fig. 2.6.In 3D, tetrahedral elements can behave in different ways. A taxonomy [70] of different types of de-generate tetrahedra is shown in Fig. 2.7 and Fig. 2.8. There are two general types of degenerate tetrahedra,tetrahedra collapsing to an edge shown in Fig. 2.7a to Fig. 2.7e and tetrahedra collapsing to a plane, whichhave a dihedral angle (angle between two planes) near 180◦ or a large circumsphere (unique sphere thatpasses through each of a tetrahedron’s four vertices) shown in Fig. 2.8a to Fig. 2.8d. Even though they areall tetrahedra, they have different characteristics: spears, spindles, spades, caps, and slivers have a dihedralangle near 180◦. Most of the degenerate tetrahedra have a dihedral angle near 0◦. Spires can have all theirdihedral angles between 60◦ and 90◦, with two edges separated by an face angle near 0◦. A more detaileddescription of the characteristics of different types of tetrahedra is given in [70]. The approach for isotropicmeshes in 3D is analogous to 2D in which equilateral tetrahedra are desired. For anisotropic meshes, thedifferent types of tetrahedra mentioned before might be acceptable.The type of mesh required depends on the problem being analyzed. Whenever the physical behav-ior of the problem is isotropic, an isotropic mesh would be preferable if the size of the cells equally dis-tributes the error per cell [56]. However, many physical phenomena of interest in CFD have regions withhighly anisotropic features such as shock waves, boundary layers, and ignition fronts, among others. Foranisotropic features a higher density mesh is required, for which strong gradients perpendicular to the fea-ture are accurately captured but the simulations might be inefficient as a consequence of over-resolved flowin the tangential direction [22]. Anisotropic meshes can resolve these features better, reducing the numberof mesh vertices. This reduction in vertices could achieve a significant reduction in computational cost forthe same accuracy as an isometric mesh.102.2. Mesh Adaptationa) b) c) d) e) (a) spire a) b c) d) e) (b) speara) b) c) d) e) (c) spindlea) b c) d) e) (d) spikea) b c) d) e) (e) splinterFigure 2.7: Classification of skinny tetrahedra collapsing to an edgeg) h) i) f) (a) wedge g) h) i) f) (b) sliver) ) i) f) (c) cap (d) spadeFigure 2.8: Classification of skinny tetrahedra collapsing to a plane2.2 Mesh AdaptationMesh adaptation strategies can usually be classified as one of two general types: moving points from oneregion to another (r-adaptation) and selectively refining/coarsening cells (h-adaptation). For robustness rea-sons, most schemes rely on local mesh modification operators such as swapping faces, deleting points,adding points, and smoothing. This section describes the general idea of local mesh operations used foradapting unstructured meshes. The description of these operations within the proposed algorithm and themetric definitions used in this research are described in Chapter 3.2.2.1 SwappingIn 2D, face swapping is a simple operation which chooses the best diagonal for the quadrilateral formedbetween two neighboring triangles. Let two neighboring triangles be defined by acb and acd; after swappingthe configuration will be abd and bcd as shown in Fig. 2.9a and Fig. 2.9b. For face swapping to be possible,the final triangulation has to fill the convex hull1 of its vertices. Examples of cases which are not suitablefor swapping are shown in Fig. 2.10a and Fig. 2.10b. To decide wether a face swap improves the mesh, ageometric predicate is generally used; more details are given in Sec. 3.1.1.1The convex hull of a set S is the smallest convex polygon that contains all the points of S112.2. Mesh Adaptationa b c d (a) Before swappinga b c d (b) After swappingFigure 2.9: Configurations suitable for swappinga b c (a) Three linear verticesa b c (b) Non-convexFigure 2.10: Configurations not suitable for swappinga) b) a b c d e a b c d e (a) Configuration with 2 tetrahedra) b) a b c d e a b c d e (b) Configuration with 3 tetrahedraFigure 2.11: Configurations for a 2-3 / 3-2 swapping122.2. Mesh Adaptationa) b) a b c d e f g a b c e f g d (a) Configuration with 5 tetrahedraa) b) a b c d e f g a b c e f g d (b) Configuration with 6 tetrahedraFigure 2.12: Configurations for a 5-6 / 6-5 swappingFor 3D, both face and edge swapping are possible. Each interior face in a tetrahedral mesh is adjacentto two tetrahedra. Face swapping reconnects the tetrahedra constructed from these five vertices. Cases forwhich the tetrahedra can be acceptably reconnected are 3-2 and 2-3 swaps, as shown in Fig. 2.11. Let twoneighboring tetrahedra be defined by abcd and acbe. After face swapping its common face abc, the twotetrahedra become three tetrahedra with configuration abed, bced, and acbe. The simplest case for edgeswapping is the inverse of this operation, where three tetrahedra abed, bced, and caed swap a shared edgeto become two tetrahedra abcd and acbe.More complex swapping operations are possible, Fig.2.12 shows a 5-6 / 6-5 edge swap. The 5 tetrahedrain Fig. 2.12a all have a common edge. By swapping this edge we create three new faces abe, bce, and cde,creating six cells as shown in Fig. 2.12b. As in the previous case, it is possible to swap backwards. Ingeneral, edge swapping reconfigures N tetrahedra incident on an edge of the mesh by removing the edgeand replacing the N tetrahedra with 2N−4 tetrahedra as explained in [18, 36].a) b) a b c d e f a b c d e f (a) Configuration before swappinga) b) a b c d e f a b c d e f (b) Configuration after swappingFigure 2.13: Configurations for a 4-4 swapping132.2. Mesh Adaptationa) b) a b c d e a b c d e (a) Three linear verticesa) b) a b c d e a b c d e (b) Non-convexFigure 2.14: Configurations not suitable for swappingFor internal boundaries, a 4-4 swapping is possible. As shown in Fig. 2.13 ace, afc, abc, and acd areswapped for bde, bfd, bcd, and bda. For external boundaries a 2-2 swap is possible following the samecriteria as for internal boundaries. As in 2D, in 3D swapping may not be possible in some cases as shown inFig.2.14. Vertices aeb shown in Fig.2.14a are linear. If swapping occurs the final configuration will containa degenerate tetrahedron. If the configuration of any two tetrahedra is not convex as in Fig. 2.14b, swappingis not possible. To decide whether a face swap or an edge swap improves the mesh, a geometric predicate isgenerally used; more details are given in Sec. 3.4.1.2.2.2 InsertionThis operation adds vertices to a mesh; Fig. 2.15 shows the possibles cases for 2D. Let two triangles bedefined by vertices abc and bad, and the new vertex to be added be defined as A. Cells aAc and Abc will becreated if we split a boundary (Fig. 2.15a). Adding a point anywhere inside the cell creates three new cellsabA, bcA, and caA (Fig. 2.15b). If A splits an edge from an interior cell four new cells adA, dbA, bcA andcaA are created (Fig. 2.15c).a b c A (a) Insertion on a boundaryfacea b c A (b) Insertion inside a celld a b c d a b c A (c) Insertion on an interior faceFigure 2.15: Insertion of a vertex in 2D142.2. Mesh AdaptationFor 3D, possible insertion cases are shown in Fig. 2.16. Let a tetrahedron be defined as abcd and thevertex to insert be defined as A. Inserting a vertex on an edge splits the tetrahedron into two tetrahedra aAcdand bcAd (Fig. 2.16a). If a vertex is inserted on a face three tetrahedra abcA, acdA and bdcA are created(Fig. 2.16b). If a vertex is inserted inside the cell, the tetrahedron splits into four tetrahedra abcA, acdA,bdcA, and abdA (Fig. 2.16c).a) b) a b c d a b c d A a) b) a b c d a b c d A (a) Insertion on an edgec) d) a b c d a b c d A A c) d) a b c d a b c d A A (b) Insertion on a facec) d) a b c d a b c d A A c) d) a b c d a b c d A A (c) Insertion inside a cellFigure 2.16: Insertion of a vertex in 3DThe previous description for insertion on a face in 3D also applies for neighboring tetrahedra. Splittinga face turns two incident tetrahedra into six tetrahedra and splitting an edge turns n incident tetrahedra into2n tetrahedra. An example of edge splitting with four incident tetrahedra is shown in Fig. 2.17.a) b) a b c d e f c) d) a b c d e f A (a) Original configuration c) d) a b c d e f A c) d) a b c d e f A (b) Configuration after insertionFigure 2.17: Edge insertion: 4 to 8 tetrahedra152.2. Mesh Adaptation2.2.3 CoarseningCoarsening removes a vertex from the mesh and reconfigures the neighboring cells. Figure2.18 shows somepossible cases for coarsening in 2D. Removing a vertex is simple if it has only three cells incident, as shownin Fig.2.18a. Removing vertex d will create a new cell abc. The same case applies for boundaries, as shownin Fig. 2.18c. If the vertex to be removed has more incident cells as shown in Fig. 2.18b, deleting the vertexwill delete all its neighboring cells so re-triangulating the hull will be required. Edge contraction is one easyway to do this. By contracting edge ae from Fig. 2.18b until both points overlap, cells aed and abe will bedeleted, thus leaving only cells abc and acd which end up stretched as shown.a b c d e J a b c a b c d a b c d a b c d a b c (a) From three incident cellsa b c d e J a b c a b c d a b c d a b c d a b c (b) From four incident cellsa b c d e J a b c a b c d a b c d a b c d a b c (c) Boundary vertexFigure 2.18: Deleting a vertex in 2Dc) d) a b c d e f A (a) Vertex with eight adjacent cells a) b) a b c d e f (b) After coarseningFigure 2.19: Deleting a vertex from interior cells in 3D162.2. Mesh Adaptationc) d) a b c d a b c d A A (a) Vertex with four adjacent cellsc) d) a b c d a b c d A A (b) Boundary vertex with 3 adjacent cellsa) b) a b c d a b c d A (c) Boundary vertex with two adjacent cellsa) b) a b c d a b c d A (d) After coarseningFigure 2.20: Deleting a vertex from a boundary cells in 3DAs in 2D, in 3D edge contraction is one easy way of doing this. Examples of edge contraction are shownin Fig.2.19 and Fig.2.20. By contracting edge eA from Fig.2.19a until both points overlap, cells abAe, bcAe,cdAe, and daAe will be deleted, thus leaving only cells adAf, dcAf, cbAf, and baAf which end up stretchedas shown in Fig. 2.19b. Further details on coarsening unstructured meshes by edge contraction are given in[54].2.2.4 SmoothingSmoothing repositions vertices to improve the quality of the mesh, as shown in Fig. 2.21. Let four trianglesbe defined as abe, bce, cde, and dae, in which vertex e is incident to all the cells. If we consider anisotropic quality in which an equilateral triangle is desired, repositioning vertex e would increase the qualityof the adjacent cells. The procedure for 3D is analogous to the description given for 2D. Figure 2.22 showsan example of smoothing an interior vertex adjacent to eight tetrahedra. An increase in quality in all thetetrahedra is obtained by smoothing the vertex. There are several smoothing techniques that can be appliedfor unstructured meshes using optimization techniques; more details for smoothing can be found in [16, 17].172.2. Mesh Adaptationa) b) a b c d e a b c d e (a) Initial configurationa) b) a b c d e a b c d e (b) After smoothingFigure 2.21: Smoothing a vertex in 2D(a) Exterior view of eight tetrahedra (b) Configuration before(c) Configuration afterFigure 2.22: Smoothing a vertex in 3D182.3. Metric-Based Adaptation2.3 Metric-Based AdaptationA metric is a function Rn×Rn→ R that satisfies the following conditions for all x,y,z:d(x,y)> 0⇐⇒ x 6= y (2.1)d(x,y) = d(y,x) (2.2)d(x,z)≤ d(x,y)+d(y,z) (2.3)A Riemannian metric space (Rn,M) in n dimensions is a vector space of finite dimension where the dotproduct is defined by a symmetric definite positive tensor M called a metric tensor.u ·Mv = 〈u,v〉M = 〈u,Mv〉= utMv, ∀(u,v) ∈ Rn×Rn (2.4)The standard Euclidean distance is a metric for which M = I, where I is an identity matrix.Given a metric tensor Mp, a deformation tensor Fp can be defined as any d×d matrix satisfying:Mp = FTp Fp and detFp > 0 (2.5)where Fp maps the physical space into a metric space as in Fig. 2.23 in which distances and angles aremeasured as seen from p’s point of view [35].Figure 2.23: Deformation of tensors Fp and FqThe ellipses around the points p and q represent its isocontours which are equidistant from each otherin the metric space. The ellipses specify three major components to define anisotropy: two perpendicularradius r1 and r2 and an angle θ as shown in Fig. 2.24; more details are given in [1, 35, 56, 68] . The metric192.3. Metric-Based Adaptationis defined as a function of these three components:M =ΘΓΘT =(cosθ −sinθsinθ cosθ) 1r21 00 1r22( cosθ sinθ−sinθ cosθ)(2.6)This sets the desired cell size, shape and orientation at each point in the domain.Figure 2.24: Distance and angle measured for a certain pointThe main idea of metric-based adaptation is to generate a unit mesh in a prescribed Riemmanian metricspace. It is possible to work in a Riemmanian metric space as the notions of length, volume and angle arewell-posed. The metric tensor is commonly computed from an error-estimation method, which defines thedesired anisotropy, often generated from the Hessian of a numerically computed solution. George, Hechtand Vallet [20] introduced the use of metrics in a Delaunay mesh generator. They showed that the absolutevalue of the Hessian is a metric and proposed a Delaunay-based mesh generator by computing edge lengthin the metric space. The idea of using the Hessian as a metric was popularized by results of D’Azevedo andSimpson [14] on the linear interpolation error in a function and its gradient. Since then, the advantageousidea of metrics has been widely exploited for anisotropic mesh adaptation [7, 40, 41, 42, 56, 67].2.3.1 Metric OperationsA metric is set at each vertex of the mesh but we require a metric per cell and a metric per edge by the metricalgorithm. For this reason a relation between them must be defined. The most mathematically rigorous wayof doing this is based on the Log-Euclidean framework described in [1, 2, 3, 42].Since the metric tensor M is symmetric, it is diagonalizable:M = RΛRT (2.7)where R is a unitary matrix of eigenvectors (~Vi)i=1,...,n and Λ = diag(λi)i=1,...,n is a diagonal matrix of itseigenvalues. For a Log-Euclidean framework, we will define the notion of a metric logarithm and exponen-tial by:log(M) := R log(Λ)RT (2.8)202.3. Metric-Based Adaptationexp(M) := Rexp(Λ)RT (2.9)where ln(Λ) = diag(ln(λi)) and exp(Λ) = diag(exp(λi)). Define logarithmic addition and logarithmic scalarmultiplication as:M1⊕M2 := exp(log(M1)+ log(M2)) (2.10)αM := exp(α · log(M)) (2.11)Let (xi)i=1...k be a set of vertices in k−1 dimensions and (M(xi))i=1...k their corresponding metric tensors.For a point x of the domain such that x =k∑i=1αixi with ∑ki=1αi = 1 the interpolated metric in the Log-Euclidean framework is defined by:M(x) =⊕ki=1αiM(xi) = exp(k∑i=1αi log(M(xi)))(2.12)2.3.2 Metric LengthLet E be an edge of the mesh formed from two vertices Vi=1,2 in Euclidean space, with M(V1) and M(V2)their respective metrics. The metric length will be denoted:lMi(E) =√EtM(Vi)E (2.13)If we assume lM1(E)> lM2(E) and set a =lM1 (E)lM2 (E), we can define the metric length of E as:lM(E) = lM1(E)a−1a log(a)(2.14)Considering other linear variations along the edge, other similar approximations are analyzed in [1].2.3.3 Metric per CellA metric is also required per cell. Since the metric is set in each vertex of a cell, a relation between themmust be defined. Based on the Log-Euclidean framework previously described we will define the metric fora cell as:M(C) = exp(1kk∑i=1log(M(Vi)))(2.15)Where k is the number of vertices in the cell and M(Vi) the respective metric of vertex Vi.A more complete description of metric operations can be found in [42, 3, 1, 2]212.4. Target Matrix Optimization Paradigm2.4 Target Matrix Optimization ParadigmThe Target-Matrix Optimization Paradigm (TMOP) is an approach for mesh optimization via node-movementwhich can be applied to a wide variety of mesh and element types as shown in [32, 33, 34, 31]. TMOP isbased around the affine mappings that relates a reference element to another element through for local map-ping, a target element which describes the element desired in the mesh, and an actual element which is a cellin the mesh as shown in Fig. 2.25. Actual Element Target ElementReference ElementT k=A kW k−1Ak W kFigure 2.25: Weighted Jacobian matrix definitionLet ζ = (ξ ,η) and κ = {ζ | ξ ≥ 0,η ≥ 0,ξ +η ≤ 1} be a reference element triangle shown in Fig.2.25and x0, x1, and x2 be three points of a right handed cell in the mesh. For linear triangles in the xy-plane, themapping from the reference element to the actual element is:x(ζ ) = x0+(x1− x0)ξ +(x2− x0)η (2.16)Then xξ = x1− x0 and xη = x2− x0, so the Jacobian matrix Ak =[xξ ,xη].The Jacobian matrix Ak derived from the actual mesh relates the reference element with an actual elementin the mesh. The Target matrix Wk is mapping between the reference element and the target element, whilethe weighted Jacobian matrix Tk is a relation between the Jacobian matrix Ak and the Target matrix Wk. TheTarget matrix Wk is defined for every element the mesh. Since the construction of targets is user defined,det(W ) 6= 0 and W−1 exists. The weighted non-dimensional matrix Tk = AkW−1k provides a convenientscaling of Ak. The objective function that is optimized is defined as:F =1N∑kϒ(Tk) (2.17)where k is the element index, N the total number of elements in the mesh, and ϒ is a quality metric that is afunction of the weighted Jacobian matrix. The objective function is minimized by moving the free verticesin the mesh.Using the QR factorization of A we can identify factors with geometric meanings:A = OU = ΞOS = ΞOΦ∆ (2.18)where Ξ is a scalar local size parameter, O is related to local orientation, Φ is related to local angle or skew,∆ is related to local aspect ratio, S is related to local shape, and U is related to local shape + size. Setting ai222.4. Target Matrix Optimization Paradigmas vectors from the columns of A and α = det(A) 6= 0, these factors are defined in 2D/3D by:Ξ(A)2D = |α|1/2 Ξ(A)3D = |α|1/3 (2.19)O(A)2D =[a1|a1| ,−(a1 ·a2)a1+ |a1|2 a2|α| |a1|]O(A)3D =[a1|a1| ,−(a1 ·a2)a1+ |a1|2 a2|a1| |a1×a2| ,α (a1×a2)|α| |a1×a2|](2.20)Φ(A)2D =√|a1| |a2|α(1 a1·a2|a1||a2|0 |α||a1||a2|)Φ(A)3D =3√|a1| |a2| |a3|α1 a1·a2|a1||a2|a1·a3|a1||a3|0 |a1×a2||a1||a2|(a1×a2)·(a1×a3)|a1×a2||a1||a3|0 0 |α||a1×a2||a3|(2.21)∆(A)2D =1√|a1| |a2|(|a1| 00 |a2|)∆(A)3D =13√|a1| |a2| |a3| |a1| 0 00 |a2| 00 0 |a3| (2.22)2.4.1 Target MatrixTo define the target matrix for an element, we set a metric per element as mentioned in Sec. 2.3.3 and findthe eigen decomposition of the metric. Let ~Vi=1,...,d be the normalized eigenvectors of a d× d matrix andλi=1,...,d its eigenvalues; then the target matrix in 2D is defined by:W 2D =(~V1,~V2)( 1√λ100 1√λ2)(2.23)Two examples of 2D target matrices are shown in Eq. 2.24, one isotropic and the other anisotropic. Bysetting the columns of the target matrix as edge vectors of the target element, we can define our targetelement as shown in Fig. 2.26.W 2Diso =(1 120√32)W 2Dani =(10 00 1)(2.24)Following the same approach as in 2D, we can define the target matrix for 3D as:W 3D =(~V1, ~V2, ~V3)1√λ10 00 1√λ2 00 0 1√λ3 (2.25)232.4. Target Matrix Optimization Paradigm𝑥1 𝑥2 𝑥0 𝑉1 𝑉0 𝑥1 𝑥2 𝑥0 𝑉1 𝑉0 a) b) 𝑥0 = 0 , 0 𝑥1 = 1 , 0 𝑥2 = 12 ,32 𝑉1 = 12 ,32 𝑉0 = ( 1 , 0 ) 𝑥0 = 0 , 0 𝑥1 = 10 , 0 𝑥2 = 0 , 1 𝑉1 = 0 , 1 𝑉0 = ( 10 , 0 ) (a) For an isotropic element𝑥1 𝑥2 𝑥0 𝑉1 𝑉0 𝑥1 𝑥2 𝑥0 𝑉1 𝑉0 a) b) 𝑥0 = 0 , 0 𝑥1 = 1 , 0 𝑥2 = 1 ,32 𝑉1 = 12,32 𝑉0 = ( 1 , 0 ) 𝑥0 = 0 , 0 𝑥1 = 10 , 0 𝑥2 = 0 , 1 𝑉1 = 0 , 1 𝑉0 = ( 10 , 0 ) (b) For an anisotropic elementFigure 2.26: Target matrix in 2D𝑥1 𝑥0 𝑉0 = 2 , 2 , 0 a) b) 𝑥2 𝑥3 = 2 , 0 , 2 𝑥0 = 0 , 0 , 0 𝑥1 = 2 , 2 , 0 𝑥2 = 0 , 2 , 2 𝑥3 𝑉0 𝑉2 𝑉1 𝑉1 = 0 , 2 , 2 𝑉2 = 2 , 0 , 2 𝑥1 𝑥0 𝑥2 𝑥3 𝑉0 𝑉2 𝑉1 𝑉0 = 30 , 30 , 0 𝑥3 = 1 , 0 , 1 𝑥0 = 0 , 0 , 0 𝑥1 = 30,30 , 0 𝑥2 = 0 , 10 , 10 𝑉1 = 0 , 10, 10 𝑉2 = 1 , 0 , 1 (a) For an isotropic element𝑥1 𝑥0 𝑉0 = 2 , , 0 a) b) 𝑥2 𝑥3 = 2 , 0 , 2 𝑥0 = 0 , , 𝑥1 = 2 , , 0 𝑥2 = 0 , 2 , 𝑥3 𝑉0 𝑉2 𝑉1 𝑉1 = 0 , 2 , 𝑉2 = 2 , 0 , 2 𝑥1 𝑥0 𝑥2 𝑥3 𝑉0 𝑉2 𝑉1 𝑉0 = 30 , 30 , 𝑥3 = 1 , 0 , 1 𝑥0 = 0 , , 𝑥1 = 30,30 , 𝑥2 = 0 , 10 , 10 𝑉1 = 0 , 10, 10 𝑉2 = 1 , 0 , 1 (b) For an anisotropic elementFigure 2.27: Target matrix in 3D242.4. Target Matrix Optimization ParadigmTwo examples of 3D target matrices are shown in Eq. 2.26, one isotropic and the other anisotropic. Bysetting the columns of the matrix as edge vectors we can define a target element, some examples are shownin Fig. 2.27.W 3DIso = 2 0 22 2 00 2 2 W 3DAni = 30 0 130 10 00 10 1 (2.26)For a given metric per cell in 2D:M =(mxx mxymxy myy)(2.27)A tolerance for defining an isotropic target matrix for such cell is set as:|mxx−myy|< 10−6,−10−6 < mxy < 10−6 (2.28)While for a given metric per cell in 3D:M = mxx mxy mxzmxy myy myzmxz myz mzz (2.29)A tolerance for defining an isotropic target matrix for such cell is set as:|mxx−myy||myy−mzz||mzz−mxx|< 10−6,−10−6 <mxymxzmyz< 10−6 (2.30)2.4.2 Orientation and OrderOnce a target is computed as mentioned in Sec.2.4.1, it can have different orientations. Let five points in 2Dbe defined as a, b, c, d, and e. The possible configurations of the target matrix depending on the orientationare as shown in Fig. 2.28b:W 2D1 =(bx−ax cx−axby−ay cy−ay)W 2D2 =(dx−ax ex−axdy−ay ey−ay)(2.31)W 2D3 =(cx−ax dx−axcy−ay dy−ay)W 2D4 =(ex−ax bx−axey−ay by−ay)The weighted Jacobian Tk defined in Fig.2.25 is a relationship between our target matrix and the Jacobianof the physical element, thus the order of the configuration from our element also must be considered. Thepossible configurations depending on the order for a right-handed triangle element are as shown in Fig.2.28a:252.4. Target Matrix Optimization ParadigmA2D1 =(bx−ax cx−axby−ay cy−ay)(2.32)A2D2 =(cx−bx ax−bxcy−by ay−by)A2D3 =(ax− cx bx− cxay− cy by− cy)Element a b c abc bca cab Order (a) Element orderTarget a b c d e a b c a c d a d e a b e abc acd ade aeb (b) Target orientationFigure 2.28: Possible 2D configurationsElement Order a b c d abcd bcad cabd adbc dbac badc acdb cdab dacb bdca dcba cbda (a) Element ordere b c d f g a Target c c c c d d d d f f f f e e e e b b b b g g g g a a a a a a a a acdf abcf adcg acbg abeg aedg aebf adef (b) Target orientationFigure 2.29: Possible 3D configurations262.4. Target Matrix Optimization ParadigmAs in 2D, the same approach is handled for 3D. Let seven points in 3D be defined as a, b, c, d, e, f , andg, the possible configurations of the target matrix depending on the orientation are as shown in Fig. 2.29b:W 3D1 = dx−ax ex−ax fx−axdy−ay ey−ay fy−aydz−az ez−az fz−az W 3D2 = cx−ax dx−ax fx−axcy−ay dy−ay fy−aycz−az dz−az fz−az (2.33)W 3D3 = ex−ax bx−ax fx−axey−ay by−ay fy−ayez−az bz−az fz−az W 3D4 = bx−ax cx−ax fx−axby−ay cy−ay fy−aybz−az cz−az fz−azW 3D6 = ex−ax dx−ax gx−axey−ay dy−ay gy−ayez−az dz−az gz−az W 3D6 = dx−ax cx−ax gx−axdy−ay cy−ay gy−aydz−az cz−az gz−azW 3D7 = bx−ax ex−ax gx−axby−ay ey−ay gy−aybz−az ez−az gz−az W 3D8 = cx−ax bx−ax gx−axcy−ay by−ay gy−aycz−az bz−az gz−azThe same relation exists in 3D to compute the weighted Jacobian between the target and the Jacobianfrom the physical element. The possible configurations depending on the order for a right-handed tetrahe-dron element are as shown in Fig. 2.29a:A3D1 = bx−ax cx−ax dx−axby−ay cy−ay dy−aybz−az cz−az dz−az A3D2 = dx−ax bx−ax cx−axdy−ay by−ay cy−aydz−az bz−az cz−az A3D3 = cx−ax dx−ax bx−axcy−ay dy−ay by−aycz−az dz−az bz−az (2.34)A3D4 = dx−bx cx−bx ax−bxdy−by cy−by ay−bydz−bz cz−bz az−bz A3D5 = cx−bx ax−bx dx−bxcy−by ay−by dy−bycz−bz az−bz dz−bz A3D6 = ax−bx dx−bx cx−bxay−by dy−by cy−byaz−bz dz−bz cz−bzA3D7 = dx− cx ax− cx bx− cxdy− cy ay− cy by− cydz− cz az− cz bz− cz A3D8 = ax− cx bx− cx dx− cxay− cy by− cy dy− cyaz− cz bz− cz dz− cz A3D9 = bx− cx dx− cx ax− cxby− cy dy− cy ay− cybz− cz dz− cz az− czA3D10 = bx−dx ax−dx cx−dxby−dy ay−dy cy−dybz−dz az−dz cz−dz A3D11 = cx−dx bx−dx ax−dxcy−dy by−dy ay−dycz−dz bz−dz az−dz A3D12 = ax−dx cx−dx bx−dxay−dy cy−dy by−dyaz−dz cz−dz bz−dz272.4. Target Matrix Optimization Paradigm2.4.3 Quality MetricsQuality metrics can be defined with the weighted Jacobian matrix A, which relates the Jacobian matrix Aand the target matrix W . For each cell C in 2D, there are three possible Jacobian (mapping) matrices Ai andfour possible orientations of the target matrix Wj. We evaluate the element quality as:Qϒ (C) = mini, jϒ(AiW−1j)(2.35)where ϒ is a quality metric criterion such as Shape, Size, ShapeOrient, ShapeSize, and ShapeSizeOrient[31, 28, 29, 30]. These metrics are optimized by minimizing.Using a mesh of square domain [10×0],[0×10] with 174 cells and 104 vertices, the quality of everyelement was computed with different quality metrics as shown in Fig.2.30. In all the cases an ideal isotropicequilateral target matrix as in Eq. 2.24 was set for the entire mesh.• Fig. 2.30a: Shape 2D/3DThe quality depends only on the shape of the target, hence the elements with better quality have ahigher resemblance to an equilateral triangle.2D :|T |2F2det(T )−1 3D : |T |3F3√3det(T )−1 (2.36)where the Frobenius matrix norm is defined as |A|F =√∑mi=1∑nj=1∣∣ai j∣∣2 for Am×n.• Fig. 2.30b: SizeThe quality depends only on the size of the target.det(T )+1det(T )−2 (2.37)• Fig. 2.30c: ShapeOrientThe quality depends on a combination of shape and orientation of the target.|T |F − tr(T )√22det(T )(2.38)where tr is the trace of a matrix.• Fig. 2.30d: ShapeSizeThe quality depends on a combination of shape and size of the target.∣∣∣T −T−det(T )∣∣∣2F(2.39)282.4. Target Matrix Optimization Paradigm• Fig. 2.30e: ShapeSizeOrientThe quality depends on a combination of the three criteria: size, shape, and orientation of the target.|T − I|2F2det(T )(2.40)0.100 0.200 0.300 00.351Shape(a) Shape0.100 0.200 0.300 0.400 0.500 1.1e-050.585Size(b) Size0.02500.05000.07500.100 00.124ShapeOrient(c) Shape Orientation0.400 0.800 1.20 0.07221.55ShapeSize(d) Shape Size0.100 0.200 0.300 0.008860.356ShapeSizeOrient(e) Shape Size OrientationFigure 2.30: Quality with different metrics29Chapter 3MethodologyThis chapter describes how the basic operations to modify a mesh such as swapping, insertion, coarseningand smoothing from Sec. 2.2, are combined with metric ideas to produce an algorithm that anisotropicallyadapts a mesh with good alignment. Section 3.1 describes the mesh optimization techniques on whichthe proposed algorithm is based on in 2D. A new criteria for insertion which enhances cell alignment andcoarsening which deleted points without loosing alignment based on metric length and quality are introducedin this section. Swapping is done in metric space (which is essential for the insertion process) and smoothingwhich is defined by target elements. This operations are combined together to form a new algorithm in whicha description of how the mesh is adapted is provided in Sec. 3.2. An advantage of having an aligned quasistructured mesh is the ability of merging cells into quadrilaterals; a simple way of merging triangular cellsinto quad cells is explained in Sec. 3.3. Section 3.4 describes the extended version to 3D of the new 2Dinsertion process, coarsening with a new approach based on a conflict graph in metric space, swapping inmetric space and an extended version of the 2D smoothing. This operations are combined together to forma new algorithm in Sec. 3.5.The metric in each cell defined as M(C), is computed as in Eq. 2.15. The metric length of an edge in acell will be defined as lM(Ei), where i represents the edge number of the cell, computed as in Eq. 2.14. Wewill define our minimum and maximum parameters as lmaxM and lminM , where lM represents metric length.3.1 2D Mesh Optimization TechniquesIn this section the mesh operations described in section Sec. 2.2 such as swapping in metric space, insertionto improve cell alignment, smoothing defined by target elements and coarsening based on metric length andquality are defined to adapt a mesh to conform a metric. The criteria on which these operations are basedare explained with further detail for 2D.3.1.1 Swapping in 2DThe swapping decision is based on minimizing the maximum angle in the metric space for 2D as shown inFig. 3.1. Let two triangles be defined as A and B representing our initial configurations abc and acd. Twoother triangles A′ and B′ represent the configuration after swapping abd and bcd. The top triangles representthe physical space configurations before and after swapping, while the bottom triangles represent the samecells in metric space before and after swapping. The angle between two faces of a cell in metric space isdefined using the law of cosines as:303.1. 2D Mesh Optimization Techniquesθ(Ei,i+1) = arccos(lM(Ei)2+ lM(Ei+1)2− lM(Ei−1)22 · lM(Ei) · lM(Ei+1))(3.1)(a) Before swapping (b) After swappingFigure 3.1: 2D swapping in metric spaceThe two triangles swap their shared face only if:max[(θα ,θβ ,θγ), (θδ ,θε ,θζ )]> max[(θα ′ ,θβ ′ ,θγ ′), (θδ ′ ,θε ′ ,θζ ′)](3.2)3.1.2 Insertion in 2DWe insert points into the mesh with two different goals at different stages in the algorithm using a similarapproach at the boundary as in the interior. Figure 3.2a shows a triangle with boundary edge ab, whileFig. 3.3a shows two interior triangles A and B that share edge ab.In our primary insertion we seek to improve cell shape and orientation. That is, in the boundary case,we insert if worst-case quality is improved by bisecting the edge:Qϒ(A)> max(Qϒ(B), Qϒ(C)) (3.3)while for interior edges, we insert if:313.1. 2D Mesh Optimization Techniquesmax(Qϒ(A), Qϒ(B))> max(Qϒ(C),Qϒ(D),Qϒ(E),Qϒ(F)) (3.4)In a separate pass, we also insert to eliminate long edges in the metric space. For both boundary andinterior edges, our criterion is to split edges where:lmaxM < lM(E) (3.5)To prioritize edges for splitting we use an insertion queue. The position for insertion on the edge iscomputed after all the edges are evaluated and queued. To avoid creating cells and faces smaller thandesired, an edge will not be queued if any new edge created from the insertion of the vertex would be shorterthan lminM .The queue uses a different priority for each of the insertion criteria. The quality insertion priority isbased on quality improvement, where prospective splits that give higher quality are queued with higherpriority. The long edge insertion priority is based on metric length where edges with larger metric lengthget higher priority.a b c a b d A B C c P2 P1 (a) Initial bdry cella b c a b d A B C c P2 P1 (b) Add to queuea b c a b d A B C c P2 P1 (c) InsertionFigure 3.2: Insertion for a long boundary edgeWhen a face reaches the top of the queue, a position along the face for the new vertex must be defined.Consider the boundary triangle A defined by vertices abc as shown in Fig. 3.2a. A vertex assumed to beinserted at the middle of the edge ab is represented by vertex d with a metric as defined in Eq. 2.12. We biasthe location of the new vertex towards a or b depending on which vertex’s metric gives a shorter length foredge ab. We define a limiting position for the new vertex:∆M =lMa (ab)lMa (ab)+ lMb (ab)(3.6)323.1. 2D Mesh Optimization Techniquesd a b c d a b c B A D F E C e P2 P1 a) b) c) (a) Initial int cellsd a b c d a b c B A D F E C e P2 P1 a) b) c) (b) Add to queued a b c d a b c B A D F E C e P2 P1 a) b) c) (c) InsertionFigure 3.3: Insertion for a long interior edgeStarting from the middle of the face, the new vertex position moves along the face a distance 4l =l (E)× 0.05 if Qϒ is improved. Regardless, we do not move the point P2 located ∆M of the way along theedge from a towards b. Since the desired target is a right triangle, only a couple of steps is required in mostcases. For interior edges, exactly the same procedure is followed. After all the points from the queue havebeen inserted, cells affected by the insertion are swapped if necessary to optimize mesh quality.3.1.3 Coarsening in 2DOur goal in removing a vertex is to eliminate edges that are too short in the metric space. A coarseningqueue is used to tag the edges for which lM(E) < lminM . The priority in the queue is set to queue first theedges with the shortest metric length. If a vertex only has one adjacent edge which is smaller than lminM , thatvertex is not marked to be deleted. Vertices which only have one edge smaller than lminM where handled withbetter alignment results with other mesh operations of the algorithm, such as smoothing and swapping. Ifneither of the vertices of the edge has at least two incident edges which are smaller than lminM , the edge isremoved from the queue.Let an edge be defined by two vertices a and b and its neighboring cells A, B, C, D, E, F , and G asshown in Fig. 3.4a. In case an edge has both vertices marked for deletion, the vertex with the worst averagequality Qϒ from its adjacent cells will be deleted.A simple way to delete a vertex is by edge collapse. Vertex a is moved onto vertex b. This way we caneliminate cells C and F along with vertex a, as shown in Fig.3.4b. Local swapping within the modified cellsis done to minimize the maximum angle in the metric space as defined in Sec. 3.1.1.333.2. 2D General AlgorithmA B G C D E B F A B E C D a b f a e b c c f e b c d (a) Add edge to queueA B G C D E B F A B E C D a b f a e b c c f e b c d (b) Re triangulating hullFigure 3.4: Deletion of a vertex3.1.4 Smoothing in 2DFor smoothing, we are using the Target-Matrix Optimization Paradigm (TMOP) as implemented in Mesquite[19], setting a ShapeOrient quality metric. One of our objectives is optimizing for a quasi-structured mesh,so in 2D we are setting the target elements as right triangles as explained in Sec.2.4 using a similar approachas defined in [67]. Due to the capabilities of Mesquite, only interior vertices are smoothed; boundary verticesare not affected by smoothing.3.2 2D General AlgorithmThe general algorithm in 2D to adapt a mesh based on a metric field is shown in Alg. 3.1. A metric mustbe previously defined by an analytical function or an error estimate. The metric is set at each vertex of themesh which will define how the mesh will be adapted. An example of how the algorithm modifies a mesh isshown in Fig. 3.6 and the mesh size during is shown in Table 3.1. This case, which uses an analytic metric,is described in detail in Sec. 4.1.1. A representation of the behavior of the analytic function from Eq. 4.1 isshown in Fig. 3.5. As the cells approach to the center of the anisotropic feature the aspect ratio of the cellsincreases.Once the metric is set, a first insertion of points based on quality is made. All the edges in the meshare evaluated to determine which ones will be split as described in Sec. 3.1.2. Every new cell created in thisprocess is re-evaluated for addition to the insertion queue. This process continues until the insertion queueis empty. All cells modified by the quality insertion process are evaluated for swapping. Since the qualityinsertion process is focused only on shape and orientation, a relatively low number of vertices is added butorientation is improved near areas where high aspect ratio is required producing good mesh alignment evenbefore the smoothing process as seen in Fig. 3.6b.Smoothing is done with a right angle target based on shape and orientation as described in Sec. 2.2.4.Since global smoothing is being used, all the interior vertices in the mesh are affected by the smoothing343.2. 2D General Algorithmprocess, therefore all the faces in the mesh are evaluated for swapping (Fig.3.6c). A second quality insertionstep follows the smoothing and swapping step (Fig. 3.6d).All the previous steps described are based in an eigen-decomposition of the metric which is used toadapt the mesh for shape and orientation. The next two steps improve the mesh size. An insertion based ona maximum metric length is done as described in Sec. 3.1.2. The result of this process is shown in Fig. 3.6e.A coarsening step based on a minimum metric length as described in Sec. 3.1.3 is done. This step willeliminate all the vertices where the level of refinement is more than the metric specifies, while keeping goodalignment due to the quality criteria used within the coarsening process. Next, all the edges are evaluated forswapping to improve the results from coarsening the mesh (Fig.3.6f). The process from the second step untilthe eight step is repeated until the number of added and deleted points in a pass drops below a threshold,which is defined by a percentage of the mesh. The objective of the threshold is to keep iterating until itconverges into a final mesh. A final insertion based on quality is done to reduce any big cells produced fromcoarsening.Figure 3.5: Representation of quarter circle metric function behaviorProcess # Vertices # CellsStep 1.- Set Metric (Initial Mesh) 264 462Step 2.- Quality Insertion 832 1474Step 5.- Quality Insertion 929 1668Step 6.- Long Edge Insertion 2067 3944Step 8.- Swapping after coarsening 1827 3491Table 3.1: Mesh size during 2D adaptation process353.2. 2D General AlgorithmAlgorithm 3.1 2D general algorithmStep 1 A metric previously defined is set in all the vertices Step 2 Insertion based on quality, every face of a new cell is reevaluated for insertion and swapping Step 3 Smoothing with right angle targets based on shape and orientation Step 4 All the faces are evaluated for swapping recursively Step 5 2nd Insertion based on quality, reevaluating new faces for insertion and swapping Step 6 Faces larger than a max previously defined are split Step 7 Faces smaller than a min are deleted, decision for which vertex to delete based on quality Step 8 All the faces are evaluated for swapping recursively Quality Insertion Smooth Long Edge Insertion Coarsen Quality Insertion Swap Queue Swap Swap Set metric Step 9 Total of vertices added/deleted Ω is higher than α keep iterating, otherwise step 2 is repeated only Repeat step 2 Compute metric No Yes Ω >α 363.2. 2D General Algorithm(a) Initial mesh (b) 1st quality insertion (step 2)(c) Smoothing and swapping (step 4) (d) 2nd quality insertion (step 5)(e) Long edge insertion (step 6) (f) Swapping after coarsening (step 8)Figure 3.6: 2D adaptation process 373.3. Tri-to-Quad Conversion3.3 Tri-to-Quad ConversionQuad meshes as shown in Fig.2.1 are preferred over triangular meshes for some application for severalreasons. In some numerical analysis, such as finite element methods with highly elastic and plastic domains,a reduction in approximation error and number of elements can be obtained with quad meshes comparedwith triangular meshes. Many geometries have two dominant local directions, which can be associated withprincipal curvature directions, for which quads can be aligned. More details are reviewed in [6].To convert a triangular mesh into a quad mesh, a simple approach is to pair original adjacent trianglesinto a quad. Several methods for combining triangular cells have been develop:• SQuad [23] is designed to improve internal representation of meshes.• BlossomQuad [64] is an algorithm that finds matches from combinatorial global optimization throughlocal operations on the connectivity.• A greedy approach is presented in [71] where best candidates are identified and the remaining trianglesare brought together through local operations.One of our objectives with the proposed algorithm in this research is to optimize for a quasi-structured mesh,which is a good candidate for combining triangular cells into quad cells. We will only use a greedy algorithmwhich matches the two neighboring cells based on quality without applying any local modifications to showthe effectiveness of this algorithm for creating aligned cells. To combine two triangular cells into a quad,quality based with a ShapeOrient metric is computed using an analogous approach as in Sec. 2.4.3. Thequality for every pair of neighboring cells is computed and stored for every edge. Let a cell be defined byA with neighboring cells B, E, I and a second cell B with neighboring cells A, C, G as shown in Fig. 3.7.For cell A the best merge option is the neighbor producing the best quality (in this case, cell B). If for cellB the best option for merging into a quad is also A, those two cells will merge. In case it is not the bestoption, it will continue to search for best matches. When it there are any best matches left, the cells willstart merging prioritizing on higher quality. The algorithm runs until there are no two neighboring triangularcells together or the remaining cells do not create a convex quad. Figure 3.8 shows a second example ofcombining triangular cells into quad cells, merging cells C and D into a quad cell.A B C D E F G H I J K L M N A C D E F I J K L M N H B G A C D E F I J K L M N H B G (a) Before and after merging cellsA C B G B G A B C B A E I B A I A B A E C D F B C F C D C B C D E C D D E (b) Possible combinations for cells A and BFigure 3.7: 1st example of combining neighboring triangular cells into a quad cell383.4. 3D Mesh Optimization TechniquesA B C D E F G H I J K L M N A C D E F I J K L M N H B G A C D E F I J K L M N H B G (a) Before and after merging cellsA C B G B G A B C B A E I B A I A B A E C D F B C F C D C B C D E C D D E (b) Possible combinations for cells C and DFigure 3.8: 2nd example of combining neighboring triangular cells into a quad cell3.4 3D Mesh Optimization TechniquesIn this section the mesh operations described in section Sec. 2.2 such as swapping in metric space, insertionto improve cell alignment, smoothing defined by target elements and coarsening based on a metric lengthconflict graph are defined to adapt a mesh to conform a metric in 3D. The criteria which these operationsare based on is explained with further detail.3.4.1 Swapping in 3DSwapping is based on face swapping and edge swapping as mentioned in Sec. 2.2.1. Swapping a face fromtwo neighboring cells turns them into three tetrahedra, while swapping an edge from three neighboringcells turns them into two tetrahedra. Swapping the adjacent face of four neighboring tetrahedra or of twoboundary neighboring tetrahedra is also considered, as shown in Fig. 3.9.a b c d e f a b c d e f a b c d e a b c d e a b c d e a b c d e 2 3 4 4 2 2 Number of Cells Figure 3.9: 3D swapping configurations393.4. 3D Mesh Optimization TechniquesLet five points a, b, c, d, and e define two different sets of neighboring tetrahedra. A set A defined bytwo neighboring tetrahedra abcd and cabe represents the actual configuration of a set of cells in the mesh,and a set B defined by three neighboring tetrahedra abed, bced, and caed represent the cells after swappingas shown at the top of Fig. 3.9. Swapping face abc in set A will result in set B while swapping edge ed in setB will result in set A.The swapping decision is based on two different swapping criteria. The first swapping decision is basedon maximizing the minimum sine of the dihedral angle. This measure penalizes both large and small dihedralangles if computed in the physical space; Freitag and Ollivier-Gooch [18] find it to be the most effectivemeasure they considered.The dihedral angle is computed by cross multiplying the faces unit normals:nˆ =(b−a)× (c−a)‖(b−a)× (c−a)‖ (3.7)sinθFi,Fi+1 = ‖nˆi× nˆi+1‖ (3.8)A comparison between sinθFi,Fi+1 computed in the physical space and sinM θFi,Fi+1 computed in the metricspace is shown in Fig. 3.10 and Fig. 3.11. If the minimum dihedral angle is maximized in the configurationafter swapping in the metric space, the cells will swap:min(sinM θFi,Fi+1)SetB > min(sinM θFi,Fi+1)SetA (3.9)The same swapping criteria applies for evaluating swapping between 2-2 and 4-4 cell configurations asshown in Fig. 3.9.A B C D a b c d = Sin𝜃𝐴𝐶 = Sin𝜃𝐶𝐷 = Sin𝜃𝐴𝐵 = Sin𝜃𝐴𝐷 = Sin𝜃𝐵𝐶 = Sin𝜃𝐵𝐷 Physical space element Figure 3.10: 3D maxmin sine in physical space403.4. 3D Mesh Optimization TechniquesFigure 3.11: 3D maxmin sine in metric spaceThe second criterion is based on the volume-length quality measure. This measure was suggested byParthasarathy, Graichen, and Hathaway [61], which will be used in metric space. A comparison betweenvolume and length in physical space and in metric space is shown in Fig. 3.12 and Fig. 3.13.Let a tetrahedron be defined by abcd as shown in Fig. 3.13 where a, b, c, and d are vertices mappedinto the metric space. The quality measure in Eq. 3.10 is computed as the volume of a tetrahedron in metricspace from Eq. 3.11 space divided by the cube of its root-mean-squared edge metric length from Eq. 3.12.We set it so the highest quality is minimized.QV LM =∣∣∣∣1−(6√2)(VolMl3rmsM)∣∣∣∣ (3.10)VolM = bx−axby−aybz−az · cx−axcy−aycz−az× dx−axdy−aydz−az6(3.11)lrmsM =√√√√√(l2aM + l2bM + l2cM + l2dM + l2eM + l2fM)6 (3.12)This quality is computed for the all the tetrahedra in set A before swapping and in set B after swapping,if the quality measured is minimized after swapping, the cells will be swapped:min(QV LM)SetB < min(QV LM)SetA (3.13)413.4. 3D Mesh Optimization Techniques= 𝑙𝑎𝑐 = 𝑙𝑐𝑑 = 𝑙𝑎𝑏 = 𝑙𝑏𝑐 = 𝑙𝑏𝑑 Physical space element a b c d a b c d a c b d b c = 𝑙𝑎𝑑 a d Figure 3.12: 3D volume-length in physical spacea b c d Metric space element a c b d b c a d a b c d = 𝑙𝑀𝑎𝑐 = 𝑙𝑀𝑐𝑑 = 𝑙𝑀𝑎𝑏 = 𝑙𝑀𝑏𝑐 = 𝑙𝑀𝑏𝑑 = 𝑙𝑀𝑎𝑑 Figure 3.13: 3D volume-length in metric space423.4. 3D Mesh Optimization Techniques3.4.2 Insertion in 3DInsertion in 3D uses an analogous approach to insertion in 2D. examples of adding a vertex at the midpointof an interior edge and at the midpoint of a boundary edge are shown in Fig. 3.14 and Fig. 3.15.A B A B A B C A B C Figure 3.14: 3D insertion on an interior edgeA B A B C A B A B C A B A B C A B A B C Figure 3.15: 3D insertion on a boundary edgeWe insert points into the mesh with two different objectives at different stages of the algorithm, using asimilar approach in the interior of the mesh as at the boundary. A point is added near the midsection of anedge to improve a quality measure. In our primary insertion we seek to improve cell shape and orientation.Let an edge be defined by points AB and let C be the new point to add. Whether the edge is an interioredge or a boundary edge, adding a point will split n neighboring tetrahedra into 2n tetrahedra. The quality iscomputed for every neighboring cell of the edge and for the two new cells that are created from it if a pointwas added at the midpoint of the edge. If the quality of the neighboring cell before the insertion is worsethan the worst quality of the two new cells created, this counts as a favorable split (Eq. 3.14). Otherwise,this counts as a non-favorable split (Eq. 3.15). We insert a point at the midpoint of an edge if the number offavorable insertions is higher than the non favorable insertions Eq. 3.16. That is, we insert for interior andboundary edges if:Qϒ(Before)> max(Qϒ(NewCells)) Favorable Insertion (3.14)433.4. 3D Mesh Optimization TechniquesQϒ(Before)< max(Qϒ(NewCells)) NonFavorable Insertion (3.15)Total Favorable Insertion > Total NonFavorable Insertion (3.16)in which we use a ShapeOrient quality metric. We also insert to eliminate long edges in the metric space.For both boundary and interior edges, our criterion is to split edges where:lmaxM < lM(E) (3.17)To prioritize edges and faces for splitting we use insertion queues. For insertion based on quality thequeue is prioritized based on how much the quality improves by adding a point: higher improvement getsa higher priority in the queue. For insertion based on metric length the queue is prioritized based on size:edges with a larger metric length get a higher priority in the queue. To avoid creating cells and faces smallerthan desired, an edge will not be queued if any new edge created from the insertion of the vertex would beshorter than lminM .In contrast to the 2D insertion process, where the location of a new point is computed along the edgebased on a quality criterion in a further step, the location of a new point in 3D is set to the midpoint of theedge. In 3D we found that computing the location of a new point along the edge based on a quality criteriondid not alter the location significantly.3.4.3 Coarsening in 3DThe criterion for coarsening in 3D is based on eliminating edges that are too short in the metric space byusing a conflict graph. The conflict graph is designed to choose which vertices in a mesh should remain afterthe mesh is coarsened. The conflict graph is constructed by identifying, for each vertex Vi all the verticesthat are too close to it in the metric space; that is, vertices Vj for which lM(Vi,Vj) < lminM . We then select amaximal set of vertices that are not in conflict; this means that no two vertices marked to be kept will be inconflict, but if we were to keep any other additional vertex, we would create a conflict.The selection of final vertices to keep is hierarchical, working first on boundary curves, then on theentire boundary and in the interior of the domain. For each level of this hierarchical process, the vertices atthat level are marked as active, meaning only the specified type of vertices will be considered.An example of the conflict graph process on a small portion of a mesh is shown in Fig. 3.16. Let thevertices A, B, C, and D be marked to keep. Using the conflict graph, we can easily identify all vertices thatare too close to any of these and mark those conflicting vertices for deletion. Then we add another unmarkedvertex at this level to keep and mark its neighbors for deletion. This process continues until all vertices at agiven level are marked, then proceeds to the next level. Further details on the conflict graph are explained in[54].The vertices marked for deletion will be removed from the mesh by edge collapse, as described inSec. 2.2.3.443.4. 3D Mesh Optimization TechniquesA B C D (a) Vertices tagged to keepA B C D (b) Checking for conflicts in first layerA B C D (c) Checking for conflicts in second layerFigure 3.16: Process of searching for conflicts3.4.4 Smoothing in 3DSmoothing is analogous to the 2D approach, using Target-Matrix Optimization Paradigm (TMOP) as im-plemented in Mesquite [19] and setting a ShapeOrient quality metric. We are using global smoothing andsetting the target elements as right tetrahedra as explained in Sec. 2.4. Due to the limitations in Mesquite,only interior vertices are smoothed, thus surface vertices are not affected by smoothing.453.5. 3D General Algorithm3.5 3D General AlgorithmThe general algorithm in 3D to adapt a mesh based on a metric field is shown in Alg.3.2, which is analogousto the 2D algorithm. A metric, which will define how the mesh will be adapted, must be previously definedby an analytical function or an error estimate and set at every vertex of the mesh. A frontal view and anisometric view of the surface and of a crinkle cut of the interior of the mesh per step of the process is shownin Fig. 3.17 to Fig. 3.21. The mesh size in the adaptation process is shown in Table 3.2. The metric used toadapt the mesh in this example described in Sec.4.4.1. It represents a quarter circle anisotropic feature alongthe z axis, as the cells approach to the center of the anisotropic feature the aspect ratio of the cells increases.Once the metric is set for the initial mesh shown in Fig. 3.17, point insertion based on quality improve-ment is made. All edges in the mesh are evaluated to determine which ones will be split as described inSec. 3.4.2. Every new cell created in this process is re-evaluated for addition to the queue. This processcontinues until the insertion queue is empty. All new cells created are evaluated for swapping as describedin Sec. 3.4.1; the mesh after swapping is shown in Fig. 3.18. Swapping is done first with a volume-lengthcriterion in the metric space to capture a higher number of cells in this process; this improves the qualitycells of some cells but also creates bad quality cells. The bad cells created are treated by the second swap-ping criterion, maximizing the minimum sine of the dihedral angle in the metric space. This is intended tofix the bad cells created in the previous swapping. Both swapping criteria are only used for the queued cellscreated by the quality insertion process, every other swapping process in the algorithm only maximizes theminimum sine of the dihedral angle in the metric space. The previous process is only enhancing the shapeand orientation of the cells in the mesh, for which a relative low number of points is added. Smoothingis done with right angle tetrahedral targets based on a shape and orientation quality metric as described inSec. 2.22. Since the algorithm uses a global smoothing, all faces are then evaluated for swapping maximiz-ing the minimum sine of the dihedral angle in metric space. A second quality insertion step follows thesmoothing; the results are shown in Fig. 3.19.All the previous steps described are based on an eigen-decomposition of the metric which is used toadapt the mesh to improve the shape and orientation of the cells. For the next steps, the main goal is toimprove the size of the cells. An insertion based on metric length to reduce the size of edges larger than amaximum metric length is done (Fig. 3.20), as described in Sec. 3.4.2. Coarsening based on metric lengthincreases the size of edges smaller than a minimum metric length is done after, as described in Sec. 3.4.3.All cells are evaluated for swapping after coarsening (Fig. 3.21), as described in Sec. 3.4.1. A final insertionbased on quality is done to reduce any big cells produced from coarsening.Process # Cells # VerticesStep 1.- Set Metric (Initial Mesh) 5053 1332Step 2.- Quality Insertion 14008 3584Step 5.- Quality Insertion 19300 4588Step 6.- Long Edge Insertion 95071 18933Step 8.- Swapping after coarsening 83468 16605Table 3.2: Mesh size during 3D adaptation process463.5. 3D General AlgorithmAlgorithm 3.2 3D general algorithmStep 1 A metric previously defined is set in all the vertices Step 2 Insertion based on quality, every edge of a new cell is reevaluated for insertion and swapping with both criteria Step 3 Smoothing with right angle targets based on shape and orientation Step 4 All the cells are evaluated for swapping recursively by maximizing the minimum angle Step 5 2nd Insertion based on quality, reevaluating new edges for insertion and swapping with both criteria Step 6 Edges larger than a max previously defined are split Step 7 Edges smaller than a min are deleted, decision for which vertex to delete based on quality Step 8 All the cells are evaluated for swapping recursively by maximizing the minimum angle Quality Insertion Smooth Long Edge Insertion Coarsen Quality Insertion Swap Queue both criteria Swap Swap Set metric Repeat step 2 Compute metric Swap Queue max min angle 473.5. 3D General Algorithm(a) Frontal view (b) Frontal crinkle cut view(c) Isometric view (d) Isometric crinkle cut viewFigure 3.17: 3D adaptation process: initial mesh483.5. 3D General Algorithm(a) Frontal view (b) Frontal crinkle cut view(c) Isometric view (d) Isometric crinkle cut viewFigure 3.18: 3D adaptation process: 1st quality insertion (2nd step)493.5. 3D General Algorithm(a) Frontal view (b) Frontal crinkle cut view(c) Isometric view (d) Isometric crinkle cut viewFigure 3.19: 3D adaptation process: 2nd quality insertion (5th step)503.5. 3D General Algorithm(a) Frontal view (b) Frontal crinkle cut view(c) Isometric view (d) Isometric crinkle cut viewFigure 3.20: 3D adaptation process: long edge insertion (6th step)513.5. 3D General Algorithm(a) Frontal view (b) Frontal crinkle cut view(c) Isometric view (d) Isometric crinkle cut viewFigure 3.21: 3D adaptation process: swapping after coarsening (8th step)52Chapter 4ResultsIn this chapter the proposed algorithm is used to adapt a coarse mesh to conform to a metric. A coarsesquare domain mesh is adapted to a metric defined analytically by two different given functions for 2D inSec. 4.1. In Sec. 4.2 the adaptation effects of this algorithm are illustrated with a subsonic viscous flow witha boundary layer on a NACA 0012 airfoil coarse mesh, in which the metric is computed using a 2nd orderreconstruction error metric based on the solution. The meshes adapted by this algorithm are merged intoquad cells to show how this can be facilitated by the alignment obtained in Sec. 4.3. A coarse cube domainmesh is adapted to a metric defined analytically by given functions for 3D in Sec. 4.4, the functions used todefine the metric are analogous to the 2D functions used in Sec. 4.1.4.1 2D Analytic Metric ExamplesIn this section the metric is defined analytically on an initial mesh by a given function which was definedin [42]. Each time a point is inserted, the metric is interpolated as shown in Eq. 2.12. After the terminationof the algorithm, the metric is reset with the given function for the next iteration. The metric length rangeallowed will be bounded by lminM =√22 and lmaxM =√2. In all the examples in this section a square mesh witha domain of [1,0]× [0,1], 462 cells and 264 vertices as shown in Fig. 4.1 was used as an initial mesh.Figure 4.1: Initial mesh534.1. 2D Analytic Metric Examples4.1.1 2D Quarter Circle FunctionThe metric function in this example represents a quarter circle curved anisotropic feature:Mcircle2D(x,y) =(h−21 cos2 θ +h−22 sin2 θ(h−21 −h−22)cosθ sinθ(h−21 −h−22)cosθ sinθ h−21 sin2 θ +h−22 cos2 θ)(4.1)where h1 = min(0.002× 5α ,hmax), h2 = min(0.05× 2α ,hmax), hmax = 0.1, θ = arctan(x,y), α = 10×∣∣∣0.75−√x2+ y2∣∣∣.The results of the first two iterations of the proposed algorithm are shown in Fig. 4.2 and Fig. 4.3 and amore detailed description of the first iteration for this example is shown in Fig. 3.6 at Sec. 3.2. The resultingmesh size after every adaptation process is shown in Table 4.1. Analyzing the results obtained from thefirst and second adaptation, our proposed algorithm has captured the desired curved anisotropic feature withgood alignment after the first adaptation. A better defined curved anisotropic feature is obtained after thesecond adaptation without losing the alignment. Figure 4.2: 1st adaptation for a quarter circle curved anisotropic feature in 2D with a close up Figure 4.3: 2nd adaptation for a quarter circle curved anisotropic feature in 2D with a close up544.1. 2D Analytic Metric ExamplesA Shape Orientation quality criterion as described in Sec. 2.4.3 is applied after the adaptations and ahistogram of the quality in a log scale are shown in Fig.4.4 and Fig.4.5. Since the criterion is minimized, anideal quality would be zero. Results from both iterations show a good approximation of the desired quartercircle with good alignment as shown in the close ups of the highlighted regions and a high percent of cellswith good shape orientation quality according to the metric.Figure 4.4: 1st adaptation, Shape Orientation quality and histogram for a quarter circle curved anisotropicfeature in 2DFigure 4.5: 2nd adaptation, Shape Orientation quality and histogram for a quarter circle curved anisotropicfeature in 2DMesh Cells VerticesInitial mesh 462 264Mesh after 1st Adapt 3537 1853Mesh after 2nd Adapt 2957 1557Table 4.1: Mesh size in adaptation for a quarter circle curved anisotropic feature in 2D554.1. 2D Analytic Metric Examples4.1.2 2D Cross FunctionThe metric function in this example represents two anisotropic features; two lines cross perpendicular toeach other creating a cross shape at the center of the mesh:Mcross2D(x,y) =(h−2x 00 h−2y)(4.2)where hx =min(2αx×hmin,hmax), hy =min(2αy×hmin,hmax), hmin = 0.005, hmax = 0.1, αx = 20×|x−0.5|,αy = 20×|y−0.5|.The results from the second and fourth iterations of the proposed algorithm are shown in Fig. 4.6 andFig.4.7 and quality results are shown in Fig.4.8 for the first adaptation and Fig.4.9 for the second adaptation.The resulting mesh size after every adaptation process is shown in Table 4.2. Figure 4.6: 1st adaptation for perpendicular anisotropic features in 2D with a close up Figure 4.7: 2nd adaptation for perpendicular anisotropic features in 2D with a close up564.1. 2D Analytic Metric ExamplesResults from both iterations show a good approximation of the desired anisotropic features, having highaspect ratio cells in the arms of the cross and smaller cells at the center. Better defined anisotropic featuresare obtained as the mesh keeps adapting with the proposed algorithm. A high percent of cells have shapeorientation quality approaching to zero, thus having a good quality according to this metric.Figure 4.8: 1st adaptation, Shape Orientation quality and histogram for perpendicular anisotropic featuresin 2DFigure 4.9: 2nd adaptation, Shape Orientation quality and histogram for perpendicular anisotropic featuresin 2DMesh Cells VerticesInitial mesh 462 264Mesh after 2nd Adapt 3061 1594Mesh after 4th Adapt 3127 1626Table 4.2: Mesh size in adaptation for perpendicular anisotropic features in 2D574.2. 2D Flow Around an Airfoil Examples4.2 2D Flow Around an Airfoil ExamplesThe adaptation effects of this algorithm are illustrated with a subsonic viscous flow with a boundary layer.The flow is at Ma = 0.5 , with an angle of attack of α = 0◦ and a Reynolds number of Re = 5000. This caseshows physical behavior that commonly occurs in aerodynamics; an anisotropic feature parallel to the flowwithin the boundary layer and wake. The solution is computed using a second order vertex-centered finitevolume solver as described in [53]. The solver is based on least squares reconstruction [55], Roe’s scheme[65], and Newton-GMRES [48, 49] for convergence. The viscous terms are discretized as described in [51].(a) Farfield view(b) Closeup viewFigure 4.10: Initial mesh around a NACA 0012 airfoil with 2543 vertices and 4941 cellsThe initial mesh is constructed using the isotropic mesh generator in GRUMMP [52] based on Delaunaytriangulation. The farfield boundary is a circle of radius 500 chords centered at the leading edge of the584.2. 2D Flow Around an Airfoil Examplesthe airfoil, thus the errors due to boundary placement can be neglected. A farfield view and a closeup onthe initial mesh is shown in Fig. 4.10. This mesh is clearly ill-suited to compute viscous flow around theairfoil due to anisotropic features parallel to the flow. The mesh is adapted by using the proposed algorithmpreviously described. The metric is computed using the 2nd order reconstruction error metric proposedin [58], which is based on the Hessian of the solution. The metric length range used for this example islminM = ldes− τ and lmaxM = ldes + τ (computed as proposed in [58]), where ldes is the edge length from thedesired equilateral triangle with an average metric area of the cells in the mesh, and τ is a value defined bya lower bound on minimum angle in the metric space.4.2.1 Viscous Subsonic FlowThe metric is calculated based on the initial solution and used for the adaptation algorithm. After everyadaptation, the solution is re-computed, and the metric is computed based on the new solution at everyvertex of the mesh.To prove this algorithm improves the adaptation efficiency and to show the effectiveness of the adapta-tion scheme in producing aligned meshes without relying on smoothing, the solution for the initial mesh aswell as the first three adaptations shown in Fig. 4.11 were obtained without anisotropic smoothing. Afterevery adaptation the accuracy of the solution is improved and mesh alignment near the boundary layer andthe wake is obtained. A close up view of the highlighted region of the resulting mesh after the third adap-tation Fig. 4.11d is shown in Fig. 4.12, the airfoil top/bottom surface and downstream of the trailing edge.Our scheme produces anisotropic cells with good alignment even with real-world metrics that are not aswell-behaved as the analytic tests in the previous section. Also, note that the mesh anisotropy is highest atthe edge of the boundary layer, at the center of the wake, and in the recirculation region. At all of these lo-cations the curvature of the solution — notably the x-component of the velocity — is at a maximum locally.Figure 4.4 and Fig. 4.5 show a Shape Orientation quality criterion after the third adaptation and a histogramof the quality in a log scale, in which good quality has been obtained in portions of the anisotropic feature.A comparison between a smoothed mesh and a non-smoothed mesh after the 3rd adaptation is shown inFig. 4.14. Both meshes reached similar results after the 3rd adaptation having small variation between themas shown in a closeup from both meshes in Fig. 4.14b and Fig. 4.14c. These small variations do not alter thesolution significantly, thus the solution from the non-smoothed mesh will be analyzed for this section.For a zero angle of attack, the pressure distribution and the wall shear stress ideally should behavesymmetrically, and the lift coefficient should be zero; this value is a measure of how well we capture thesymmetry. Due to the non-symmetric mesh, numerical results are expected to have a small variation. Acomparison of drag and lift coefficients with other schemes is shown in Table 4.3 showing similar resultswith fewer vertices for our algorithm. Fig. 4.15 shows the convergence for the drag coefficient due topressure and viscosity which remains without any major variation through out several iterations.The pressure distribution and the wall shear stress per adaptation for the top surface are shown inFig. 4.16. After every iteration, a more consistent distribution is obtained. The comparison between thetop surface and the bottom surface is shown in Fig. 4.17; as a result of mesh asymmetry, these are quitesimilar but not identical.594.2. 2D Flow Around an Airfoil Examples0.20 0.40 0.050 0.60 Mach(a) Solution of initial mesh (2543 vertices)0.20 0.40 0.050 0.60 Mach(b) Solution of non-smoothed mesh after 1st adaptation (4235 vertices)0.20 0.40 0.050 0.60 Mach(c) Solution of non-smoothed mesh after 2nd adaptation(6688 vertices)0.20 0.40 0.050 0.60 Mach(d) Solution of non-smoothed mesh after 3rd adaptation(10335 vertices)Figure 4.11: Mesh and solutions for a subsonic viscous flow over a NACA 0012, Ma= 0.5, Re= 5000,α =0◦604.2. 2D Flow Around an Airfoil Examples(a) Top surface(b) Bottom surface(c) Wake near the trailing edgeFigure 4.12: Closeup of airfoil NACA 0012 after 3rd adaptation614.2. 2D Flow Around an Airfoil ExamplesSmoothed meshNon-smoothed mesh(a) Smoothed and non-smoothed mesh(b) Non-smoothed closeup(c) Smoothed closeupFigure 4.14: Comparison between smoothed and non-smoothed mesh after 3rd adaptation624.2. 2D Flow Around an Airfoil Examples0.400 0.800 1.20 1.60 2.00 2.40 2.80 0 3ShapeOrient(a) ShapeOrient quality(b) Quality histogramFigure 4.13: ShapeOrient quality after 3rd adaptation634.2. 2D Flow Around an Airfoil ExamplesScheme Mesh Size CL CD,P CD,νXSeparationxchordInitial mesh 2543 0.01449 0.03179 0.02939 —-Non-smoothed 1st Adapt 4235 0.01644 0.006784 0.02806 —-Non-smoothed 2nd Adapt 6688 0.001424 0.02237 0.03258 0.7971Non-smoothed 3rd Adapt 10335 0.004885 0.02212 0.03231 0.8092Sharbatdar / Pagnutti [68]2543 0.015074 0.6108 0.2925 —-3899 0.009244 0.5571 0.2067 0.82369140 0.000203 0.022429 0.03412 0.817619964 0.000108 0.02138 0.03374 0.8141Pagnutti [57]3885 0.005786 0.4858 0.1723 0.85087645 0.021304 0.1297 0.03797 0.848617386 0.006265 0.05612 0.03643 0.8449ARC2D 320 X 128 —- 0.0221 0.0321 0.824Mavriplis [45] 320 X 64 —- 0.02290 0.0332 0.814Radespiel [63] 512 X 128 —- 0.02242 0.0330 0.8144th order Spectral Volume [27] —- —- 0.02206 0.0324 0.813Table 4.3: Comparison of lift and drag for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦(a) Drag coefficient due to pressure(b) Drag coefficient due to viscosityFigure 4.15: Drag convergence for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦644.2. 2D Flow Around an Airfoil Examples(a) Wall shear stress(b) Pressure coefficientFigure 4.16: Distribution behavior per adaptation for NACA 0012, Ma = 0.5, Re = 5000,α = 0◦(a) Wall shear stress(b) Pressure coefficientFigure 4.17: Distribution behavior for top and bottom surfaces of 3rd Adaptation for NACA 0012, Ma =0.5, Re = 5000,α = 0◦654.3. Tri-to-Quad Conversion Results4.3 Tri-to-Quad Conversion ResultsAn advantage of adapting a triangular mesh for a quasi structure is the ability to merge triangular cellsinto quadrilateral cells. To show how the results from the proposed algorithm can facilitate this process,the adapted meshes modified from Sec. 4.1 and Sec. 4.2 will be merged into a quad mesh. This is doneby a simple procedure of matching the best pair of cells into a quad cell without any mesh modificationoperations, as described in Sec. 3.3. The adapted mesh with a curve anisotropic feature after the seconditeration from Sec. 4.1.1 is merged through this process. The results are shown in Fig. 4.18a and a close upof the highlighted region of the curved anisotropic feature is shown in Fig. 4.18b. The mesh contains 3199triangular cells before the conversion and 1788 triangular/quadrilateral cells after the conversion from which78.9% are quadrilateral cells and 21.1% are triangular cells. (a) Mesh after the 2nd adapt(b) Close up of curved anisotropic featureFigure 4.18: Tri-to-quad for curved anisotropic feature664.3. Tri-to-Quad Conversion ResultsThe adapted mesh with two perpendicular anisotropic features after the fourth iteration from Sec. 4.1.2is converted into a quadrilateral/triangular mesh. The results are shown in Fig. 4.19 and a close up of thehighlighted region of a straight anisotropic feature is shown in Fig.4.19b. The mesh contains 3127 triangularcells before the conversion and 1743 triangular/quadrilateral cells after the conversion from which 79.4%are quadrilateral cells and 20.6% are triangular cells. (a) Mesh after 4th adapt(b) Close up of straight anisotropic featureFigure 4.19: Tri-to-quad for perpendicular anisotropic featuresThe adapted mesh from the subsonic viscous flow example after the third iteration from Sec. 4.2.1 ismerged through the same process. The results are shown in Fig. 4.20 and a close up of the top surface of theairfoil is shown in Fig. 4.20b and a close up of the bottom surface of the airfoil is shown in Fig. 4.20c. Themesh contains 20274 triangular cells before the conversion and 11537 triangular/quadrilateral cells after theconversion from which 75.73% are quadrilateral cells and 24.27% are triangular cells.674.3. Tri-to-Quad Conversion Results(a) Mesh after 3rd adapt(b) Close up of top surface of airfoil(c) Close up of bottom surface of airfoilFigure 4.20: Tri-to-quad for subsonic viscous flow over a NACA 0012, Ma = 0.5, Re = 5000,α = 0◦Using a simple operation which chooses the best quadrilateral between two neighboring triangular cellswithout any mesh modification operations as described in Sec. 3.3 was enough to reach up to a 75% conver-sion of the mesh into quads with a reliable alignment in all the cases above.684.4. 3D Analytic Metric Examples4.4 3D Analytic Metric ExamplesIn this section the metric is defined analytically by a given function on an initial mesh, these functions areanalogous to the ones used in 2D in Sec. 4.1. Each time a point is inserted, the metric is interpolated asshown in Eq. 2.12. After the termination of the algorithm, the metric is reset for the entire mesh with thegiven function for the next iteration. The metric length range allowed will be bounded by lminM =√22 andlmaxM =√2. In all the examples in this section a cube mesh with a domain of [1,0,0]× [0,1,0]× [0,0,1],5053 cells and 1332 vertices as shown in Fig. 4.21 was used as an initial mesh.(a) Surface view(b) Crinkle cut viewFigure 4.21: Initial mesh694.4. 3D Analytic Metric Examples4.4.1 3D Quarter Circle FunctionThe metric function in this example is an extension to 3D of Eq.4.1, which represents a quarter circle curvedanisotropic feature along the z axis:Mcircle3D(x,y,z) = h−21 cos2 θ +h−22 sin2 θ(h−21 −h−22)cosθ sinθ 0(h−21 −h−22)cosθ sinθ h−21 sin2 θ +h−22 cos2 θ 00 0 100 (4.3)where h1 = min(0.002× 5α ,hmax), h2 = min(0.05× 2α ,hmax), hmax = 0.1, θ = arctan(x,y), α = 10×∣∣∣0.75−√x2+ y2∣∣∣.The results of the first two iterations of the proposed algorithm in 3D are shown in Fig.4.22 and Fig.4.23.A close up of the highlighted region of the frontal view shows how the anisotropic curve feature is capturedas the mesh is adapted. This anisotropic feature is captured with good alignment at the surface of the mesh.A better defined anisotropic curved feature is obtained after the second adaptation without affecting thealignment of the cells. Figure 4.24 and Fig. 4.25 show a frontal view and an isometric view of both thesurface of the mesh and the interior with a crinkle cut. Even though alignment is obtained at the surface ofthe mesh, and in the interior of the mesh the anisotropic curve feature is captured, alignment of the mesh inthe interior is poor compared with the surface.Quality is measured after every adaptation using a Shape Orientation quality criterion as described inSec. 2.4.3 and a histogram on a log scale is shown in Fig. 4.26 and Fig. 4.27. Since the quality is minimized,an ideal quality would be zero. Quality results from both iterations show a good quality of the cells atthe anisotropic curved feature at the surface, but poor quality in the interior according to the metric dueto the lack of alignment. After the second iteration, the quality histogram does not show any significantimprovement, meaning another iteration of the proposed algorithm would not increase the general quality ofthe mesh. To get the appropriate alignment a more rigorous criterion is used for swapping which providesgood results for the surface of the mesh, but due to the higher complexity in 3D at the interior of the mesh,it also creates bad quality cells which affect the quality and alignment of the cells. The resulting mesh sizeafter adaptation cycle is shown in Table 4.4.Mesh Cells VerticesInitial mesh 5053 1332Mesh after 1st Adapt 89964 17811Mesh after 2nd Adapt 146967 27993Table 4.4: Mesh size in adaptation for a quarter circle curved anisotropic feature in 3D704.4. 3D Analytic Metric Examples Figure 4.22: Close up of 1st adaptation quarter circle curved anisotropic feature in 3D Figure 4.23: Close up of 2nd adaptation quarter circle curved anisotropic feature in 3D714.4. 3D Analytic Metric Examples(a) Frontal view(b) Isometric view(c) Crinkle cut of isometric viewFigure 4.24: 1st adaptation for a quarter circle curved anisotropic feature in 3D724.4. 3D Analytic Metric Examples(a) Frontal view(b) Isometric view(c) Crinkle cut of isometric viewFigure 4.25: 2nd adaptation for a quarter circle curved anisotropic feature in 3D734.4. 3D Analytic Metric Examples0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(a) Frontal view0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(b) Isometric view(c) Quality histogramFigure 4.26: ShapeOrient quality for 1st adaptation of quarter circle curved anisotropic feature in 3D744.4. 3D Analytic Metric Examples0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(a) Frontal view0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(b) Isometric view(c) Quality histogramFigure 4.27: ShapeOrient quality for 2nd adaptation of quarter circle curved anisotropic feature in 3D754.4. 3D Analytic Metric Examples4.4.2 3D Cross FunctionThe metric function in this example is an extension to 3D of Eq. 4.2, which represents two anisotropicfeatures; two lines cross perpendicular to each other creating a cross shape at the center of the mesh alongthe z axis:Mcross3D(x,y) = h−2x 0 00 h−2y 00 0 100 (4.4)where hx =min(2αx×hmin,hmax), hy =min(2αy×hmin,hmax), hmin = 0.005, hmax = 0.1, αx = 20×|x−0.5|,αy = 20×|y−0.5|.The results from the second and fourth iterations of the proposed algorithm in 3D are shown in Fig.4.28and Fig. 4.30; This images show a frontal view and an isometric view of both the surface of the mesh andthe interior with a crinkle cut. A close up of the highlighted region of the isometric view shows how theanisotropic features are captured as the mesh is adapted in Fig. 4.29 and Fig. 4.31. The images show a closeup view of the front surface, the side, and the top of the mesh. The anisotropic features are captured withalignment at the surface and definition of the anisotropic features improves with more adaptation cycles.Similar to the previous example, alignment is obtained at the surface. In the interior of the mesh, theanisotropic features are captured but with poor alignment compared with the surface.Quality is measured after every adaptation using a Shape Orientation quality criterion and a histogramon a log scale is shown in Fig.4.32 and Fig.4.33. Poor quality of the interior cells according to the metric forboth iterations is obtained due to the lack of alignment in the interior cells. As in the previous example, asthe mesh keeps adapting the quality histogram does not show any significant improvement, meaning moreiterations of the proposed algorithm would not increase the general quality of the mesh. A similar behaviorto the previous example is obtained, getting good results for the surface of the mesh and bad quality cellswhich affect the quality and alignment of the interior cells. The resulting mesh size after every adaptationcycle is shown in Table 4.5.Mesh Cells VerticesInitial mesh 5053 1332Mesh after 2nd Adapt 67961 13359Mesh after 4th Adapt 97978 18850Table 4.5: Mesh size in adaptation for perpendicular anisotropic features in 3D764.4. 3D Analytic Metric Examples(a) Frontal view(b) Isometric viewFigure 4.28: 2nd adaptation for perpendicular anisotropic features in 3D774.4. 3D Analytic Metric Examples (a) Highlighted regions of close up(b) Frontal anisotropic feature close up(c) Lateral anisotropic feature close up(d) Top anisotropic feature close upFigure 4.29: Close up of 2nd adaptation perpendicular anisotropic features in 3D784.4. 3D Analytic Metric Examples(a) Frontal view(b) Isometric viewFigure 4.30: 4th Adaptation for perpendicular anisotropic features in 3D794.4. 3D Analytic Metric Examples (a) Highlighted regions of close up(b) Frontal anisotropic feature close up(c) Lateral anisotropic feature close up(d) Top anisotropic feature close upFigure 4.31: Close up of 4th adaptation perpendicular anisotropic features in 3D804.4. 3D Analytic Metric Examples0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(a) Frontal view0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(b) Isometric view(c) Quality histogramFigure 4.32: ShapeOrient quality for 2nd adaptation of perpendicular anisotropic features in 3D814.4. 3D Analytic Metric Examples0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(a) Frontal view0.250 0.500 0.750 0 1ShapeOrient0.250 0.500 0.750 0 1ShapeOrient(b) Isometric view(c) Quality histogramFigure 4.33: ShapeOrient quality for 4th adaptation of perpendicular anisotropic features in 3D82Chapter 5Concluding Remarks5.1 Summary and ConclusionsAn existing anisotropic adaptation scheme by Pagnutti [57], which produces meshes without any regularstructure, was extended by Sharbatdar by adding a anisotropic node-shifting based with Target Matrix Op-timization Paradigm in order to achieve a quasi-structured mesh [68]. The previous adaptation schemedepended on smoothing to obtain an aligned-quasi structured mesh in the anisotropic parts of the domain.Even though an improvement was obtained, a higher alignment was still desired.This thesis has presented a new algorithm for anisotropic mesh adaptation with good alignment proper-ties. The effectiveness of the new algorithm in producing aligned meshes does not only rely on smoothing toimprove alignment of the cells, thus a higher cell alignment is obtained. The new algorithm is based on fourmesh optimization techniques: point insertion designed to improve mesh alignment as well as conformingto a metric, swapping in the metric space, point movement defined by target elements from a metric, andpoint deletion based on quality and metric length. This algorithm improves mesh alignment by improvingdifferent aspects of the mesh in separate steps. First, the shape and orientation of the cells in the meshis improved by inserting points in strategic places based on a shape orientation criterion, which matchesanisotropic elements in the mesh with the anisotropic features of the flow properties. Once the shape andorientation of the cells is improved, the size of the cells is adapted according to a maximum and minimummetric length.The proposed algorithm is used to modify coarse meshes in 2D according to metrics defined analyticallyby a given function and according to a second order reconstruction error metric for a subsonic viscous flowwith a boundary layer. The 3D extension to the algorithm is used to modify a coarse mesh accordingto metrics defined analytically, analogous to the 2D test cases. The adapted meshes include quantitativemeasures of anisotropic mesh quality based on a shape orientation criterion. An advantage of creating aquasi-structured triangular mesh is the ability to merge triangular cells into well-shaped quadrilateral cells.The adapted meshes in 2D are merged with a simple merging criterion into quad cells to show an advantageof the cell alignment obtained.The ability of the proposed algorithm to create aligned anisotropic meshes has been shown by adaptingcoarse isotropic meshes to conform a metric, which is computed by an analytical function or with an errormetric based on the Hessian of a solution. Two different examples where shown with analytic metrics for2D in a square domain. In each of these examples, the algorithm captured the anisotropic features withgood quality and reliable alignment. With additional adaptation cycles, better refined anisotropic featureswere obtained without affecting the cell alignment. To show that our scheme produces anisotropic cells with835.2. Recommendationsgood alignment even with real-world metrics that are not as well-behaved as the analytic metrics, a subsonicviscous flow around an NACA-0012 airfoil was examined. The subsonic viscous flow leads to anisotropicbehavior at the surface of the airfoil due to the presence of boundary layers, which was captured by theadapted mesh created with our algorithm. A comparison with the previous schemes shows that our newalgorithm can increase the solution accuracy with fewer points even without anisotropic smoothing. Theadapted meshes were merged into quad cells by a simple process of matching cells for which a high per-centage of the mesh merged into quad cells, this to show an advantage of the effectiveness of the algorithmin obtaining well aligned meshes.The algorithm is extended to 3D for which two different cases were shown. The meshes were adaptedby an analytic metric function analogous to the 2D metric functions. For both cases, the anisotropic featuresare captured by our algorithm with good alignment at the surface of the mesh. In the interior of the mesh,the anisotropy is captured but with a lack of alignment; thus interior cells have a lower quality for a shapeorientation quality metric.5.2 RecommendationsThe proposed algorithm was developed to adaptively produce anisotropic metric-orthogonal meshes con-forming a metric, which can be applied to any metric-based CFD algorithm with a minimum adjustment;however, this approach was tested for analytical metric functions cases and for a two-dimensional finitevolume subsonic viscous flow around a NACA-0012 airfoil. Testing the algorithm with different parameterssuch as a transonic inviscid flow in 2D or any three-dimensional flow might provide a better idea of theweaknesses of the algorithm and how to improve them. Another potential improvement could be obtainedby using adjoint-based error estimation methods instead of solution approximation methods for computingthe metric, combining the algorithm with an adjoint-based adaptation could provide greater accuracy.Even though the extended 3D algorithm captured the anisotropic features with good alignment at thesurface for the test cases, there is a lack of alignment at the interior of the mesh which results in bad qualitycells. This is caused by the higher complexity of 3D unstructured meshes. Once a bad cell is created inthe interior, it affects how the neighboring cells are adapted and so forth. Adding specific operations toeliminate or avoid creating these bad cells while obtaining good alignment is a potential improvement to thealgorithm.The mesh operations used in the proposed algorithm are also a potential improvement. Smoothing isdone by setting right angle targets. Using a variety of targets based on the metric could improve the resultsobtained, especially in 3D where a wide variety of different cells can be obtained. Smoothing boundariesin 2D and surfaces in 3D is also another potential improvement. Insertion is currently done only along anexisting edge; a wider range of specific cases for bad cells might improve the capacity for removing thesecells from the mesh.The goal of this work is to improve solution accuracy by generating a higher quality mesh, even thoughthe obtained solutions are more accurate, there is potential improvement in the optimization of the process.One key point in this algorithm is based on computing quality with TMOP for which a wide variety of845.2. Recommendationspossibilities is analyzed to obtain the quality per cell. Defining how to measure quality while minimizingthe different possibilities per cell without affecting the results is a potential improvement for optimization,specially in 3D.As further research is done in improving high-order CFD solvers and anisotropic mesh generators, theadvantages of obtaining a well aligned mesh while maintaining good accuracy at a low computational costwill increase, for which a new approach with potential improvement was presented in this thesis.85Bibliography[1] Frederic Alauzet. Metric-based anisotropic mesh adaptation. Course material, CEA-EDF-INRIASchools, June 2010.[2] Frederic Alauzet and Pascal Frey. Estimate d’erreur geometrique et metriques anisotropes pourl’adaptation de maillage: Partie 1: Aspects theoriques. RR-4789,INRIA, March 2003.[3] Frederic Alauzet and Pascal Frey. Estimate d’erreur geometrique et metriques anisotropes pourl’adaptation de maillage: Partie 2: Exemples d’ applications. RR-4789,INRIA, March 2003.[4] Frederic Alauzet, Adrien Loseille, Alain Dervieux, and Pascal Frey. Multi-dimensional continuousmetric for mesh adaptation. 15th International Meshing Roundtable, pages 191–214, September 2006.[5] R. Balasubramanian and J.C. Newman III. Comparison of adjoint-based and feature-based grid adap-tation for functional outputs. International Journal for Numerical Methods in Fluids, 53:1541–1569,2007.[6] David Bommes, Bruno Levy, Nico Pietroni, Enrico Puppo, Claudio Silva, Marco Tarini, and DenisZorin. State of the art in quad meshing. Computer Graphics Forum, pages 1–24, 2012.[7] Frank J. Bossen and Paul S. Heckbert. A pliant method for anisotropic mesh generation. In Proceedingsof the Fifth International Meshing Roundtable, pages 63–74. Sandia National Laboratories, October1996.[8] Carlo L. Bottasso. Anisotropic mesh adaption by metric-driven optimization. International Journalfor Numerical Methods in Engineering, 60:597 – 639, 2004.[9] Philip Caplan, Robert Haimes, David L. Darmofal, and Steven R. Allmaras. Anisotropic advancingfront quadrangulation with optimization-based smoothing. 22nd International Meshing Roundtable,October 2013.[10] G. Carpentieri, B. Koren, and M.J.L. van Tooren. Adjoint-based aerodynamic shape optimization onunstructured meshes. Journal of Computational Physics, 224:267–287, 2007.[11] A. Choudhary and Christopher J. Roy. A truncation error-based approach to understanding and im-proving mesh quality in cfd. 50th AIAA Aerospace Sciences Meeting including the New HorizonsForum and Aerospace Exposition, 2012.86Bibliography[12] Francois Courty, David Leservoisier, Paul-Louis George, and Alain Dervieux. Continuous metrics andmesh adaptation. Applied Numerical Mathematics, 56:117–145, 2006.[13] David Lee Davidson. The role of computational fluid dynamics in process industries. The Bridge,Linking Engineering and Society, 32(4):9–14, 2002.[14] Eduardo F. D’Azevedo and R. Bruce Simpson. On optimal triangular meshes for minimizing thegradient error. Numerische Mathematik, 59:321–348, 1991.[15] R.P. Dwight. Heuristic a posteriori estimation of error due to dissipation in finite volume schemes andapplication to mesh adaptation. Journal of Computational Physics, 227:2845–2863, 2008.[16] Mehdi Falsafioon, Sina Arabi, Francois Guibault, and Ricardo Camarero. Smoothing of unstructuredgrids. 27th IAHR Symposium on Hydraulic Machinery and Systems, 2014.[17] Lori Freitag, Mark Jones, and Paul Plassmann. An efficient parallel algorithm for mesh smoothing. InProceedings of the Fourth International Meshing Roundtable, pages 47–58. Sandia National Labora-tories, October 1995.[18] Lori A. Freitag and Carl F. Ollivier-Gooch. Tetrahedral mesh improvement using swapping andsmoothing. International Journal for Numerical Methods in Engineering, 40(21):3979–4002, Novem-ber 1997.[19] Lori Freitag-Diachin, Patrick Knupp, and Jason Kraftcheck. Mesquite: Mesh Quality ImprovementToolkit User’s Guide. Technical report, Center for Computation, Computers, Information, and Math-ematics, 2009.[20] Paul L. George, F. Hecht, and Marie G. Vallet. Creation of internal points in voronoi’s type method.control and adaptation. Adv. Eng. Software, 13(5-6):303-312, 1991.[21] Paul-Louis George and Houman Borouchaki. Delaunay Triangulation and Meshing: Application toFinite Elements. Hermes, Paris, 1998.[22] Serge Gosselin. Delaunay refinement mesh generation of curve-bounded domains. PhD thesis, TheUniversity of British Columbia, Department of Mechanical Engineering, 2009.[23] Topraj Gurung, Daniel Laney, Peter Lindstrom, and Jarek Rossignac. Squad: Compact representationfor triangle meshes. EUROGRAPHICS, 30(2), 2011.[24] M. J. Harris. Flow featured aligned mesh generation and adaptation. Master’s thesis, University ofSheffield, Department of Mechanical Engineering, April 2013.[25] Weizhang Huang. Metric tensors for anisotropic mesh generation. Journal of Computational Physics,204:633–665, 2005.[26] Weizhang Huang. Mathematical principles of anisotropic mesh adaptation. Communications in Com-putational Physics, 1(2):276–310, 2006.87Bibliography[27] Ravi Kannan and Zhijian Wang. Improving the high order spectral volume formulation using a diffu-sion regulator. Communications in Computational Physics, 12(No. 1):247–260, 2012.[28] Patrick Knupp. Algebraic mesh quality metrics. SIAM Journal of Scientific Computing, 23(1):133–149, 2001.[29] Patrick Knupp. Local 2D metrics for mesh optimization in the target-matrix paradigm. TechnicalReport SAND2006-7382J, Sandia National Laboratories, 2006.[30] Patrick Knupp. Analysis of 2D rotational-invariant non-barrier metrics in the target-matrix paradigm.Technical Report SAND2008-8219P, Sandia National Laboratories, 2008.[31] Patrick Knupp. Label-invariant mesh quality metrics. In Proceedings of the Eighteenth InternationalMeshing Roundtable, pages 139–155, 2009.[32] Patrick Knupp. Target-matrix construction algorithms. Technical report, Sandia National Laboratories,Albuquerque, 2009. Applied Mathematics & Applications Department.[33] Patrick Knupp. Introducing the target-matrix paradigm for mesh optimization via node-movement.Proc. of 19th International Meshing Roundtable, Springer-Verlag, pages 67–84, 2010.[34] Patrick Knupp and Emile van der Zee. Convexity of mesh optimization metrics using a target-matrixparadigm. Technical Report SAND2006-4975J, Sandia National Laboratories, 2006.[35] François Labelle and Jonathan Richard Shewchuk. Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation. In Proceedings of the Nineteenth Annual Symposium on Compu-tational Geometry, pages 191–200. Association for Computing Machinery, 2003.[36] Burkhard Lehner, Bernd Hamann, and Georg Umlauf. Generalized swap operation for tetrahedriza-tions. In Scientific Visualization: Advanced Concepts’10, pages 30–44, 2010.[37] Xiangrong Li, Mark S. Shephard, and Mark W. Beall. 3-D anisotropic mesh adaptation by meshmodifications. Computer Methods in Applied Mechanics and Engineering, 194(48–49):4915–4950,2005.[38] S. H. Lo. 3D anisotropic mesh refinement in compliance with a general metric specification. FiniteElement Analysis and Design, 38(1):3–19, 2001.[39] Harvard Lomax, Thomas H. Pulliam, and David W. Zingg. Fundamentals of computational fluiddynamics. Springer-Verlag, Berlin, 2003.[40] Adrien Loseille. Metric-orthogonal anisotropic mesh generation. In Proc of 23rd International Mesh-ing Roundtable, volume 82, pages 403 – 415. Procedia Engineering, 2014.[41] Adrien Loseille and Rainald Lohner. Anisotropic adaptive simulations in aerodynamics. In 48th AIAAAerospace Sciences Meeting, 2010,.88Bibliography[42] David Marcum and Frederic Alauzet. Aligned metric-based anisotropic solution adaptive mesh gen-eration. In Proc. of 23rd International Meshing Roundtable, volume 82, pages 428–444. ProcediaEngineering, 2014.[43] Mark M. Mathis and Darren J. Kerbyson. A general performance model of stuctured and unstructuredmesh particle transport computations. The Journal of SuperComputing, (34):181–199, 2005.[44] D. Mavriplis. A discrete adjoint-based approach for optimization problems on three-dimensional un-structured meshes. 2006. AIAA 2006-0050.[45] Dimitri Mavriplis and Antony Jameson. Multigrid solution of the Euler equations on unstructured andadaptive meshes. ICASE Report No. 87-53, NASA Langley Research Center, 1987.[46] Dimitri J. Mavriplis. Unstructured mesh related issues in computational fluid dynamics (cfd) basedanalysis and design. In Springer-Verlag, editor, 11th International Meshing Roundtable, page 143,September 2002.[47] Fotis Mavriplis. CFD in aerospace in the new millenium. Canadian Aeronautics and Space Journal,46(4):167–179, 2000.[48] Christopher Michalak and Carl Ollivier-Gooch. Matrix-explicit GMRES for a higher-order accurateinviscid compressible flow solver. In Proceedings of the Eighteenth AIAA Computational Fluid Dy-namics Conference, 2007.[49] Amir Nejat and Carl Ollivier-Gooch. A high-order accurate unstructured finite volume Newton-Krylovalgorithm for inviscid compressible flows. Journal of Computational Physics, 227(4):2592–2609,2008.[50] Marian Nemec, Michael J. Aftosmis, and Mathias Wintzer. Adjoint-based adaptive mesh refinementfor complex geometries. In Forty-sixth AIAA Aerospace Sciences Meeting, 2008. AIAA 2008-725.[51] Hiroaki Nishikawa. Two ways to extend diffusion schemes to navier-stokes schemes: Gradient formulaor upwind flux. In In Proceedings of the 20th AIAA Computational Fluid Dynamics Conference, 2011.[52] Carl F. Ollivier-Gooch. An Unstructured Mesh Improvement Toolkit with Application to Mesh Im-provement, Generation and (De-)Refinement. AIAA 98-0218, January 1998.[53] Carl F. Ollivier-Gooch. ANSLib: A scientific computing toolkit supporting rapid numerical solutionof practically any PDE. In Proceedings of the Eighth Annual Conference of the Computational FluidDynamics Society of Canada, pages 21–28. Société canadienne de CFD / CFD Society of Canada, June2000.[54] Carl F. Ollivier-Gooch. Coarsening unstructured meshes by edge contraction. International Journalfor Numerical Methods in Engineering, 57(3):391–414, May 2003.89Bibliography[55] Carl F. Ollivier-Gooch and Michael Van Altena. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 181(2):729–752, 2002.[56] Doug Pagnutti. Anisotropic adaptation: Metrics and meshes. Master’s thesis, The University of BritishColumbia, Department of Mechanical Engineering, 2008.[57] Doug Pagnutti and Carl Ollivier-Gooch. Delaunay-based anisotropic mesh refinement. In Proceedingsof the Seventeenth International Meshing Roundtable, 2008.[58] Doug Pagnutti and Carl Ollivier-Gooch. A generalized framework for high order anisotropic meshadaptation. Computers and Structures, 87:670–679, 2009.[59] Michael Park. Adjoint-based, three-dimensional error prediction and grid adaptation. AIAA 32ndFluid Dynamics Conference, 2002.[60] Michael A. Park and David L. Darmofal. Parallel anisotropic tetrahedral adaptation. Aerospace Sci-ences Meeting and Exhibit, 2008. AIAA 2008-917.[61] V. N. Parthasarathy, C. M. Graichen, and A. F. Hathaway. A comparison of tetrahedron quality mea-sures. Finite Elements in Analysis and Design, (9):309–320, 1991.[62] Shahyar Pirzadeh. Three-dimensional unstructured viscous grids by the advancing-layers method.American Institute of Aeronautics and Astronautics Journal, 34(1):43–49, January 1996.[63] R. Radespiel. A cell-vertex multigrid method for the Navier-Stokes equations. NASA TM-101557,1989.[64] J.F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, A. Johnen, and C. Geuzaine. Blossom-quad: anon-uniform qualdrilateral mesh generator using a minimum cost perfect matching algorithm. Inter-national Journal for Numerical Methods in Fluids, (89):1102–1119, 2012.[65] Philip L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal ofComputational Physics, 43:357–372, 1981.[66] Christopher Roy. Strategies for driving mesh adaptation in CFD. In Proceedings of the Forty-SeventhAIAA Aerospace Sciences Meeting. American Institute of Aeronautics and Astronautics, 2009.[67] Mahkame Sharbatdar. Anisotropic mesh adaptation: Recovering quasi-structured meshes. Master’sthesis, The University of British Columbia, Department of Mechanical Engineering, 2012.[68] Mahkame Sharbatdar and Carl Ollivier-Gooch. Anisotropic mesh adaptation: Recovering quasi-structured meshes. In Proceedings of the Fifty-First AIAA Aerospace Sciences Meeting, 2013.[69] Jonathan R. Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Computa-tional Geometry: Theory and Applications, 22(1–3):21–74, May 2002.90Bibliography[70] Jonathan Richard Shewchuk, Tamal Krishna Dey, and Siu-Wing Cheng. Delaunay Mesh Generation.Chapman & Hall / CRC Press, 2013.[71] Marco Tarini, Nico Pietroni, Paolo Cignoni, Daniele Panozzo, and Enrico Puppo. Practical quad meshsimplification. EUROGRAPHICS, 29(2), 2010.[72] David A. Venditti and David L. Darmofal. Grid adaptation for functional outputs: Application totwo-dimensional inviscid flows. Journal of Computational Physics, 175(1):40–69, February 2002.[73] David A. Venditti and David L. Darmofal. Grid adaptation for functional outputs: Application totwo-dimensional viscous flows. Journal of Computational Physics, 187:22–46, 2003.[74] Nigel P. Weatherill, Joe F. Thompson, and Bharat K. Soni, editors. Handbook of Grid Generation.CRC Press, 1999.[75] Masayuki Yano and David Darmofal. An optimization framework for anisotropic simplex mesh adap-tation: Application to aerodynamic flows. In Proceedings of the Fiftieth AIAA Aerospace SciencesMeeting. American Institute of Aeronautics and Astronautics, 2012.91
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Adaptive creation of orthogonal anisotropic triangular...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Adaptive creation of orthogonal anisotropic triangular meshes using target matrices Zuniga Vazquez, Jose David 2015
pdf
Page Metadata
Item Metadata
Title | Adaptive creation of orthogonal anisotropic triangular meshes using target matrices |
Creator |
Zuniga Vazquez, Jose David |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | We present a new procedure to adaptively produce anisotropic metric-orthogonal meshes in this thesis. The approach is based on mesh optimization techniques: point insertion designed to improve mesh alignment as well as conforming to a metric, swapping in the metric space, point movement defined by target elements from a metric, and point deletion based on quality and metric length. These techniques are intended to produce quasi-structured meshes which have the advantages of flexibility for complex geometries of unstructured meshes and the directional accuracy of structured meshes. The result is reliable alignment for anisotropic meshes reducing our previous work's reliance on smoothing for good alignment. The methodology is implemented in 2D and extended to 3D. Examples of analytical metrics and error estimation metrics on a numerical simulation of a flow are shown. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-12-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0220768 |
URI | http://hdl.handle.net/2429/55660 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2016-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2016_february_zunigavazquez_josedavid.pdf [ 31.61MB ]
- Metadata
- JSON: 24-1.0220768.json
- JSON-LD: 24-1.0220768-ld.json
- RDF/XML (Pretty): 24-1.0220768-rdf.xml
- RDF/JSON: 24-1.0220768-rdf.json
- Turtle: 24-1.0220768-turtle.txt
- N-Triples: 24-1.0220768-rdf-ntriples.txt
- Original Record: 24-1.0220768-source.json
- Full Text
- 24-1.0220768-fulltext.txt
- Citation
- 24-1.0220768.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0220768/manifest