source: inundation/parallel/pmesh_divide.py @ 3236

Last change on this file since 3236 was 3095, checked in by ole, 19 years ago

Got parallel example RunParallelSwMerimbulaMetis? to run from an arbitrary location using packages relative to PYTHONPATH

File size: 5.9 KB
RevLine 
[2850]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, ArrayType
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# Divide the mesh using a call to metis, through pymetis.
75#
76# -------------------------------------------------------
77#
78# *)  The nodes, triangles, boundary, and quantities are
79# returned. triangles_per_proc defines the subdivision.
80# The first triangles_per_proc[0] triangles are assigned
81# to processor 0, the next triangles_per_proc[1] are
82# assigned to processor 1 etc. The boundary and quantites
83# are ordered the same way as the triangles
84#
85#########################################################
86
[3095]87#path.append('..' + sep + 'pymetis')
88from pymetis.metis import partMeshNodal
[2850]89
90def pmesh_divide_metis(domain, n_procs):
91   
92    # Initialise the lists
93    # List, indexed by processor of # triangles.
94   
95    triangles_per_proc = []
96   
97    # List of lists, indexed by processor of vertex numbers
98   
99    tri_list = []
100   
101    # List indexed by processor of cumulative total of triangles allocated.
102   
103    proc_sum = []
104    for i in range(n_procs):
105        tri_list.append([])
106        triangles_per_proc.append(0)
107        proc_sum.append([])
108
109    # Prepare variables for the metis call
110   
111    n_tri = len(domain.triangles)
112    if n_procs != 1: #Because metis chokes on it...
113        n_vert = len(domain.coordinates)
114        t_list = domain.triangles.copy()
115        t_list = reshape(t_list, (-1,))
116   
117        # The 1 here is for triangular mesh elements.
118        edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
119        # print edgecut
120        # print npart
121        # print epart
122        del edgecut
123        del npart
124
125        # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero
126        # dimensional arrays. Correct this.
127        if type(epart[0]) == ArrayType:
128            epart_new = zeros(len(epart), Int)
129            for i in range(len(epart)):
130                epart_new[i] = epart[i][0]
131            epart = epart_new
132            del epart_new
133        # Assign triangles to processes, according to what metis told us.
134       
135        # tri_index maps triangle number -> processor, new triangle number
136        # (local to the processor)
137       
138        tri_index = {}
139        triangles = []       
140        for i in range(n_tri):
141            triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
142            tri_list[epart[i]].append(domain.triangles[i])
143            tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
144       
145        # Order the triangle list so that all of the triangles belonging
146        # to processor i are listed before those belonging to processor
147        # i+1
148
149        for i in range(n_procs):
150            for t in tri_list[i]:
151                triangles.append(t)
152           
153        # The boundary labels have to changed in accoradance with the
154        # new triangle ordering, proc_sum and tri_index help with this
155
156        proc_sum[0] = 0
157        for i in range(n_procs - 1):
158            proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]
159
160        # Relabel the boundary elements to fit in with the new triangle
161        # ordering
162
163        boundary = {}
164        for b in domain.boundary:
165            t =  tri_index[b[0]]
166            boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]
167
168        quantities = reorder(domain.quantities, tri_index, proc_sum)
169    else:
170        boundary = domain.boundary.copy()
171        triangles_per_proc[0] = n_tri
172        triangles = domain.triangles.copy()
173       
174        # This is essentially the same as a chunk of code from reorder.
175       
176        quantities = {}
177        for k in domain.quantities:
178            quantities[k] = zeros((n_tri, 3), Float)
179            for i in range(n_tri):
180                quantities[k][i] = domain.quantities[k].vertex_values[i]
181       
182    # Extract the node list
183   
184    nodes = domain.coordinates.copy()
185   
186    # Convert the triangle datastructure to be an array type,
187    # this helps with the communication
188
189    ttriangles = zeros((len(triangles), 3), Int)
190    for i in range(len(triangles)):
191        ttriangles[i] = triangles[i]
192       
193    return nodes, ttriangles, boundary, triangles_per_proc, quantities
Note: See TracBrowser for help on using the repository browser.