source: inundation-numpy-branch/parallel/pmesh_divide.py @ 3031

Last change on this file since 3031 was 2193, checked in by jack, 19 years ago

Corrected potential error in pmesh_divide when examining pmesh_divide_metis
run_parallel_sw_rectangle was asking the visualiser to colour the stage, but was not specifying a colour. It now colours the stage at (0.0, 0.0, 0.8)

File size: 11.6 KB
Line 
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
18from os import sep
19from sys import path
20from math import floor
21
22from 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
43def 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
90def 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#########################################################
194def 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
293path.append('..' + sep + 'pymetis')
294
295from metis import partMeshNodal
296
297def 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        # print edgecut
328        # print npart
329        # print epart
330        del edgecut
331        del npart
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
Note: See TracBrowser for help on using the repository browser.