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, 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 | |
---|
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 | # 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 | |
---|
87 | #path.append('..' + sep + 'pymetis') |
---|
88 | from pymetis.metis import partMeshNodal |
---|
89 | |
---|
90 | def 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 |
---|