1 | ######################################################### |
---|
2 | # |
---|
3 | # |
---|
4 | # Read in a data file and subdivide the triangle list |
---|
5 | # |
---|
6 | # |
---|
7 | # The final routine, pmesh_divide_metis, does automatic |
---|
8 | # grid partitioning. Once testing has finished on this |
---|
9 | # routine the others should be removed. |
---|
10 | # |
---|
11 | # Authors: Linda Stals and Matthew Hardy, June 2005 |
---|
12 | # Modified: Linda Stals, Nov 2005 |
---|
13 | # Jack Kelly, Nov 2005 |
---|
14 | # |
---|
15 | # |
---|
16 | ######################################################### |
---|
17 | |
---|
18 | from os import sep |
---|
19 | from sys import path |
---|
20 | from math import floor |
---|
21 | |
---|
22 | from Numeric import zeros, Float, Int, reshape, argsort |
---|
23 | |
---|
24 | |
---|
25 | ######################################################### |
---|
26 | # |
---|
27 | # If the triangles list is reordered, the quantities |
---|
28 | # assigned to the triangles must also be reorded. |
---|
29 | # |
---|
30 | # *) quantities contain the quantites in the old ordering |
---|
31 | # *) proc_sum[i] contains the number of triangles in |
---|
32 | # processor i |
---|
33 | # *) tri_index is a map from the old triangle ordering to |
---|
34 | # the new ordering, where the new number for triangle |
---|
35 | # i is proc_sum[tri_index[i][0]]+tri_index[i][1] |
---|
36 | # |
---|
37 | # ------------------------------------------------------- |
---|
38 | # |
---|
39 | # *) The quantaties are returned in the new ordering |
---|
40 | # |
---|
41 | ######################################################### |
---|
42 | |
---|
43 | def reorder(quantities, tri_index, proc_sum): |
---|
44 | |
---|
45 | # Find the number triangles |
---|
46 | |
---|
47 | N = len(tri_index) |
---|
48 | |
---|
49 | # Temporary storage area |
---|
50 | |
---|
51 | index = zeros(N, Int) |
---|
52 | q_reord = {} |
---|
53 | |
---|
54 | # Find the new ordering of the triangles |
---|
55 | |
---|
56 | for i in range(N): |
---|
57 | bin = tri_index[i][0] |
---|
58 | bin_off_set = tri_index[i][1] |
---|
59 | index[i] = proc_sum[bin]+bin_off_set |
---|
60 | |
---|
61 | # Reorder each quantity according to the new ordering |
---|
62 | |
---|
63 | for k in quantities: |
---|
64 | q_reord[k] = zeros((N, 3), Float) |
---|
65 | for i in range(N): |
---|
66 | q_reord[k][index[i]]=quantities[k].vertex_values[i] |
---|
67 | del index |
---|
68 | |
---|
69 | return q_reord |
---|
70 | |
---|
71 | |
---|
72 | ######################################################### |
---|
73 | # |
---|
74 | # Do a simple coordinate based division of the grid |
---|
75 | # |
---|
76 | # *) n_x is the number of subdivisions in the x direction |
---|
77 | # *) n_y is the number of subdivisions in the y direction |
---|
78 | # |
---|
79 | # ------------------------------------------------------- |
---|
80 | # |
---|
81 | # *) The nodes, triangles, boundary, and quantities are |
---|
82 | # returned. triangles_per_proc defines the subdivision. |
---|
83 | # The first triangles_per_proc[0] triangles are assigned |
---|
84 | # to processor 0, the next triangles_per_proc[1] are |
---|
85 | # assigned to processor 1 etc. The boundary and quantites |
---|
86 | # are ordered the same way as the triangles |
---|
87 | # |
---|
88 | ######################################################### |
---|
89 | |
---|
90 | def pmesh_divide(domain, n_x = 1, n_y = 1): |
---|
91 | |
---|
92 | # Find the bounding box |
---|
93 | |
---|
94 | x_coord_min = domain.xy_extent[0] |
---|
95 | x_coord_max = domain.xy_extent[2] |
---|
96 | y_coord_min = domain.xy_extent[1] |
---|
97 | y_coord_max = domain.xy_extent[3] |
---|
98 | |
---|
99 | # Find the size of each sub-box |
---|
100 | |
---|
101 | x_div = (x_coord_max-x_coord_min)/n_x |
---|
102 | y_div = (y_coord_max-y_coord_min)/n_y |
---|
103 | |
---|
104 | # Initialise the lists |
---|
105 | |
---|
106 | tri_list = [] |
---|
107 | triangles_per_proc = [] |
---|
108 | proc_sum = [] |
---|
109 | for i in range(n_x*n_y): |
---|
110 | tri_list.append([]) |
---|
111 | triangles_per_proc.append([]) |
---|
112 | proc_sum.append([]) |
---|
113 | tri_list[i] = [] |
---|
114 | |
---|
115 | # Subdivide the triangles depending on which sub-box they sit |
---|
116 | # in (a triangle sits in the sub-box if its first vectex sits |
---|
117 | # in that sub-box) |
---|
118 | |
---|
119 | tri_index = {} |
---|
120 | N = domain.number_of_elements |
---|
121 | for i in range(N): |
---|
122 | t = domain.triangles[i] |
---|
123 | |
---|
124 | x_coord = domain.centroid_coordinates[i][0] |
---|
125 | bin_x = int(floor((x_coord-x_coord_min)/x_div)) |
---|
126 | if (bin_x == n_x): |
---|
127 | bin_x = n_x-1 |
---|
128 | |
---|
129 | y_coord = domain.centroid_coordinates[i][1] |
---|
130 | bin_y = int(floor((y_coord-y_coord_min)/y_div)) |
---|
131 | if (bin_y == n_y): |
---|
132 | bin_y = n_y-1 |
---|
133 | |
---|
134 | bin = bin_y*n_x + bin_x |
---|
135 | tri_list[bin].append(t) |
---|
136 | tri_index[i] = ([bin, len(tri_list[bin])-1]) |
---|
137 | |
---|
138 | # Find the number of triangles per processor and order the |
---|
139 | # triangle list so that all of the triangles belonging to |
---|
140 | # processor i are listed before those belonging to processor |
---|
141 | # i+1 |
---|
142 | |
---|
143 | triangles = [] |
---|
144 | for i in range(n_x*n_y): |
---|
145 | triangles_per_proc[i] = len(tri_list[i]) |
---|
146 | for t in tri_list[i]: |
---|
147 | triangles.append(t) |
---|
148 | |
---|
149 | # The boundary labels have to changed in accoradance with the |
---|
150 | # new triangle ordering, proc_sum and tri_index help with this |
---|
151 | |
---|
152 | proc_sum[0] = 0 |
---|
153 | for i in range(n_x*n_y-1): |
---|
154 | proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i] |
---|
155 | |
---|
156 | # Relabel the boundary elements to fit in with the new triangle |
---|
157 | # ordering |
---|
158 | |
---|
159 | boundary = {} |
---|
160 | for b in domain.boundary: |
---|
161 | t = tri_index[b[0]] |
---|
162 | boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b] |
---|
163 | |
---|
164 | quantities = reorder(domain.quantities, tri_index, proc_sum) |
---|
165 | |
---|
166 | # Extract the node list |
---|
167 | |
---|
168 | nodes = domain.coordinates.copy() |
---|
169 | ttriangles = zeros((len(triangles), 3), Int) |
---|
170 | for i in range(len(triangles)): |
---|
171 | ttriangles[i] = triangles[i] |
---|
172 | |
---|
173 | return nodes, ttriangles, boundary, triangles_per_proc, quantities |
---|
174 | |
---|
175 | |
---|
176 | ######################################################### |
---|
177 | # |
---|
178 | # Do a coordinate based division of the grid using the |
---|
179 | # centroid coordinate |
---|
180 | # |
---|
181 | # *) n_x is the number of subdivisions in the x direction |
---|
182 | # *) n_y is the number of subdivisions in the y direction |
---|
183 | # |
---|
184 | # ------------------------------------------------------- |
---|
185 | # |
---|
186 | # *) The nodes, triangles, boundary, and quantities are |
---|
187 | # returned. triangles_per_proc defines the subdivision. |
---|
188 | # The first triangles_per_proc[0] triangles are assigned |
---|
189 | # to processor 0, the next triangles_per_proc[1] are |
---|
190 | # assigned to processor 1 etc. The boundary and quantites |
---|
191 | # are ordered the same way as the triangles |
---|
192 | # |
---|
193 | ######################################################### |
---|
194 | def pmesh_divide_steve(domain, n_x = 1, n_y = 1): |
---|
195 | |
---|
196 | # Find the bounding box |
---|
197 | x_coord_min = domain.xy_extent[0] |
---|
198 | x_coord_max = domain.xy_extent[2] |
---|
199 | y_coord_min = domain.xy_extent[1] |
---|
200 | y_coord_max = domain.xy_extent[3] |
---|
201 | |
---|
202 | |
---|
203 | # Find the size of each sub-box |
---|
204 | |
---|
205 | x_div = (x_coord_max-x_coord_min)/n_x |
---|
206 | y_div = (y_coord_max-y_coord_min)/n_y |
---|
207 | |
---|
208 | |
---|
209 | # Initialise the lists |
---|
210 | |
---|
211 | tri_list = [] |
---|
212 | triangles_per_proc = [] |
---|
213 | proc_sum = [] |
---|
214 | for i in range(n_x*n_y): |
---|
215 | tri_list.append([]) |
---|
216 | triangles_per_proc.append([]) |
---|
217 | proc_sum.append([]) |
---|
218 | tri_list[i] = [] |
---|
219 | |
---|
220 | # Subdivide the triangles depending on which sub-box they sit |
---|
221 | # in (a triangle sits in the sub-box if its first vectex sits |
---|
222 | # in that sub-box) |
---|
223 | |
---|
224 | tri_index = {} |
---|
225 | N = domain.number_of_elements |
---|
226 | |
---|
227 | # Sort by x coordinate of centroid |
---|
228 | |
---|
229 | sort_order = argsort(argsort(domain.centroid_coordinates[:,0])) |
---|
230 | |
---|
231 | x_div = float(N)/n_x |
---|
232 | |
---|
233 | for i in range(N): |
---|
234 | t = domain.triangles[i] |
---|
235 | |
---|
236 | bin = int(floor(sort_order[i]/x_div)) |
---|
237 | if (bin == n_x): |
---|
238 | bin = n_x-1 |
---|
239 | |
---|
240 | tri_list[bin].append(t) |
---|
241 | tri_index[i] = ([bin, len(tri_list[bin])-1]) |
---|
242 | |
---|
243 | # Find the number of triangles per processor and order the |
---|
244 | # triangle list so that all of the triangles belonging to |
---|
245 | # processor i are listed before those belonging to processor |
---|
246 | # i+1 |
---|
247 | |
---|
248 | triangles = [] |
---|
249 | for i in range(n_x*n_y): |
---|
250 | triangles_per_proc[i] = len(tri_list[i]) |
---|
251 | for t in tri_list[i]: |
---|
252 | triangles.append(t) |
---|
253 | |
---|
254 | |
---|
255 | # The boundary labels have to changed in accoradance with the |
---|
256 | # new triangle ordering, proc_sum and tri_index help with this |
---|
257 | |
---|
258 | proc_sum[0] = 0 |
---|
259 | for i in range(n_x*n_y-1): |
---|
260 | proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i] |
---|
261 | |
---|
262 | # Relabel the boundary elements to fit in with the new triangle |
---|
263 | # ordering |
---|
264 | |
---|
265 | boundary = {} |
---|
266 | for b in domain.boundary: |
---|
267 | t = tri_index[b[0]] |
---|
268 | boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b] |
---|
269 | |
---|
270 | quantities = reorder(domain.quantities, tri_index, proc_sum) |
---|
271 | |
---|
272 | # Extract the node list |
---|
273 | |
---|
274 | nodes = domain.coordinates.copy() |
---|
275 | |
---|
276 | # Convert the triangle datastructure to be an array type, |
---|
277 | # this helps with the communication |
---|
278 | |
---|
279 | ttriangles = zeros((len(triangles), 3), Int) |
---|
280 | for i in range(len(triangles)): |
---|
281 | ttriangles[i] = triangles[i] |
---|
282 | |
---|
283 | return nodes, ttriangles, boundary, triangles_per_proc, quantities |
---|
284 | |
---|
285 | |
---|
286 | |
---|
287 | |
---|
288 | ### |
---|
289 | # |
---|
290 | # Divide the mesh using a call to metis, through pymetis. |
---|
291 | |
---|
292 | |
---|
293 | path.append('..' + sep + 'pymetis') |
---|
294 | |
---|
295 | from metis import partMeshNodal |
---|
296 | |
---|
297 | def pmesh_divide_metis(domain, n_procs): |
---|
298 | |
---|
299 | # Initialise the lists |
---|
300 | # List, indexed by processor of # triangles. |
---|
301 | |
---|
302 | triangles_per_proc = [] |
---|
303 | |
---|
304 | # List of lists, indexed by processor of vertex numbers |
---|
305 | |
---|
306 | tri_list = [] |
---|
307 | |
---|
308 | # List indexed by processor of cumulative total of triangles allocated. |
---|
309 | |
---|
310 | proc_sum = [] |
---|
311 | for i in range(n_procs): |
---|
312 | tri_list.append([]) |
---|
313 | triangles_per_proc.append(0) |
---|
314 | proc_sum.append([]) |
---|
315 | |
---|
316 | # Prepare variables for the metis call |
---|
317 | |
---|
318 | n_tri = len(domain.triangles) |
---|
319 | if n_procs != 1: #Because metis chokes on it... |
---|
320 | n_vert = len(domain.coordinates) |
---|
321 | t_list = domain.triangles.copy() |
---|
322 | t_list = reshape(t_list, (-1,)) |
---|
323 | |
---|
324 | # The 1 here is for triangular mesh elements. |
---|
325 | |
---|
326 | edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs) |
---|
327 | del edgecut |
---|
328 | del npart |
---|
329 | # print edgecut |
---|
330 | # print npart |
---|
331 | # print epart |
---|
332 | # Assign triangles to processes, according to what metis told us. |
---|
333 | |
---|
334 | # tri_index maps triangle number -> processor, new triangle number |
---|
335 | # (local to the processor) |
---|
336 | |
---|
337 | tri_index = {} |
---|
338 | triangles = [] |
---|
339 | for i in range(n_tri): |
---|
340 | triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1 |
---|
341 | tri_list[epart[i]].append(domain.triangles[i]) |
---|
342 | tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1]) |
---|
343 | |
---|
344 | # print tri_list |
---|
345 | # print triangles_per_proc |
---|
346 | |
---|
347 | # Order the triangle list so that all of the triangles belonging |
---|
348 | # to processor i are listed before those belonging to processor |
---|
349 | # i+1 |
---|
350 | |
---|
351 | for i in range(n_procs): |
---|
352 | for t in tri_list[i]: |
---|
353 | triangles.append(t) |
---|
354 | |
---|
355 | # The boundary labels have to changed in accoradance with the |
---|
356 | # new triangle ordering, proc_sum and tri_index help with this |
---|
357 | |
---|
358 | proc_sum[0] = 0 |
---|
359 | for i in range(n_procs - 1): |
---|
360 | proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i] |
---|
361 | |
---|
362 | # Relabel the boundary elements to fit in with the new triangle |
---|
363 | # ordering |
---|
364 | |
---|
365 | boundary = {} |
---|
366 | for b in domain.boundary: |
---|
367 | t = tri_index[b[0]] |
---|
368 | boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b] |
---|
369 | |
---|
370 | quantities = reorder(domain.quantities, tri_index, proc_sum) |
---|
371 | else: |
---|
372 | boundary = domain.boundary.copy() |
---|
373 | triangles_per_proc[0] = n_tri |
---|
374 | triangles = domain.triangles.copy() |
---|
375 | |
---|
376 | # This is essentially the same as a chunk of code from reorder. |
---|
377 | |
---|
378 | quantities = {} |
---|
379 | for k in domain.quantities: |
---|
380 | quantities[k] = zeros((n_tri, 3), Float) |
---|
381 | for i in range(n_tri): |
---|
382 | quantities[k][i] = domain.quantities[k].vertex_values[i] |
---|
383 | |
---|
384 | # Extract the node list |
---|
385 | |
---|
386 | nodes = domain.coordinates.copy() |
---|
387 | |
---|
388 | # Convert the triangle datastructure to be an array type, |
---|
389 | # this helps with the communication |
---|
390 | |
---|
391 | ttriangles = zeros((len(triangles), 3), Int) |
---|
392 | for i in range(len(triangles)): |
---|
393 | ttriangles[i] = triangles[i] |
---|
394 | |
---|
395 | return nodes, ttriangles, boundary, triangles_per_proc, quantities |
---|