source: anuga_core/source/anuga_parallel/pmesh_divide.py @ 3591

Last change on this file since 3591 was 3591, checked in by linda, 18 years ago

Fixed some bugs in build_submesh

File size: 5.8 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]
[3591]192   
[2850]193    return nodes, ttriangles, boundary, triangles_per_proc, quantities
Note: See TracBrowser for help on using the repository browser.