[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 | |
---|
| 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 | |
---|
[3095] | 87 | #path.append('..' + sep + 'pymetis') |
---|
| 88 | from pymetis.metis import partMeshNodal |
---|
[2850] | 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 |
---|