Multithread parallelization of lepp-bisection algorithms

Elsevier Editorial System(tm) for Applied Numerical Mathematics Manuscript Draft Manuscript Number: APNUM-D-10-00228R2 Title: Multithread parallelization of lepp-bisection algorithms Article Type: Special Issue - WONAPDE 2010 Keywords: Longest edge bisection, Triangulation refinement, Lepp-bisection algorithm, Parallel multithread refinement, Parallel lepp-bisection algorithm, Finite element method Corresponding Author: Maria-Cecilia Rivara, Ph.D. Corresponding Author's Institution: Universidad de Chile First Author: Maria-Cecilia Rivara, Ph.D. Order of Authors: Maria-Cecilia Rivara, Ph.D.; Pedro Rodriguez; Rafael Montenegro, Ph.D.; Gaston Jorquera Manuscript Region of Origin: CHILE Abstract: Abstract Longest edge (nested) algorithms for triangulation refinement in two dimensions are able to produce hierarchies of quality and nested irregular triangulations as needed both for adaptive finite element methods and for multigrid methods. They can be formulated in terms of the longest edge propagation path (Lepp) and terminal edge concepts, to refine the target triangles and some related neighbors. We discuss a parallel multithread algorithm, where every thread is in charge of refining a triangle t and its associated Lepp neighbors. The thread manages a changing Lepp(t) (ordered set of increasing triangles) both to find a last longest (terminal) edge and to refine the pair of triangles sharing this edge. The process is repeated until triangle t is destroyed. We discuss the algorithm, related synchronization issues, and the properties inherited from the serial algorithm. We present an empirical study that shows that a reasonably efficient parallel method with good scalability was obtained. Multithread parallelization of lepp-bisection algorithms Maria-Cecilia Rivaraa, Pedro Rodrigueza,b, Rafael Montenegroc, Gaston Jorqueraa aDepartment of Computer Science, University of Chile, Santiago, Chile bDepartment of Information Systems, University of Bio-Bio, Concepci´on, Chile cUniversity of Las Palmas de Gran Canaria, University Institute for Intelligent Systems and Numerical Applications for Engineering, Spain Abstract Longest edge (nested) algorithms for triangulation refinement in two dimensions are able to produce hierarchies of quality and nested irregular triangulations as needed both for adaptive finite element methods and for multigrid methods. They can be formulated in terms of the longest edge propagation path (Lepp) and terminal edge concepts, to refine the target triangles and some related neighbors. We discuss a parallel multithread algorithm, where every thread is in charge of refining a triangle t and its associated Lepp neighbors. The thread manages a changing Lepp(t) (ordered set of increasing triangles) both to find a last longest (terminal) edge and to refine the pair of triangles sharing this edge. The process is repeated until triangle t is destroyed. We discuss the algorithm, related synchronization issues, and the properties inherited from the serial algorithm. We present an empirical study that shows that a reasonably efficient parallel method with good scalability was obtained. Keywords: Longest Edge Bisection, Triangulation Refinement, Parallel Multithread Refinement, Lepp-bisection Algorithm, Finite Element Method 1. Introduction Triangular mesh generation for finite element methods has been extensively studied and addressed by engineers and numerical analysts since the seventies. Finite element methods are widely used numerical techniques for the practical analysis of complex physical problems modeled by partial differential equations, which require appropriate discretizations of the associated geometries. Because of their flexibility, triangulations are preferred practical tools. Pioneer works on triangulations for finite element methods and related issues are due to Babuska Email addresses: mcrivar[email protected]e.cl (Maria-Cecilia Rivara), prodri[email protected]le.cl (Pedro Rodriguez), [email protected]lpg.es (Rafael Montenegro), gjorqu[email protected]ile.cl (Gaston Jorquera) Preprint submitted to Elsevier July 14, 2011 *Manuscript Click here to view linked References and Aziz [4], Lawson [27], Sibson [51]. Since then, intensive research on practical mesh generation has been performed. See e.g. the papers of Baker [6, 7], Bouraki and George [9], Jones and Plassmann [23, 24], Williams [55]. Computational methods for generating and refining triangular and tetrahedral finite element meshes in 2 and 3-dimensions can be roughly classified as Delaunay based methods, [6, 19, 7, 17, 50, 47] and methods based on the partition of triangles and tetrahedra [39, 33, 34, 28]. Longest edge refinement algorithms for triangulations, based on the longest edge bisection of triangles (obtained by joining the midpoint of the longest edge with the opposite vertex) were especially designed to deal with adaptive multigrid finite element methods [33, 34, 35, 36]. They are able to perform iterative local refinement by essentially maintaining the geometric quality of the input mesh as needed in finite element applications; they produce hierarchies of quality, smooth and nested irregular triangulations as required for non-structured multigrid methods. The properties of these algorithms are inherited from the non-degeneracy properties of the iterative longest edge bisection of triangles, [45, 52, 53, 1, 21] and are summarized in section 2.2. The longest edge algorithms were generalized for 3-dimensional mesh refinement [38, 30], as well as for the derefinement or coarsening of the mesh [37]. Improved longest edge algorithms based on using the concepts of terminal edges and longest edge propagating paths were also developed [40, 39, 43]. The algorithms have been used for developing software for partial differential equations. See e.g. the applications of Nambiar et al [31], Muthukrishnan et al. [30]. Based on the longest edge idea over Delaunay meshes, Lepp-Delaunay algorithms for triangulation improvement have been also developed [39, 41, 44] 2. Longest edge refinement algorithms in 2-dimensions. Roughly speaking the problem is the following: given a conforming, nondegenerate triangulation, construct a locally refined triangulation, with a desired resolution and such that the smallest (or the largest angle) is bounded. To simplify we introduce a subregion Rto define the refinement area; and a condition over the longest-edge of the triangles to fix the desired resolution. We can consider the following subproblems: Area Refinement. Given a quality acceptable triangulation (a triangulation with angles greater than or equal to an angle α) of a polygonal region D, construct a locally refined triangulation such that the longest edge of the triangles that intersect the refinement region Rare less than δ. Point / Edge Refinement. Here the refinement is performed around one vertex or along a boundary side, until all the triangles that intersect the vertex or the boundary edge are less than δ. Adaptive finite element refinement. In the adaptive finite element context, the refinement region is repeatedly defined as a subset of triangles Sref of the current triangulation (not necessarily connected) where the error of the finite element solution is too big to be acceptable [36]. 2 Given the input mesh, the algorithm locally and iteratively refines the triangles of a changing Sref set (or those intersecting the refinement region R) and some neighboring triangles. The new points introduced in the mesh are midpoints of the longest edge of some triangles of either of the input mesh or of some refined nested meshes. The longest edge bisection guarantees in a natural way the construction of non-degenerate and smooth irregular triangulations whose geometrical properties only depend on the initial mesh. In practice, for adaptive triangulation refinement, at each step longest-edge algorithms produce a refined conforming and guaranteed-quality output triangulation by performing selective longest edge bisection of the triangles of Sref and some (longest edge) related neighbors. 2.1. A serial Lepp-bisection algorithm An edge Eis called a terminal edge in triangulation τif Eis the longest edge of every triangle that shares E, while the triangles that share Eare called terminal triangles [39, 40]. Note that in 2-dimensions either Eis shared by two terminal triangles t1,t2if Eis an interior edge, or Eis contained in a single terminal triangle t1if Eis a boundary (constrained) edge. See Figure 1(a) where edge AB is an interior terminal edge shared by two terminal triangles t2, t3. For any triangle t0in τ, the longest edge propagating path of t0, called Lepp(t0), is the ordered sequence {tj}N+1 0, where tjis the neighbor triangle on a longest edge of tj−1, and longest-edge (tj)>longest-edge (tj−1), for j=1,... N. Edge E= longest-edge(tN+1) = longest-edge(tN) is an interior terminal edge in τand this condition determines N. Consequently either Eis shared by a couple of terminal triangles (tN, tN+1) if Eis an interior edge in τ, or Eis shared by a unique terminal triangle tNwith boundary (constrained) longest edge. See Figure 1(a) for an illustration of these ideas, where Lepp(t0)=(t0, t1, t2, t3). The Lepp-bisection algorithm can be simply described as follows: each triangle tin Sref is refined by finding Lepp(t), a pair of terminal triangles t1, t2and associated terminal edge l. Then the longest edge bisection of t1, t2is performed by the midpoint of l. The process is repeated until tis destroyed (refined) in the mesh. An efficient formulation of the algorithm where Lepp(t0) is not repeatedly recomputed, but repeatedly updated starting from the non-modified part of the previous Lepp(t0), is presented below. To this end we use a dynamic ordered list that stores pointers to the increasing triangles of (partial and full) Lepp(t0), while t0remains in the changing mesh. This will be the basis to develop a parallel algorithm in section 4. Lepp-Bisection Algorithm Input : a quality triangulation, τ, and a set Sref of triangles to be refined for each t in Sref do Insert-Lepp-Points(τ, t) end for Insert-Lepp-Points(τ, t0) Initialize Ordered-List (associated dynamically to Lepp(t0)) with t0 3 while Ordered-List is not empty do Find last triangle tNin Ordered-List Find longest edge neighbor tN+1 of tNand add it to Ordered-List (tN+1 can be null if longest edge of tNis over the boundary) if tN,tN+1 share a terminal edge or tN+1 is null then Perform longest-edge bisection of tN,tN+1 by midpoint of common terminal edge Eliminate tN,tN+1 from Ordered-List end if end while Note that the refinement task is performed when each Lepp(t0) is fully computed and a terminal edge is identified, by using a very local refinement operation (the bisection of pairs of terminal triangles that share a common longest edge). This guarantees that the mesh is conforming throughout the whole refinement process. This improves previous longest edge algorithms [34] that produced intermediate non-conforming meshes. Also the algorithm is free of non-robustness issues, since this do not depend of complex computations, and the selected points are midpoints of existing previous edges. Figure 1 illustrates the refinement of triangle t0in the input triangulation (a). Triangulation (b) shows the first point inserted, while triangulation (c) corresponds to the final triangulation obtained where the vertices are numbered in the creation order. In order to achieve this work, the full Lepp computation is performed three times to respectively insert points P1, P2, P3. The first Lepp computation includes triangles t0, t1, t2, t3, being t2, t3terminal triangles. Once the refinement of these triangles is performed, Lepp(t0) is partially recomputed starting from t1(Figure 1(b)). This now includes triangles t0, t1,˜ t2being t1,˜ t2 terminal triangles, which are then refined. This time the new Lepp (t0) computation starts from t0and includes t0,˜ t1which are in turn terminal triangles. The processing of t0concludes after refinement of t0,˜ t1. t0 t2 t3 t1B A (a) t0 t2 t1B A (b) 1 (c) 1 2 3 ~ Figure 1: (a) Lepp(t0) = {t0, t1, t2, t3}and AB is terminal edge; (b) For refining triangle t, a first vertex 1 is added by bisection of the terminal triangles sharing AB. (c) Final triangulation obtained for refining t. 4 2.2. Properties of the 2-dimensional refinement algorithms The non-degeneracy properties of the longest edge algorithms are summarized in lemmas 1 to 4 [21, 43] Lemma 1 (a) The iterative and arbitrary use of the algorithms only produces triangles whose smallest interior angles are always greater than or equal to α/2, where αis the smallest interior angle of the initial triangulation. Also every triangle generated is similar to one of a finite number of reference triangles. (b) Furthermore, for any triangle t generated throughout the refinement process, its smallest angle αtis greater than or equal to α0/2where α0is the smallest angle of the triangle t0of the initial triangulation which contains t. Also t belongs to a finite number of similar triangles associated to t0. Lemma 2 Longest-edge refinement algorithms always terminate in a finite number of steps with the construction of a conforming triangulation. Lemma 3 Any triangulation τgenerated by means of the iterative use of the algorithms satisfies the following smoothness condition: for any pair of sideadjacent triangles t1, t2∈τ(with respective diameters h1, h2) it holds that min(h1,h2) max(h1,h2)≥k > 0, where kdepends on the smallest angle of the initial triangulation. Lemma 4 For any triangulation τ, the global iterative application of the algorithm (the refinement of all the triangles in the preceding iteration) covers, in a monotonically increasing form, the area of τwith quasi-equilateral triangles (with smallest angles ≥30◦). The proof of Lemma 2 is based both on the fact that the refinement propagation is always performed towards bigger triangles in the current mesh, and on the fact that every mesh has bounded smallest angle. The smoothness property of Lemma 3 follows directly from the bound on the smallest angle of part (a) of Lemma 1. Lemma 4 states that the algorithm tends to isolate the worst angles. 2.3. Algorithm costs in two dimensions The practical (adaptive) use of the refinement algorithm, requires of a number of Krefinement iterations which produces a final refined mesh having a not a-priori known number of vertices N=Nref +Nprop, where Nref is the number of triangles iteratively marked for triangle refinement and Nprop is the sum of the number of points introduced by refinement propagation throughout the K iterations [43]. Thus the cost study requires of an amortized cost analysis [54] based on asymptotically studying the behavior of the algorithm throughout the refinement iterations, instead of the classical worst case study [29]. The amortized cost analysis must take into account the fact that, in one iteration the algorithm can introduce propagation points over all the triangles of the mesh (usually when the mesh is small), while for the remaining iterations a very small number of propagation points is introduced by each iteration. 5 To illustrate these ideas consider the mesh of Figure 2(a), where the single longest edge refinement of one triangle (triangle ABC) reaches the complete mesh as shown in Figure 2(b), which corresponds to a worst case behavior for a single refinement step. More importantly, note that after this step, the arbitrary iterative refinement of the triangles of vertex C, produces a very local refinement (1 or 2 vertices by refinement step) that approaches vertex C. Figure 2(c) shows the subtriangulation ADBC obtained after 4 iterative refinement of the triangles of vertex C. CB A (a) (b) (c) A D B B C C AD Figure 2: (a) Initial triangulation; (b) Refinement of triangle ABC reaches the complete mesh in one refinement step; (c) Further refinement around Cis very local (triangulation detail ADBC) Thus the iterative refinement of sets of target triangles introduces both a number Nref of new vertices mandatorily generated to get the required triangle size, and a number Nprop of new propagation vertices which are introduced to keep the geometric quality and the smoothness of the refined mesh, as stated in lemmas 1 and 3. We need: (1) To bound the numbers Nref and Nprop, of new vertices inserted in the mesh; and (2) To study the computational cost of inserting them. The following lemma deals with the second item, which is a non obvious result since an algorithm can be able to select an optimal number of points Kfor point insertion, but the cost of inserting them can be higher than linear as happens with Delaunay point insertion strategy. Lemma 5 Consider any mesh M1obtained by iterative application of the Leppbisection algorithm over an initial mesh M0, implying the introduction of Knew vertices. Then the computational cost of finding and inserting the Kvertices is linear in K, independently of the size of the triangulations and the number of iterations performed. The proof of this lemma follows directly from the fact that each triangle is added one time to the Ordered-List of the Insert-Lepp-Points function, and that the triangle refinement is performed in constant time. 6 Bounds on the number of points Nref and Nprop introduced by the algorithm, as a function of L, the longest interior distance in the geometry D, and the required triangle size δ, are presented in Lemmas 6, 7, 8. We will consider the following two simple problems: (P1) Vertex refinement problem: Iteratively refine the mesh around a vertex Q until the adjacent triangles have longest edge less than or equal to a length parameter δ. (P2) Circle area refinement: Iteratively refine the triangles that intersect a circular region Rcuntil every triangle in Rchas longest edge less than or equal to a length parameter δ. Lemma 6 For solving (P1), a finite number of points N is added to the mesh, by longest edge bisection of pairs of terminal triangles, where N < K(Log(L/δ)), K is a constant such that K= 2π/α, and L is the longest interior distance in the polygonal geometry D measured over the smallest rectangle that contains D. The following lemma states that, for a big enough refinement area and for small δ, the number of propagation vertices is smaller than the number of vertices needed to get the desired mesh refinement. Lemma 7 For solving (P2), finite number of points Niand Neneed to be respectively added in the interior and the exterior of Rcwhere Ni< K1((r δ)2),Ne< K2(r δ)Log(L δ), L is equal to the longest distance from the boundary of Rcto the boundary of D, r is the radius of Rc, and the constants are K1= 4πand K2= 2π. Lemma 8 (a) For solving (P1), the algorithm is linear in N defined in Lemma 6. (b) For solving (P2), the algorithm is linear in (Ni+Ne), the number of points inserted in the mesh. In addition if r >> δ, then the algorithm is linear in Ni. 3. Previous results on parallel refinement algorithms Distributed memory parallel algorithms for the refinement of huge meshes have been studied and used for complex practical applications for the last 20 years [55, 13, 23]. They are mainly based on the partition of basic geometric elements which produce nested refined meshes: algorithms based on the longest edge bisection of triangles / tetrahedra, or algorithms based on quadtree / octree techniques [48, 49]. As far as we know, parallel Delaunay algorithms have not been used in practice, because the serial Delaunay algorithms are less robust and more difficult to parallelize. However, research on the study of different aspects of parallel Delaunay methods, centered in the reuse of serial Delaunay codes, have been performed [2, 3, 15, 16] in recent years. Distributed memory algorithms are based on partitioning the mesh and distributing its pieces to the computers of the cluster. To develop efficient and scalable distributed memory algorithms for mesh refinement, it is necessary to 7 deal with the following issues: (a) to use efficient methods for mesh partitioning and related strategies such as dynamic mesh redistribution to equilibrate the load between processors throughout the computations; (b) to develop wise and efficient strategies for assuring coherent mesh refinement in the submeshes interfaces; and (c) to develop efficient methods for the movement of information between processors, which should be minimized. The development of distributed parallel longest edge algorithms for the parallel refinement of triangulations have been studied and used for complex practical applications related with finite element methods. In a first review paper, Williams [55] recommends the use of parallel 4-triangles longest edge algorithm for the refinement of huge triangulations, for fluid dynamics applications. Later Jones and Plassmann [23, 24, 22] discussed parallel 4-triangles refinement algorithms for distributed memory finite element computations. They find independent sets of triangles to ensure that neighboring triangles are never simultaneously refined on different processors. The independent sets are chosen in parallel by a randomized strategy. A technique that dynamically repartitions the mesh to maintain adequate load distribution, which require to move triangles between processors, is also used. More recently Casta˜nos and Savage [12, 14, 13] proposed a distributed memory parallelization of the original longest edge algorithm in 3-dimensions, and use this method for developing a parallel finite element code; they use a tree data structure to allow the refinement and derefinement of the meshes. Finally, Rivara et al [42] studied a Lepp-based algorithm for uniform refinement of tetrahedral meshes. 4. Multithread Lepp-bisection algorithm In what follows we discuss a shared memory multithread algorithm that takes advantage both of the properties of the Lepp-bisection algorithm, and of the current multicore computers, which have several cores (light processors that perform the reading and instructions execution tasks) and dispose of a great amount of available memory to deal with big meshes. To discuss the multithread Lepp-bisection algorithm, consider one step of the (adaptive finite element) refinement problem introduced in section 2. Given an input triangulation τand a set S⊂τof Ntriangles to be refined, we want to produce a conforming refined triangulation τfsuch that all the triangles of Sare refined in τf. To this end we use a shared memory multicore computer having pphysical cores with p << N. To perform this task each core Pi(i= 1, ..p) is in charge of the processing of an individual triangle tin Sand its associated changing Lepp sequence until the triangle tis refined in the mesh. Once the refinement of tis performed, the associated core will pick up another triangle of Sto continue the refinement task. To design the multithread algorithm, we need to deal with the following synchronization problems: S1 To avoid processing collisions associated to the parallel processing of triangles whose Lepp polygons overlap. 8 6.2. Refinement of an L-shaped domain In order to simulate the adaptive refinement associated to finite element methods, we have considered the L-shaped domain of the Figure 8, with reentrant vertex Bof coordinates (5,5). We have performed refinement by using a circle refinement region of center Band radius r(see Figure 8), for different values of the parameter r. The initial mesh is shown in Figure 8. We have performed iterative refinement of all the triangles that intersect Rcuntil obtaining triangles of size δ= 0.001 (longest edge ≤δ), by starting with the initial mesh of 6 triangles of Figure 8. Note that all the refined meshes only include right isosceles triangles. Tables 14,15,16 show some statistics obtained for the iterative refinement for the last refinement steps, while Tables 17,18,19 summarize the computing times and the values of speedup and efficiency associated to these problems. We can see that the algorithm shows an almost ideal behavior with efficiency higher than 0.86 for all the iterations included. Note that for these problems the refinement concentrates in the interior of Rcwhere the number of triangles refined by propagation remains very low throughout the refinement iterations. This is in complete agreement with the results of Lemmas 7 and 8. A small refined mesh of the L shaped domain is shown in Figure 9. 15 Table 1: Statistics on iterative refinement, 10% smallest triangles, 4 threads case. Mesh Triangles Final Execution Accum Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 19997 390056 444 444 3.62137 17 390056 39006 648271 596 1040 3.62780 16 648271 64827 929824 665 1705 3.71453 17 929824 92982 1244975 856 2561 3.75525 18 1244975 124498 1601045 850 3411 3.75182 19 1601045 160105 2016463 1145 4556 3.71457 21 2016463 201646 2499121 1388 5944 3.67123 23 2499113 249911 3078594 1534 7478 3.61935 23 Table 2: Statistics on iterative refinement, 10% biggest triangles, 4 threads case. Mesh Triangles Final Execution Accum Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 19997 229556 123 123 3.09404 12 229556 22956 255764 66 189 2.94647 14 255764 25576 284528 160 349 2.84749 11 284528 28453 316069 91 440 2.77763 11 316069 31607 350533 85 525 2.72792 12 350533 35053 388300 99 624 2.68983 13 388300 38830 429846 187 811 2.66200 13 429846 42985 475448 113 924 2.64107 13 475448 47545 525769 126 1050 2.62206 12 525769 52577 581251 137 1187 2.60169 10 581251 58125 642463 194 1381 2.58612 10 642463 64246 710053 264 1645 2.57197 10 710053 71005 784723 183 1828 2.56060 10 784723 78472 867121 261 2089 2.55001 11 Table 3: Statistics on iterative refinement, 10% random selection, 4 threads case. Mesh Triangles Final Execution Accum Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 19996 294738 254 254 3.27653 19 294738 29473 420972 303 557 3.19591 16 420972 42097 595838 426 983 3.14954 15 595838 59583 839651 594 1577 3.13082 15 839651 83965 1180479 840 2417 3.11291 14 1180479 118047 1658733 1207 3624 3.10386 15 1658733 165873 2327815 1602 5226 3.09373 13 2327815 232781 3261653 2237 7463 3.08542 14 3261653 326165 4569294 3167 10630 3.07860 17 16 Table 4: Statistics on iterative refinement, 25% smallest triangles, 4 threads case. Mesh Triangles Final Execution Accum. Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 49991 495340 755 755 3.41820 17 495340 123835 1035383 1296 2051 3.28188 16 1035383 258845 1881482 2124 4175 3.26975 17 1881482 470370 3164480 3132 7307 3.27757 17 3164480 791120 5066810 4684 11991 3.26006 19 Table 5: Statistics on iterative refinement, 25% largest triangles, 4 threads case. Mesh Triangles Final Execution Accum. Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 49991 276067 197 197 2.85749 12 276067 69016 356806 250 447 2.71562 12 356806 89201 459687 315 762 2.64974 12 459687 114921 590847 372 1134 2.60073 10 590847 147711 757825 505 1639 2.56403 10 757825 189456 970193 578 2217 2.53998 11 970193 242548 1240472 785 3002 2.52329 10 1240472 310118 1585070 955 3957 2.51098 10 1585070 396267 2023093 1108 5065 2.49905 10 2023093 505773 2580320 1573 6638 2.48278 9 2580320 645080 3288422 1761 8399 2.47567 9 3288422 822105 4188240 2254 10653 2.47121 10 Table 6: Statistics on iterative refinement, 25% random selection, 4 threads case. Mesh Triangles Final Execution Accum. Avrge Longest size to be refined Mesh Size Time [ms] Time Lepp Lepp 199967 49992 389426 532 532 3.16258 15 389426 97357 707114 819 1351 2.99168 13 707114 176779 1255906 1316 2667 2.89928 13 1255906 313977 2205602 2416 5083 2.84503 13 2205602 551401 3847925 4252 9335 2.80838 12 17 Table 7: Execution time and efficiency measures, 10% smallest triangles. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 1508 1007 657 451 1,50 2,30 3,34 0,75 0,77 0,84 22030 1226 800 604 1,66 2,54 3,36 0,83 0,85 0,84 32206 1449 876 759 1,52 2,52 2,91 0,76 0,84 0,73 42475 1449 1006 747 1,71 2,46 3,31 0,85 0,82 0,83 52805 1695 1256 877 1,65 2,23 3,20 0,83 0,74 0,80 63280 1860 1329 1134 1,76 2,47 2,89 0,88 0,82 0,72 73835 2249 1552 1386 1,71 2,47 2,77 0,85 0,82 0,69 84588 2649 1846 1508 1,73 2,49 3,04 0,87 0,83 0,76 Table 8: Execution time and efficiency measures, 10% largest triangles. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 239 133 139 92 1.80 1.72 2.60 0.90 0.57 0.65 2212 120 88 72 1.77 2.41 2.94 0.88 0.80 0.74 3233 130 95 100 1.79 2.45 2.33 0.90 0.82 0.58 4256 143 104 84 1.79 2.46 3.05 0.90 0.82 0.76 5281 157 113 123 1.79 2.49 2.28 0.89 0.83 0.57 6309 171 124 99 1.81 2.49 3.12 0.90 0.83 0.78 7338 188 136 116 1.80 2.49 2.91 0.90 0.83 0.73 8372 225 150 159 1.65 2.48 2.34 0.83 0.83 0.58 9411 230 166 131 1.79 2.48 3.14 0.89 0.83 0.78 10 453 267 197 138 1.70 2.30 3.28 0.85 0.77 0.82 11 500 279 233 216 1.79 2.15 2.31 0.90 0.72 0.58 12 552 307 229 173 1.80 2.41 3.19 0.90 0.80 0.80 18 Table 9: Execution Time and efficiency measures, 10% random selection. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 759 455 307 253 1.67 2.47 3 0.83 0.82 0.75 21085 590 397 294 1.84 2.73 3.69 0.92 0.91 0.92 31650 858 532 423 1.92 3.1 3.9 0.96 1.03 0.98 42438 1066 743 595 2.29 3.28 4.1 1.14 1.09 1.02 53504 1548 1032 839 2.26 3.4 4.18 1.13 1.13 1.04 65034 2153 1439 1295 2.34 3.5 3.89 1.17 1.17 0.97 77191 3070 2023 1583 2.34 3.55 4.54 1.17 1.18 1.14 810650 4147 2870 2221 2.57 3.71 4.8 1.28 1.24 1.2 914628 5903 3978 3190 2.48 3.68 4.59 1.24 1.23 1.15 Table 10: Execution time and efficiency measures, 25% smallest triangles. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 2351 1527 929 755 1,54 2,53 3,11 0,77 0,84 0,78 24236 2424 1626 1296 1,75 2,61 3,27 0,87 0,87 0,82 36691 3805 2593 2124 1,76 2,58 3,15 0,88 0,86 0,79 49988 5828 3979 3132 1,71 2,51 3,19 0,86 0,84 0,80 515008 8729 6306 4684 1,72 2,38 3,20 0,86 0,79 0,80 Table 11: Execution time and efficiency measures, 25% largest triangles. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 595 342 249 197 1,74 2,39 3,02 0,87 0,80 0,76 2628 371 266 250 1,69 2,36 2,51 0,85 0,79 0,63 3799 489 354 315 1,63 2,26 2,54 0,82 0,75 0,63 41017 676 472 372 1,50 2,15 2,73 0,75 0,72 0,68 51291 769 536 505 1,68 2,41 2,56 0,84 0,80 0,64 61639 977 679 578 1,68 2,41 2,84 0,84 0,80 0,71 72081 1232 858 785 1,69 2,43 2,65 0,84 0,81 0,66 82718 1557 1093 955 1,75 2,49 2,85 0,87 0,83 0,71 93370 1999 1392 1108 1,69 2,42 3,04 0,84 0,81 0,76 10 4402 2547 1780 1573 1,73 2,47 2,80 0,86 0,82 0,70 11 5449 3227 2560 1761 1,69 2,13 3,09 0,84 0,71 0,77 12 7091 4144 2877 2254 1,71 2,46 3,15 0,86 0,82 0,79 19 Table 12: Execution time and efficiency measures, 25% random selection. Execution Time (ms) Speed-Up Efficiency Mesh Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 1 1513 980 597 472 1.54 2.53 3.21 0.77 0.84 0.80 22522 1579 987 782 1.60 2.56 3.23 0.80 0.85 0.81 34160 2564 1713 1354 1.62 2.43 3.07 0.81 0.81 0.77 47401 4606 2955 2312 1.61 2.50 3.20 0.80 0.83 0.80 512940 7719 5143 4002 1.68 2.52 3.23 0.84 0.84 0.81 Table 13: Distribution of smallest angles throughout the iterations. Random triangles selection, 10% refinement. Distribution of smallest angles (in %) Mesh Mesh Size 0o −10o10o −20o20o −30o30o −40o40o −50o50o −60o M0 199967 6.43 17.00 24.43 25.51 19.37 7.26 M1 294738 5.55 13.41 21.97 30.05 20.50 8.55 M2 420972 4.76 10.98 20.01 33.76 21.38 9.10 M3 595838 4.11 9.15 18.34 36.85 22.12 9.43 M4 839651 3.58 7.70 16.84 39.55 22.71 9.61 M5 1180479 3.14 6.50 15.54 41.84 23.24 9.77 M6 1658733 2.77 5.54 14.39 43.85 23.63 9.86 M7 2327815 2.46 4.73 13.40 45.57 23.94 9.95 M8 3261653 2.18 4.07 12.54 47.02 24.19 10.04 M9 4700254 1.93 3.52 11.88 48.30 24.29 10.08 20 Figure 5: Distribution of smallest angles, 10% random refinement, 4 threads Figure 6: Initial triangulation Figure 7: 25% random refinement; second refined mesh 21 5 10 5 10 R c B Figure 8: L-shaped domain with refinement region Rc Table 14: Statistics on iterative refinement, L domain, circle refinement region, r=0.3, triangle size δ=0.001, 4 threads case. Refinement Mesh Triangles Final Execution Accum Avrge Longest Iteration size to be refined Mesh Size Time [ms] Time Lepp Lepp 22 39234 37530 77130 222 222 2.05207 14 23 77130 74760 152478 439 661 2.03722 14 24 152478 149094 302370 867 1528 2.02305 14 25 302370 297612 600990 1788 3316 2.01720 14 26 600990 594474 1197294 3557 6873 2.01720 14 Table 15: Statistics on iterative refinement, L domain, circle refinement region, r=0.5, triangle size δ=0.001, 4 threads case. Refinement Mesh Triangles Final Execution Accum Avrge Longest Iteration size to be refined Mesh Size Time [ms] Time Lepp Lepp 22 106452 103662 210738 677 677 2.03220 12 23 210738 206844 418512 1236 1913 2.02137 13 24 418512 413082 833040 2440 4353 2.01533 14 25 833040 825294 1660182 5001 9354 2.00983 11 26 1660182 1649400 3312420 10067 19421 2.00983 11 22 Table 16: Statistics on iterative refinement, L domain, circle refinement region, r=1.2, 24 iterations, 4 threads case. Refinement Mesh Triangles Final Execution Accum Avrge Longest Iteration size to be refined Mesh Size Time [ms] Time Lepp Lepp 20 152406 149094 302298 873 873 2.03703 14 21 302298 297612 600918 1791 2664 2.02295 14 22 600918 594468 1197198 3470 6134 2.01711 14 23 1197198 1188012 2387346 7132 13266 2.01245 14 24 2387346 2374542 4764960 14227 27493 2.00870 14 Table 17: Execution time and efficiency measures, L domain, circle refinement region, r=0.3, δ=0.001. Execution Time (ms) Speed-Up Efficiency Refinement Iteration Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 22 768 439 281 222 1,75 2,73 3,46 0,87 0,91 0,86 23 1559 823 571 439 1,89 2,73 3,55 0,95 0,91 0,89 24 3062 1666 1123 867 1,84 2,73 3,53 0,92 0,91 0,88 25 6343 3427 2300 1788 1,85 2,76 3,55 0,93 0,92 0,89 26 12543 6838 4600 3557 1,83 2,73 3,53 0,92 0,91 0,88 Table 18: Execution time and efficiency measures, L domain, circle refinement region, r = 0.5, δ=0.001 Execution Time (ms) Speed-Up Efficiency Refinement Iteration Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 22 2112 1205 784 677 1,75 2,69 3,12 0,88 0,90 0,78 23 4556 2321 1627 1236 1,96 2,80 3,69 0,98 0,93 0,92 24 8461 4650 3135 2440 1,82 2,70 3,47 0,91 0,90 0,87 25 17367 9593 6396 5001 1,81 2,72 3,47 0,91 0,91 0,87 26 34668 19002 12804 10067 1,82 2,71 3,44 0,91 0,90 0,86 23 Table 19: Execution time and efficiency measures, L domain, circle refinement region, r=1.2, 24 iterations Execution Time (ms) Speed-Up Efficiency Refinement Iteration Serial 2P 3P 4P 2P 3P 4P E2 E3 E4 20 3038 1645 1144 873 1,85 2,66 3,48 0,92 0,89 0,87 21 6213 3425 2312 1791 1,81 2,69 3,47 0,91 0,90 0,87 22 12152 6673 4481 3470 1,82 2,71 3,50 0,91 0,90 0,88 23 24777 13910 9231 7132 1,78 2,68 3,47 0,89 0,89 0,87 24 49303 27328 18363 14227 1,80 2,68 3,47 0,90 0,89 0,87 Figure 9: 24 4. The quality of the refined meshes is only presented for one case (Reviewer #1 and Reviewer #3) 5. The execution times, speed up and efficiency statistics throughout the refinement iterations are presented together in one table. The serial fraction was omitted (Reviewer #3). 6. The hardware is fully described. We have used the physical cores, not the virtual ones, included in the computer. 7. The presentations of Lemmas 1, 4 and 5 were improved by considering the comments of Reviewer #1. 8. The typos and minor errors were corrected (Reviewer #1) Maria-Cecilia Rivara [email protected]