Issues in Unstructured Mesh Generation

1. Issues in Unstructured Mesh Generation by Gopal Shevare ppts are based on Lecture notes by T J Barth (1994)

2. CONTENTS • Un-structured Grid • Unstructured Grids and Application of Graph Theory • Delaunay Triangulation and Voronoi diagrams • Methods of Delaunay Triangulations • Method of Non- Triangulation

3. Grid Generation Present Status Grid Generation is the gateway for CAE. Grids connects various simulations for designing a products Make grid generation easy and user friendly. If possible make it redundant Grid generation now encompasses each and everything from creating / importing geometry to export a single file that any solver may need irrespective what simulation code is simulating (Boundary Fitted Coordinate) BFC system is not grid generation! Grid generation has become interdisciplinary area where in CAGD, Computational Geometry, CFD, Computer Graphics experts contribute to make a package User may require execution of around 60 to 100 tasks before grid gets ready Present day simulations use millions of grid points Grids are a commodity. There is an ISO format for grids

4. Classification of Grids Adaptive Grids more points where solution changes rapidly Surface grid generation grids on curved surface Elliptic Grid Generation systems grids produced by solving elliptic PDEs Orthogonal systems grids lines are orthogonal to each other Moving grids grid points move with time Composite overset grid multiple grids overlapping Cartesian Adaptive Grids grid lines along x, y and z axes Hybrid grids more than one the of grids coexisting

5. Structured Grid Transform a given two-dimensional domain in to a rectangle (and three-dimensional domain in to a box) by a suitable affine transformation.  = f (x, y) and  = g (x, y). Mould arbitrary shaped physical domain (x, y, z) in to a rectangular shaped computational domain (, , ζ) ? Find three functions f, g and h such that  = f (x, y, z),  = g (x, y, z), and ζ = h(x,y,z). Lines corresponding to constant values of  = i  = i,  = j  = j, and ζ = k ζ = k for i = 1 to i max, j = 1 to j max and k = 1, k max in (x, y, z) plane are grid lines in (x,y, z) plane and intersection of these points gives (a) grid points (i, j) and also (b) quadrilateral cells.

6. Why Triangulation ? Why Triangulation & Tetrahedrons Theory of triagulations is highly developed Higher dimensions pose no problems Triangle is the polytope in 2D Simple to deal with in 2D And also on surfaces Easy to use in numerical calculations polygons Triangles only Mixture of triangles And quadrilaterals quadrilaterals only

7. Unstructured Grids View #2 Fill a given domain with simple shapes such as triangles, quadrilaterals, etc. so that the given domain is fully covered. Emphasis is on cells (triangles in the present case); there are grid points but no continuous lines or what can be called as grid lines. Structure imposed on grid points 6 2 3 7 5,8 1 4 10 9

8. Unstructured Meshes and Graphs • Graphs and Meshes • Graph theory directly affects the algorithms in unstructured grids • However they do not depend on the method of generating triangles • Planar graph • graph without self loop or parallel edges • Most important graph theory result is Euler formula which relates no of vertices, edges and faces of polyhedron (closed surface) • N(f) = n(e) – n(v) + 2 Euler formula • The polyhedron is flattened to 2D space – one face is at infinity 4 2 4 2 1 1 3 3 © Zeus Numerix 25th Jan -2006

9. Now there are two types of edges Internal – edges : shared by two faces Boundary edges : shared by only one face Edges for the face at infinity (boundary edges) (you may consider that face at infinity is absent) Edges binding regions, (holes) embedded in the domain (not a part of the domain) Euler formula becomes n(f) + n(v) = n(e) + 1 - n(h) Eqn. (1) But no of edges and no of faces are related to each other Unstructured Meshes and Graphs max d(f) =  i n(f) i =3 2 n(e)interior + n(e)boundary Eqn. (2) No of pointers from edges to triangles No of pointers from triangles to edges

10. Unstructured Meshes and Graphs • Subtracting 2 time Eqn (1) from Eqn (3) • n(f) = 2 n(v) - n(e)boundary -2 + 2n(h) • No of faces is approx. 2 times no of vertices • Substituting for n(f) in Eqn. (1) • n(e) = 3 n(v) - n(e)boundary - n(e) - 3 + 3n(h) • No of edges is approx. 3 time no of vertices! • What can happen if • n(f)  2 n(v) - n(e)boundary -2 + 2n(h) or • n(e) 3 n(v) - n(e)boundary - n(e) - 3 + 3n(h)

11. Basic Definitions Dirichlet Tessellation (Also called Voronoi Regions) The Dirchlet tessellation of a set of points is the pattern of convex regions, each being closer to some Point P in the point set than any other point in the set The edges of Voronoi polygons comprise the diagram. The idea of Voronoi regions naturally extends to 3D point edge

12. Basic Definitions - Duality Duality of Graph Let G be a planer graph of vertices, edges and cells then GDual the Dual graph of G exhibits the following properties: vertex in GDual corresponds to Face in G Face in G Face in GDual corresponds to vertex in G Vertex in G Edge in GDual corresponds to edge in G Edge in G • If an edge separates two faces • fi and fj in G, then the edge, GDual connects two vertices associated with faces fi and fj • GDual can be formed by (a)Medium segments, (b)Centroid segments, (c)Dirichlet tessellation

13. Triangulation triangle Voronoi diagram • Definition : The triangulation of a point set is defined as the dual of the Voronoi diagram of the set. • In 2D triangulation is formed by connecting if and only if their Voronoi regions a common border • If no four or more points are cocircular, then vertices of the Voronoi are circumcentres of triangles.

14. Let use define N0 = n(v), N1 = n(e), N2 = n(f), N3 = n() Vertices edges faces volume Generalised Euler formula is d (-1)k N k= 0 (Euler formula in Rd); by definition N-1 = 1 k = -1 For the surface of polyhedron in 3D, -1 + N0 – N1 + N2 – 1 = 0 or N(f) + N(v) = N(e) +2 For polyhedral volumes n(f) + n(v) = n(e) + n() Like 2D case this formula does not take care of boundary effects because it is valid for four-dimensional polytope Unstructured Meshes and Graphs

15. Geometric data : (x, y) coordinates of vertices Connectivity data : Vertex based v-v, v-e, v-f Edge based e-v, e-e, e-f Face based f-v, f-e, f-f All 9 relations are not stored because Large RAM is required Consistency checks must be carried out on relations Only some relations are stored. All the remaining relations are derived from the stored relations The stored connectivity relations are called data structures: Some well-known data structures: Finite element data structure. Edge structure. Out-degree structure. Quad-edge structure. Winged-edge data structure. Data Structures

16. Data Structures Storing of Structured grids in unique Edges in 3D are always formed by (i+1, j, k) and (i, j, k) etc…. Faces are always formed by (i, j, k), (i+1, j, k), (I, j+1, k) and (i+1, j+1, k)… How to Store Unstructured Meshes? Collection of ensure does not ensure edges are connected Collection of faces does not ensure faces form a 3D cell and that faces not intersecting . . . . .

17. Data Structures Edge structure Ei Fi, Fj For planar meshes Edge is list of pair of vertices Ideal for FVM Standard FE data structure Ti Vl, Vm, Vn Sequence of vertices matters f-v connectivity Naturally extend to 3D Extensively used in EFM of solids and fluids Quad-edge structure Ei-1/2 Ei+1/2, Ej, Ek Edges are stored as a pair of directed edges Stores vertex at the origin Stores previous and next edge Out-degree structure Vl  Vm, Vn, Vo But this is not v-v connectivity Max three out degree required

18. Graph Operations on Unstructured Grids • Planar Graph with Minimum Out-Degree • Theorem: Every planar graph has a 3 bounded orientation; i.e. each edge of a planar • graph can be assigned a direction such that the maximum number of edges pointing • outward from any vertex is less than three • Proof. • Find a vertex such that incident boundary edges that connected to at most • maximum two. Other boundary vertices and any number of interior edges. • Such a vertex is called reducible vertex. It can be proved that at least two vertices • always exit • The two edges connected to boundary vertices are directed outward • All the remaining edges are directed inward • The directed edges are removed from the graph • The process is repeated on the remaining graph

19. Unstructured Meshes and Graphs Planar Graph with Minimum Out-Degree : Graphically explained

20. Unstructured Meshes and Graphs Graph Ordering Techniques Ordering can have a large influence on CPU time and RAM requirements. Vertices close to each other should be stored close to each other Most of the grid generation algorithms produce poor ordering If a Boolean matrix of size [N x N] matrix stores connectivity of vertices such that (i,j)th entry is 1 if ith and jth vertices are connected, many 1s should be lying close to diagonal of the matrix As obtained form the algorithm Ordered to minimize bandwidth

21. A typical unstructured mesh Entries of Non-zero Laplacian Laplacian of reordered mesh Cuthill-McKee algorithm

22. Unstructured Meshes and Graphs Graph Ordering Techniques Rosen Algorithm Step I : Determine the band width and choose a pair of vertices (say i and j) such that i < j Step II : If increasing i or decreasing j reduces band width then (increase i or decrease j and goto Step I) Step III : If increasing i and decreasing j does not change the band width then (increase I or decrease j and goto Step I) Comments Produces good ordering but is very expensive especially for large graphs

23. Unstructured Meshes and Graphs Graph Ordering Techniques Cuthill-McKee Algorithm Step I : Find the vertex with lowest degree Step II : Find all the neighbouring vertices incident on root. Step III : Order the neghbouring vertices by increasing degree Step IV : Choose the neighbour with lowest degree if not previously not ordered. Step V : Continue with Step II till all vertices are ordered. Comments Very efficient Recommended in LD decomposition of matrices by reverse renumbering

24. Unstructured Meshes and Graphs Graph Partitioning Graph partitioning is required for load balancing in parallel computing on unstructured grids Requirements Each ‘processor’ should get as per its capability Communication between ‘processors’ should be minimum. Communication graph This is graph in which vertices are ‘processors’ and edge indicates that they need to share some information between. The colours in the graph show ‘communication cycles’ 1 4 Communication cycles 1st cycle: 1-4 & 2-3 2nd cycle: 1-2, 3-4 3rd cycle: 1-3 4 1 3 2 2 3 Communication graph Mesh partitioned

25. Unstructured Meshes and Graphs Graph Partitioning Vizing’s theorem Any graph of maximum degree  (i.e. number of edges incident on vertex) can be coloured with using n colours such that  <= n <= +1 Cost of Communication is A + Bm, where A – time required to initiate the message B – rate of data transfer and m – length of the message n Total cost =  (A + Bm) 1 To reduce cost (i) reduce  (reduce number of messages, n), but 1 <= n < = P-1 or (ii) reduce the length of message Partitioning for minimum message length Partitioning for minimum 

26. Unstructured Meshes and Graphs Methods of Graph Partitioning – Coordinate bisection Step I : Sort the face centroids either horizontally or vertically Step II : Separate the triangles such that no of faces are the same Step III : Repeat Step II recursively Comment Scheme is efficient in partitioning the triangulation. Partition is sub-optimal as the message length is not minimum

27. Unstructured Meshes and Graphs Methods of Graph Partitioning – CutHill-McKee Algorithm Step I : Produce a dual of given mesh Step II : Reorder the dual graph Step III : Set the partition at the median or closet to the median Comment Scheme is efficient in partitioning the triangulation. Reduces number of sub-graphs Tends to produce sub-optimal as the message length is not minimum

28. Methods of Coordinate Bi-section

29. Methods of CutHill-McKee Algorithm the message length is not minimum

30. Methods of Spectral Partitioning the message length is not minimum

38. 1. Uniqueness: The triangulation is unique, assuming that no four sites are circular. Properties of Triangulation A D Point free interior Incircle test for the triangle Fails if ∠DCB + ∠DAB≤π 2. The circumcircle criteria A triangulation for N > 2 site is if and only if the circumcircle of every interior triangle is point-free B C D Interior is not empty A D A Incircle test for the triangle is true if ∠DCB + ∠DAB ≥π B C Since interior angles of quadrilaterals sum to 360 degrees, swapping diagonal can be used to produce triangle B C

39. 3. Edge circle property A triangulation of sites is if and only if there exists some circle passing through the end points of each and every edge which is point free Properties of Triangulation 4. Edge circle property for constrained triangulation A triangulation of sites is constrained triangulation if for each and every edge there exists some circle passing through its end points containing other site in the triangulation which is visible to the edge Point not visible to edge Constrained edge

40. 4. Equiangularity property triangulation maximises the minimum angle of the triangulation. (Hence Delauanay triangulation  Maxmin triangulation. The property is valid globally and locally. The property can be used to produce Deluanay triangulation by swapping of diagonal of convex quadrilateral. Properties of Triangulation Minimum angle Only one edge is acceptable Non- triangulation made Edge swapping can not be used in concave quadrilaterals Minimum angle

41. Properties of Triangulation 5. Minimum Containment property triangulation minimises the maximum containment circle over entire triangulation Containment circle is the smallest triangle enclosing the three vertices of a triangle Containment circle is (i) the circumcircle if the triangle is acute (ii) circle with diameter equal to the largest side if triangle is obtuse. Containment circle  circumcircle for acute triangles Containment circle  circle with diameter equal to the longest edge

42. Properties of Triangulation 6. Nearest neighbour Property In triangulation an edge is obtained by joining the nearest neighbour But the nearest neighbour edge is one of the many edges Nearest neighbour edge Non-nearest neighbour edge 7. Minimal roughness Property Suppose a property fi is given at all vertices Vi i = 1 to N Any triangulation will produce piecewise linear interpolation surface for this data triangulation minimises the integral given by  [(f/x)2 + (f/y)2 ] dx dy for any f

43. Divide and conquer algorithm Incremental insertion algorithm Bowyer algorithm Watson algorithm Green and Sibson algorithm Global edge swapping algorithm. Generation of Triangulation Advancing front algorithm Graded triangulation algorithm

44. Generation of Triangulation Divide and conquer algorithm • Overall algorithm • Sites are assumed to be pre-specified • Partition the sites in to Left (L) and right (R) half planes • Triangulate each half recursively • Merge each half recursively R triangulation L triangulation • Points to note while merging • Only L-R or R-L edges are created • No new L-L or R-R edges do not get generated as L-triangulation or R-triangulation is Delaunany. But the merging may lead to deleting L-L and R-R edges • Merge algorithm • Connect vertices with minimum y in the left and right triangulation as the first guess to the lowest common tangent (LCT) of the L and R triangulations • contd….. First guess for lowest common tangent (LCT)

45. C B A Generation of Triangulation • Merge algorithm • Allow a ‘bubble’ passing through LCT (edge AB) to rise till hits a point (C) from either L or R triangulation. From the new edge CB or CA • Use CB as the LCT and use rising bubble to hit a new point on either L or R triangulation • Delete the L-L and/or R-R edges to complete the merge either after every creation of LCT or after all LTC are obtained. rising bubble new edge

46. Incremental Insertion Algorithms - Common Features Points to be inserted • New points are added to existing triangulation one by one • New points are assumed to lie in the existing triangulation • All incremental insertion algorithms start with a phantom triangle large enough to enclose all the given points • The edges formed by these points, if given as a part of input data is initially ignored. These edges will be built later after all the points are inserted. • Every insertion of new point requires locating (i) a triangle or (ii) circumcircles containing this point. • In mesh adaptation the location is known. In initial triangulation some search is required Phantom triangle Edge ignored

47. Incremental Insertion Algorithms • Search is is carried out by using either quad-tree or split-tree data structures as they give 0(NlogN) complexity • Alternately, ‘walking’ in the direction towards new vertex using graphs can also be used as these algorithms are simpler. A worst case in this gives 0(N2) complexity • Changing the connectivity of vertices to accommodate criterion. This may require at worst 0(N2) changes to the data structure. • Remove the phantom triangle and the recreate boundary or edges given as input • Normally the triangles thus produced require some smoothening Triangle containing the new point new point Circumcircles containing the new point

48. Incremental Triangulation - Bowyer’s Algorithm Assume some triangulation exists; Insert a new point (site) into Voronoi diagram Find out all the vertices in the Voronoi diagram closer to new point than the existing points (sites). Delete these vertices (shown as v1, v2, v3 and v4) Delete the edges connecting the vertices if both of them have been deleted Produce the new vertices and the edges by calculating the vertices Join the sites to form the triangulation by joining the sites separated by an edge Existing points Existing points new points Voronoi diagram v1 v4 v2 v3

49. Incremental Triangulation - Watson’s Algorithm Assume some triangulation exists; Insert a new point into existing triangles Find the triangles in which the new point lies. Staring with this triangles perform a tree search to find all triangles with circumcircle containing the new point This produces a polygonal cavity around the new point. Connect the new point with vertices of the polygonal cavity Swap edges if triangles are not Delaunay

50. Global edge swapping by Lawson Choose an interior edge Select the quadrilateral with the edge as a diagonal Swap to satisfy Delaunay criterion Continue swapping edges till there is no swapping