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 | # Steve Roberts, Aug 2009 (updating to numpy) |
---|
15 | # |
---|
16 | # |
---|
17 | ######################################################### |
---|
18 | |
---|
19 | |
---|
20 | import sys |
---|
21 | from os import sep |
---|
22 | from sys import path |
---|
23 | from math import floor |
---|
24 | |
---|
25 | import numpy as num |
---|
26 | |
---|
27 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
28 | |
---|
29 | ######################################################### |
---|
30 | # |
---|
31 | # If the triangles list is reordered, the quantities |
---|
32 | # assigned to the triangles must also be reorded. |
---|
33 | # |
---|
34 | # *) quantities contain the quantites in the old ordering |
---|
35 | # *) proc_sum[i] contains the number of triangles in |
---|
36 | # processor i |
---|
37 | # *) tri_index is a map from the old triangle ordering to |
---|
38 | # the new ordering, where the new number for triangle |
---|
39 | # i is proc_sum[tri_index[i][0]]+tri_index[i][1] |
---|
40 | # |
---|
41 | # ------------------------------------------------------- |
---|
42 | # |
---|
43 | # *) The quantaties are returned in the new ordering |
---|
44 | # |
---|
45 | ######################################################### |
---|
46 | |
---|
47 | def reorder(quantities, tri_index, proc_sum): |
---|
48 | |
---|
49 | # Find the number triangles |
---|
50 | |
---|
51 | N = len(tri_index) |
---|
52 | |
---|
53 | # Temporary storage area |
---|
54 | |
---|
55 | index = num.zeros(N, num.int) |
---|
56 | q_reord = {} |
---|
57 | |
---|
58 | # Find the new ordering of the triangles |
---|
59 | |
---|
60 | for i in range(N): |
---|
61 | bin = tri_index[i][0] |
---|
62 | bin_off_set = tri_index[i][1] |
---|
63 | index[i] = proc_sum[bin]+bin_off_set |
---|
64 | |
---|
65 | # Reorder each quantity according to the new ordering |
---|
66 | |
---|
67 | for k in quantities: |
---|
68 | q_reord[k] = num.zeros((N, 3), num.float) |
---|
69 | for i in range(N): |
---|
70 | q_reord[k][index[i]]=quantities[k].vertex_values[i] |
---|
71 | del index |
---|
72 | |
---|
73 | return q_reord |
---|
74 | |
---|
75 | |
---|
76 | ######################################################### |
---|
77 | # |
---|
78 | # Divide the mesh using a call to metis, through pymetis. |
---|
79 | # |
---|
80 | # ------------------------------------------------------- |
---|
81 | # |
---|
82 | # *) The nodes, triangles, boundary, and quantities are |
---|
83 | # returned. triangles_per_proc defines the subdivision. |
---|
84 | # The first triangles_per_proc[0] triangles are assigned |
---|
85 | # to processor 0, the next triangles_per_proc[1] are |
---|
86 | # assigned to processor 1 etc. The boundary and quantites |
---|
87 | # are ordered the same way as the triangles |
---|
88 | # |
---|
89 | ######################################################### |
---|
90 | |
---|
91 | #path.append('..' + sep + 'pymetis') |
---|
92 | |
---|
93 | try: |
---|
94 | from pymetis.metis import partMeshNodal |
---|
95 | except ImportError: |
---|
96 | print "***************************************************" |
---|
97 | print " Metis is probably not compiled." |
---|
98 | print " Read \anuga_core\source\pymetis\README" |
---|
99 | print "***************************************************" |
---|
100 | raise ImportError |
---|
101 | |
---|
102 | def pmesh_divide_metis(domain, n_procs): |
---|
103 | |
---|
104 | # Initialise the lists |
---|
105 | # List, indexed by processor of # triangles. |
---|
106 | |
---|
107 | triangles_per_proc = [] |
---|
108 | |
---|
109 | # List of lists, indexed by processor of vertex numbers |
---|
110 | |
---|
111 | tri_list = [] |
---|
112 | |
---|
113 | # List indexed by processor of cumulative total of triangles allocated. |
---|
114 | |
---|
115 | proc_sum = [] |
---|
116 | for i in range(n_procs): |
---|
117 | tri_list.append([]) |
---|
118 | triangles_per_proc.append(0) |
---|
119 | proc_sum.append([]) |
---|
120 | |
---|
121 | # Prepare variables for the metis call |
---|
122 | |
---|
123 | n_tri = len(domain.triangles) |
---|
124 | if n_procs != 1: #Because metis chokes on it... |
---|
125 | n_vert = domain.get_number_of_nodes() |
---|
126 | t_list = domain.triangles.copy() |
---|
127 | t_list = num.reshape(t_list, (-1,)) |
---|
128 | |
---|
129 | # The 1 here is for triangular mesh elements. |
---|
130 | edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs) |
---|
131 | # print edgecut |
---|
132 | # print npart |
---|
133 | # print epart |
---|
134 | del edgecut |
---|
135 | del npart |
---|
136 | |
---|
137 | # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero |
---|
138 | # dimensional arrays. Correct this. |
---|
139 | if type(epart[0]) == num.ndarray: |
---|
140 | epart_new = num.zeros(len(epart), num.int) |
---|
141 | for i in range(len(epart)): |
---|
142 | epart_new[i] = epart[i][0] |
---|
143 | epart = epart_new |
---|
144 | del epart_new |
---|
145 | # Assign triangles to processes, according to what metis told us. |
---|
146 | |
---|
147 | # tri_index maps triangle number -> processor, new triangle number |
---|
148 | # (local to the processor) |
---|
149 | |
---|
150 | tri_index = {} |
---|
151 | triangles = [] |
---|
152 | for i in range(n_tri): |
---|
153 | triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1 |
---|
154 | tri_list[epart[i]].append(domain.triangles[i]) |
---|
155 | tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1]) |
---|
156 | |
---|
157 | # Order the triangle list so that all of the triangles belonging |
---|
158 | # to processor i are listed before those belonging to processor |
---|
159 | # i+1 |
---|
160 | |
---|
161 | for i in range(n_procs): |
---|
162 | for t in tri_list[i]: |
---|
163 | triangles.append(t) |
---|
164 | |
---|
165 | # The boundary labels have to changed in accoradance with the |
---|
166 | # new triangle ordering, proc_sum and tri_index help with this |
---|
167 | |
---|
168 | proc_sum[0] = 0 |
---|
169 | for i in range(n_procs - 1): |
---|
170 | proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i] |
---|
171 | |
---|
172 | # Relabel the boundary elements to fit in with the new triangle |
---|
173 | # ordering |
---|
174 | |
---|
175 | boundary = {} |
---|
176 | for b in domain.boundary: |
---|
177 | t = tri_index[b[0]] |
---|
178 | boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b] |
---|
179 | |
---|
180 | quantities = reorder(domain.quantities, tri_index, proc_sum) |
---|
181 | else: |
---|
182 | boundary = domain.boundary.copy() |
---|
183 | triangles_per_proc[0] = n_tri |
---|
184 | triangles = domain.triangles.copy() |
---|
185 | |
---|
186 | # This is essentially the same as a chunk of code from reorder. |
---|
187 | |
---|
188 | quantities = {} |
---|
189 | for k in domain.quantities: |
---|
190 | quantities[k] = num.zeros((n_tri, 3), num.float) |
---|
191 | for i in range(n_tri): |
---|
192 | quantities[k][i] = domain.quantities[k].vertex_values[i] |
---|
193 | |
---|
194 | # Extract the node list |
---|
195 | |
---|
196 | nodes = domain.get_nodes().copy() |
---|
197 | |
---|
198 | # Convert the triangle datastructure to be an array type, |
---|
199 | # this helps with the communication |
---|
200 | |
---|
201 | ttriangles = num.zeros((len(triangles), 3), num.int) |
---|
202 | for i in range(len(triangles)): |
---|
203 | ttriangles[i] = triangles[i] |
---|
204 | |
---|
205 | return nodes, ttriangles, boundary, triangles_per_proc, quantities |
---|
206 | |
---|
207 | |
---|
208 | ######################################################### |
---|
209 | # |
---|
210 | # Subdivide the domain. This module is primarily |
---|
211 | # responsible for building the ghost layer and |
---|
212 | # communication pattern |
---|
213 | # |
---|
214 | # |
---|
215 | # Author: Linda Stals, June 2005 |
---|
216 | # Modified: Linda Stals, Nov 2005 (optimise python code) |
---|
217 | # Steve Roberts, Aug 2009 (convert to numpy) |
---|
218 | # |
---|
219 | # |
---|
220 | ######################################################### |
---|
221 | |
---|
222 | |
---|
223 | |
---|
224 | ######################################################### |
---|
225 | # |
---|
226 | # Subdivide the triangles into non-overlapping domains. |
---|
227 | # |
---|
228 | # *) The subdivision is controlled by triangles_per_proc. |
---|
229 | # The first triangles_per_proc[0] triangles are assigned |
---|
230 | # to the first processor, the second triangles_per_proc[1] |
---|
231 | # are assigned to the second processor etc. |
---|
232 | # |
---|
233 | # *) nodes, triangles and boundary contains all of the |
---|
234 | # nodes, triangles and boundary tag information for the |
---|
235 | # whole domain. The triangles should be orientated in the |
---|
236 | # correct way and the nodes number consecutively from 0. |
---|
237 | # |
---|
238 | # ------------------------------------------------------- |
---|
239 | # |
---|
240 | # *) A dictionary containing the full_nodes, full_triangles |
---|
241 | # and full_boundary information for each processor is |
---|
242 | # returned. The node information consists of |
---|
243 | # [global_id, x_coord, y_coord]. |
---|
244 | # |
---|
245 | ######################################################### |
---|
246 | |
---|
247 | def submesh_full(nodes, triangles, boundary, triangles_per_proc): |
---|
248 | |
---|
249 | # Initialise |
---|
250 | |
---|
251 | tlower = 0 |
---|
252 | nproc = len(triangles_per_proc) |
---|
253 | nnodes = len(nodes) |
---|
254 | node_list = [] |
---|
255 | triangle_list = [] |
---|
256 | boundary_list = [] |
---|
257 | submesh = {} |
---|
258 | node_range = num.reshape(num.arange(nnodes),(nnodes,1)) |
---|
259 | |
---|
260 | #print node_range |
---|
261 | tsubnodes = num.concatenate((node_range, nodes), 1) |
---|
262 | |
---|
263 | |
---|
264 | # Loop over processors |
---|
265 | |
---|
266 | for p in range(nproc): |
---|
267 | |
---|
268 | # Find triangles on processor p |
---|
269 | |
---|
270 | tupper = triangles_per_proc[p]+tlower |
---|
271 | subtriangles = triangles[tlower:tupper] |
---|
272 | triangle_list.append(subtriangles) |
---|
273 | |
---|
274 | # Find the boundary edges on processor p |
---|
275 | |
---|
276 | subboundary = {} |
---|
277 | for k in boundary: |
---|
278 | if (k[0] >=tlower and k[0] < tupper): |
---|
279 | subboundary[k]=boundary[k] |
---|
280 | boundary_list.append(subboundary) |
---|
281 | |
---|
282 | # Find nodes in processor p |
---|
283 | |
---|
284 | nodemap = num.zeros(nnodes, 'i') |
---|
285 | for t in subtriangles: |
---|
286 | nodemap[t[0]]=1 |
---|
287 | nodemap[t[1]]=1 |
---|
288 | nodemap[t[2]]=1 |
---|
289 | |
---|
290 | |
---|
291 | node_list.append(tsubnodes.take(num.flatnonzero(nodemap),axis=0)) |
---|
292 | |
---|
293 | # Move to the next processor |
---|
294 | |
---|
295 | tlower = tupper |
---|
296 | |
---|
297 | # Put the results in a dictionary |
---|
298 | |
---|
299 | submesh["full_nodes"] = node_list |
---|
300 | submesh["full_triangles"] = triangle_list |
---|
301 | submesh["full_boundary"] = boundary_list |
---|
302 | |
---|
303 | # Clean up before exiting |
---|
304 | |
---|
305 | del (nodemap) |
---|
306 | |
---|
307 | return submesh |
---|
308 | |
---|
309 | |
---|
310 | ######################################################### |
---|
311 | # |
---|
312 | # Build the ghost layer of triangles |
---|
313 | # |
---|
314 | # *) Given the triangle subpartion for the processor |
---|
315 | # build a ghost layer of triangles. The ghost layer |
---|
316 | # consists of two layers of neighbouring triangles. |
---|
317 | # |
---|
318 | # *) The vertices in the ghost triangles must also |
---|
319 | # be added to the node list for the current processor |
---|
320 | # |
---|
321 | # |
---|
322 | # ------------------------------------------------------- |
---|
323 | # |
---|
324 | # *) The extra triangles and nodes are returned. |
---|
325 | # |
---|
326 | # *) The node information consists of |
---|
327 | # [global_id, x_coord, y_coord]. |
---|
328 | # |
---|
329 | # *) The triangle information consists of |
---|
330 | # [triangle number, t], where t = [v1, v2, v3]. |
---|
331 | # |
---|
332 | ######################################################### |
---|
333 | |
---|
334 | def ghost_layer(submesh, mesh, p, tupper, tlower): |
---|
335 | |
---|
336 | ncoord = mesh.number_of_nodes |
---|
337 | ntriangles = mesh.number_of_triangles |
---|
338 | |
---|
339 | # Find the first layer of boundary triangles |
---|
340 | |
---|
341 | trianglemap = num.zeros(ntriangles, 'i') |
---|
342 | for t in range(tlower, tupper): |
---|
343 | |
---|
344 | n = mesh.neighbours[t, 0] |
---|
345 | |
---|
346 | if n >= 0: |
---|
347 | if n < tlower or n >= tupper: |
---|
348 | trianglemap[n] = 1 |
---|
349 | n = mesh.neighbours[t, 1] |
---|
350 | if n >= 0: |
---|
351 | if n < tlower or n >= tupper: |
---|
352 | trianglemap[n] = 1 |
---|
353 | n = mesh.neighbours[t, 2] |
---|
354 | if n >= 0: |
---|
355 | if n < tlower or n >= tupper: |
---|
356 | trianglemap[n] = 1 |
---|
357 | |
---|
358 | # Find the second layer of boundary triangles |
---|
359 | |
---|
360 | for t in range(len(trianglemap)): |
---|
361 | if trianglemap[t]==1: |
---|
362 | n = mesh.neighbours[t, 0] |
---|
363 | if n >= 0: |
---|
364 | if (n < tlower or n >= tupper) and trianglemap[n] == 0: |
---|
365 | trianglemap[n] = 2 |
---|
366 | n = mesh.neighbours[t, 1] |
---|
367 | if n >= 0: |
---|
368 | if (n < tlower or n >= tupper) and trianglemap[n] == 0: |
---|
369 | trianglemap[n] = 2 |
---|
370 | n = mesh.neighbours[t, 2] |
---|
371 | if n >= 0: |
---|
372 | if (n < tlower or n >= tupper) and trianglemap[n] == 0: |
---|
373 | trianglemap[n] = 2 |
---|
374 | |
---|
375 | # Build the triangle list and make note of the vertices |
---|
376 | |
---|
377 | nodemap = num.zeros(ncoord, 'i') |
---|
378 | fullnodes = submesh["full_nodes"][p] |
---|
379 | |
---|
380 | subtriangles = [] |
---|
381 | for i in range(len(trianglemap)): |
---|
382 | if trianglemap[i] != 0: |
---|
383 | t = list(mesh.triangles[i]) |
---|
384 | nodemap[t[0]] = 1 |
---|
385 | nodemap[t[1]] = 1 |
---|
386 | nodemap[t[2]] = 1 |
---|
387 | |
---|
388 | trilist = num.reshape(num.arange(ntriangles),(ntriangles,1)) |
---|
389 | tsubtriangles = num.concatenate((trilist, mesh.triangles), 1) |
---|
390 | subtriangles = tsubtriangles.take(num.flatnonzero(trianglemap),axis=0) |
---|
391 | |
---|
392 | |
---|
393 | # Keep a record of the triangle vertices, if they are not already there |
---|
394 | |
---|
395 | subnodes = [] |
---|
396 | for n in fullnodes: |
---|
397 | nodemap[int(n[0])] = 0 |
---|
398 | |
---|
399 | nodelist = num.reshape(num.arange(ncoord),(ncoord,1)) |
---|
400 | tsubnodes = num.concatenate((nodelist, mesh.get_nodes()), 1) |
---|
401 | subnodes = tsubnodes.take(num.flatnonzero(nodemap),axis=0) |
---|
402 | |
---|
403 | # Clean up before exiting |
---|
404 | |
---|
405 | del (nodelist) |
---|
406 | del (trilist) |
---|
407 | del (tsubnodes) |
---|
408 | del (nodemap) |
---|
409 | del (trianglemap) |
---|
410 | |
---|
411 | # Return the triangles and vertices sitting on the boundary layer |
---|
412 | |
---|
413 | return subnodes, subtriangles |
---|
414 | |
---|
415 | ######################################################### |
---|
416 | # |
---|
417 | # Find the edges of the ghost trianlges that do not |
---|
418 | # have a neighbour in the current cell. These are |
---|
419 | # treated as a special type of boundary edge. |
---|
420 | # |
---|
421 | # *) Given the ghost triangles in a particular |
---|
422 | # triangle, use the mesh to find its neigbours. If |
---|
423 | # the neighbour is not in the processor set it to |
---|
424 | # be a boundary edge |
---|
425 | # |
---|
426 | # *) The vertices in the ghost triangles must also |
---|
427 | # be added to the node list for the current processor |
---|
428 | # |
---|
429 | # *) The boundary edges for the ghost triangles are |
---|
430 | # ignored. |
---|
431 | # |
---|
432 | # ------------------------------------------------------- |
---|
433 | # |
---|
434 | # *) The type assigned to the ghost boundary edges is 'ghost' |
---|
435 | # |
---|
436 | # *) The boundary information is returned as a directorier |
---|
437 | # with the key = (triangle id, edge no) and the values |
---|
438 | # assigned to the key is 'ghost' |
---|
439 | # |
---|
440 | # |
---|
441 | ######################################################### |
---|
442 | def is_in_processor(ghost_list, tlower, tupper, n): |
---|
443 | |
---|
444 | return num.equal(ghost_list,n).any() or (tlower <= n and tupper > n) |
---|
445 | |
---|
446 | |
---|
447 | def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p): |
---|
448 | |
---|
449 | ghost_list = [] |
---|
450 | subboundary = {} |
---|
451 | |
---|
452 | |
---|
453 | for t in ghosttri: |
---|
454 | ghost_list.append(t[0]) |
---|
455 | |
---|
456 | for t in ghosttri: |
---|
457 | |
---|
458 | n = mesh.neighbours[t[0], 0] |
---|
459 | if not is_in_processor(ghost_list, tlower, tupper, n): |
---|
460 | subboundary[t[0], 0] = 'ghost' |
---|
461 | |
---|
462 | n = mesh.neighbours[t[0], 1] |
---|
463 | if not is_in_processor(ghost_list, tlower, tupper, n): |
---|
464 | subboundary[t[0], 1] = 'ghost' |
---|
465 | |
---|
466 | n = mesh.neighbours[t[0], 2] |
---|
467 | if not is_in_processor(ghost_list, tlower, tupper, n): |
---|
468 | subboundary[t[0], 2] = 'ghost' |
---|
469 | |
---|
470 | return subboundary |
---|
471 | |
---|
472 | ######################################################### |
---|
473 | # |
---|
474 | # The ghost triangles on the current processor will need |
---|
475 | # to get updated information from the neighbouring |
---|
476 | # processor containing the corresponding full triangles. |
---|
477 | # |
---|
478 | # *) The tri_per_proc is used to determine which |
---|
479 | # processor contains the full node copy. |
---|
480 | # |
---|
481 | # ------------------------------------------------------- |
---|
482 | # |
---|
483 | # *) The ghost communication pattern consists of |
---|
484 | # [global node number, neighbour processor number]. |
---|
485 | # |
---|
486 | ######################################################### |
---|
487 | |
---|
488 | def ghost_commun_pattern(subtri, p, tri_per_proc): |
---|
489 | |
---|
490 | # Loop over the ghost triangles |
---|
491 | |
---|
492 | ghost_commun = num.zeros((len(subtri), 2), num.int) |
---|
493 | |
---|
494 | for i in range(len(subtri)): |
---|
495 | global_no = subtri[i][0] |
---|
496 | |
---|
497 | # Find which processor contains the full triangle |
---|
498 | |
---|
499 | nproc = len(tri_per_proc) |
---|
500 | neigh = nproc-1 |
---|
501 | sum = 0 |
---|
502 | for q in range(nproc-1): |
---|
503 | if (global_no < sum+tri_per_proc[q]): |
---|
504 | neigh = q |
---|
505 | break |
---|
506 | sum = sum+tri_per_proc[q] |
---|
507 | |
---|
508 | # Keep a copy of the neighbour processor number |
---|
509 | |
---|
510 | ghost_commun[i] = [global_no, neigh] |
---|
511 | |
---|
512 | return ghost_commun |
---|
513 | |
---|
514 | ######################################################### |
---|
515 | # |
---|
516 | # The full triangles in this processor must communicate |
---|
517 | # updated information to neighbouring processor that |
---|
518 | # contain ghost triangles |
---|
519 | # |
---|
520 | # *) The ghost communication pattern for all of the |
---|
521 | # processor must be built before calling this processor. |
---|
522 | # |
---|
523 | # *) The full communication pattern is found by looping |
---|
524 | # through the ghost communication pattern for all of the |
---|
525 | # processors. Recall that this information is stored in |
---|
526 | # the form [global node number, neighbour processor number]. |
---|
527 | # The full communication for the neighbour processor is |
---|
528 | # then updated. |
---|
529 | # |
---|
530 | # ------------------------------------------------------- |
---|
531 | # |
---|
532 | # *) The full communication pattern consists of |
---|
533 | # [global id, [p1, p2, ...]], where p1, p2 etc contain |
---|
534 | # a ghost node copy of the triangle global id. |
---|
535 | # |
---|
536 | ######################################################### |
---|
537 | |
---|
538 | def full_commun_pattern(submesh, tri_per_proc): |
---|
539 | tlower = 0 |
---|
540 | nproc = len(tri_per_proc) |
---|
541 | full_commun = [] |
---|
542 | |
---|
543 | # Loop over the processor |
---|
544 | |
---|
545 | for p in range(nproc): |
---|
546 | |
---|
547 | # Loop over the full triangles in the current processor |
---|
548 | # and build an empty dictionary |
---|
549 | |
---|
550 | fcommun = {} |
---|
551 | tupper = tri_per_proc[p]+tlower |
---|
552 | for i in range(tlower, tupper): |
---|
553 | fcommun[i] = [] |
---|
554 | full_commun.append(fcommun) |
---|
555 | tlower = tupper |
---|
556 | |
---|
557 | # Loop over the processor again |
---|
558 | |
---|
559 | for p in range(nproc): |
---|
560 | |
---|
561 | # Loop over the ghost triangles in the current processor, |
---|
562 | # find which processor contains the corresponding full copy |
---|
563 | # and note that the processor must send updates to this |
---|
564 | # processor |
---|
565 | |
---|
566 | for g in submesh["ghost_commun"][p]: |
---|
567 | neigh = g[1] |
---|
568 | full_commun[neigh][g[0]].append(p) |
---|
569 | |
---|
570 | return full_commun |
---|
571 | |
---|
572 | |
---|
573 | ######################################################### |
---|
574 | # |
---|
575 | # Given the non-overlapping grid partition, an extra layer |
---|
576 | # of triangles are included to help with the computations. |
---|
577 | # The triangles in this extra layer are not updated by |
---|
578 | # the processor, their updated values must be sent by the |
---|
579 | # processor containing the original, full, copy of the |
---|
580 | # triangle. The communication pattern that controls these |
---|
581 | # updates must also be built. |
---|
582 | # |
---|
583 | # *) Assumes that full triangles, nodes etc have already |
---|
584 | # been found and stored in submesh |
---|
585 | # |
---|
586 | # *) See the documentation for ghost_layer, |
---|
587 | # ghost_commun_pattern and full_commun_pattern |
---|
588 | # |
---|
589 | # ------------------------------------------------------- |
---|
590 | # |
---|
591 | # *) The additional information is added to the submesh |
---|
592 | # dictionary. See the documentation for ghost_layer, |
---|
593 | # ghost_commun_pattern and full_commun_pattern |
---|
594 | # |
---|
595 | # *) The ghost_triangles, ghost_nodes, ghost_boundary, |
---|
596 | # ghost_commun and full_commun is added to submesh |
---|
597 | ######################################################### |
---|
598 | |
---|
599 | def submesh_ghost(submesh, mesh, triangles_per_proc): |
---|
600 | |
---|
601 | nproc = len(triangles_per_proc) |
---|
602 | tlower = 0 |
---|
603 | ghost_triangles = [] |
---|
604 | ghost_nodes = [] |
---|
605 | ghost_commun = [] |
---|
606 | ghost_bnd = [] |
---|
607 | |
---|
608 | # Loop over the processors |
---|
609 | |
---|
610 | for p in range(nproc): |
---|
611 | |
---|
612 | # Find the full triangles in this processor |
---|
613 | |
---|
614 | tupper = triangles_per_proc[p]+tlower |
---|
615 | |
---|
616 | # Build the ghost boundary layer |
---|
617 | |
---|
618 | [subnodes, subtri] = \ |
---|
619 | ghost_layer(submesh, mesh, p, tupper, tlower) |
---|
620 | ghost_triangles.append(subtri) |
---|
621 | ghost_nodes.append(subnodes) |
---|
622 | |
---|
623 | |
---|
624 | # Find the boundary layer formed by the ghost triangles |
---|
625 | |
---|
626 | subbnd = ghost_bnd_layer(subtri, tlower, tupper, mesh, p) |
---|
627 | ghost_bnd.append(subbnd) |
---|
628 | |
---|
629 | # Build the communication pattern for the ghost nodes |
---|
630 | |
---|
631 | gcommun = \ |
---|
632 | ghost_commun_pattern(subtri, p, triangles_per_proc) |
---|
633 | ghost_commun.append(gcommun) |
---|
634 | |
---|
635 | # Move to the next processor |
---|
636 | |
---|
637 | tlower = tupper |
---|
638 | |
---|
639 | |
---|
640 | # Record the ghost layer and communication pattern |
---|
641 | |
---|
642 | submesh["ghost_nodes"] = ghost_nodes |
---|
643 | submesh["ghost_triangles"] = ghost_triangles |
---|
644 | submesh["ghost_commun"] = ghost_commun |
---|
645 | submesh["ghost_boundary"] = ghost_bnd |
---|
646 | |
---|
647 | # Build the communication pattern for the full triangles |
---|
648 | |
---|
649 | full_commun = full_commun_pattern(submesh, triangles_per_proc) |
---|
650 | submesh["full_commun"] = full_commun |
---|
651 | |
---|
652 | # Return the submesh |
---|
653 | |
---|
654 | return submesh |
---|
655 | |
---|
656 | |
---|
657 | ######################################################### |
---|
658 | # |
---|
659 | # Certain quantities may be assigned to the triangles, |
---|
660 | # these quantities must be subdivided in the same way |
---|
661 | # as the triangles |
---|
662 | # |
---|
663 | # *) The quantities are ordered in the same way as the |
---|
664 | # triangles |
---|
665 | # |
---|
666 | # ------------------------------------------------------- |
---|
667 | # |
---|
668 | # *) The quantites attached to the full triangles are |
---|
669 | # stored in full_quan |
---|
670 | # |
---|
671 | # *) The quantities attached to the ghost triangles are |
---|
672 | # stored in ghost_quan |
---|
673 | ######################################################### |
---|
674 | |
---|
675 | def submesh_quantities(submesh, quantities, triangles_per_proc): |
---|
676 | |
---|
677 | nproc = len(triangles_per_proc) |
---|
678 | |
---|
679 | lower = 0 |
---|
680 | |
---|
681 | # Build an empty dictionary to hold the quantites |
---|
682 | |
---|
683 | submesh["full_quan"] = {} |
---|
684 | submesh["ghost_quan"] = {} |
---|
685 | for k in quantities: |
---|
686 | submesh["full_quan"][k] = [] |
---|
687 | submesh["ghost_quan"][k] = [] |
---|
688 | |
---|
689 | # Loop trough the subdomains |
---|
690 | |
---|
691 | for p in range(nproc): |
---|
692 | upper = lower+triangles_per_proc[p] |
---|
693 | |
---|
694 | # Find the global ID of the ghost triangles |
---|
695 | |
---|
696 | global_id = [] |
---|
697 | M = len(submesh["ghost_triangles"][p]) |
---|
698 | for j in range(M): |
---|
699 | global_id.append(submesh["ghost_triangles"][p][j][0]) |
---|
700 | |
---|
701 | # Use the global ID to extract the quantites information from |
---|
702 | # the full domain |
---|
703 | |
---|
704 | for k in quantities: |
---|
705 | submesh["full_quan"][k].append(quantities[k][lower:upper]) |
---|
706 | submesh["ghost_quan"][k].append(num.zeros( (M,3) , num.float)) |
---|
707 | for j in range(M): |
---|
708 | submesh["ghost_quan"][k][p][j] = \ |
---|
709 | quantities[k][global_id[j]] |
---|
710 | |
---|
711 | lower = upper |
---|
712 | |
---|
713 | return submesh |
---|
714 | |
---|
715 | ######################################################### |
---|
716 | # |
---|
717 | # Build the grid partition on the host. |
---|
718 | # |
---|
719 | # *) See the documentation for submesh_ghost and |
---|
720 | # submesh_full |
---|
721 | # |
---|
722 | # ------------------------------------------------------- |
---|
723 | # |
---|
724 | # *) A dictionary containing the full_triangles, |
---|
725 | # full_nodes, full_boundary, ghost_triangles, ghost_nodes, |
---|
726 | # ghost_boundary, ghost_commun and full_commun and true boundary polygon is returned. |
---|
727 | # |
---|
728 | ######################################################### |
---|
729 | |
---|
730 | def build_submesh(nodes, triangles, edges, quantities, |
---|
731 | triangles_per_proc): |
---|
732 | |
---|
733 | # Temporarily build the mesh to find the neighbouring |
---|
734 | # triangles and true boundary polygon |
---|
735 | |
---|
736 | mesh = Mesh(nodes, triangles) |
---|
737 | boundary_polygon = mesh.get_boundary_polygon() |
---|
738 | |
---|
739 | |
---|
740 | # Subdivide into non-overlapping partitions |
---|
741 | |
---|
742 | submeshf = submesh_full(nodes, triangles, edges, \ |
---|
743 | triangles_per_proc) |
---|
744 | |
---|
745 | # Add any extra ghost boundary layer information |
---|
746 | |
---|
747 | submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc) |
---|
748 | |
---|
749 | # Order the quantities information to be the same as the triangle |
---|
750 | # information |
---|
751 | |
---|
752 | submesh = submesh_quantities(submeshg, quantities, \ |
---|
753 | triangles_per_proc) |
---|
754 | |
---|
755 | submesh["boundary_polygon"] = boundary_polygon |
---|
756 | return submesh |
---|
757 | |
---|
758 | ######################################################### |
---|
759 | # |
---|
760 | # Given the subdivision of the grid assigned to the |
---|
761 | # current processor convert it into a form that is |
---|
762 | # appropriate for the GA datastructure. |
---|
763 | # |
---|
764 | # The main function of these modules is to change the |
---|
765 | # node numbering. The GA datastructure assumes they |
---|
766 | # are numbered consecutively from 0. |
---|
767 | # |
---|
768 | # The module also changes the communication pattern |
---|
769 | # datastructure into a form needed by parallel_advection |
---|
770 | # |
---|
771 | # Authors: Linda Stals and Matthew Hardy, June 2005 |
---|
772 | # Modified: Linda Stals, Nov 2005 (optimise python code) |
---|
773 | # Steve Roberts, Aug 2009 (updating to numpy) |
---|
774 | # |
---|
775 | # |
---|
776 | ######################################################### |
---|
777 | |
---|
778 | |
---|
779 | ######################################################### |
---|
780 | # Convert the format of the data to that used by ANUGA |
---|
781 | # |
---|
782 | # |
---|
783 | # *) Change the nodes global ID's to an integer value, |
---|
784 | #starting from 0. |
---|
785 | # |
---|
786 | # *) The triangles and boundary edges must also be |
---|
787 | # updated accordingly. |
---|
788 | # |
---|
789 | # ------------------------------------------------------- |
---|
790 | # |
---|
791 | # *) The nodes, triangles and boundary edges defined by |
---|
792 | # the new numbering scheme are returned |
---|
793 | # |
---|
794 | ######################################################### |
---|
795 | |
---|
796 | def build_local_GA(nodes, triangles, boundaries, tri_index): |
---|
797 | |
---|
798 | Nnodes =len(nodes) |
---|
799 | Ntriangles = len(triangles) |
---|
800 | |
---|
801 | # Extract the nodes (using the local ID) |
---|
802 | |
---|
803 | GAnodes = num.take(nodes, (1, 2), 1) |
---|
804 | |
---|
805 | # Build a global ID to local ID mapping |
---|
806 | |
---|
807 | NGlobal = 0 |
---|
808 | for i in range(Nnodes): |
---|
809 | if nodes[i][0] > NGlobal: |
---|
810 | NGlobal = nodes[i][0] |
---|
811 | index = num.zeros(int(NGlobal)+1, num.int) |
---|
812 | num.put(index, num.take(nodes, (0,), 1).astype(num.int), \ |
---|
813 | num.arange(Nnodes)) |
---|
814 | |
---|
815 | # Change the global IDs in the triangles to the local IDs |
---|
816 | |
---|
817 | GAtriangles = num.zeros((Ntriangles, 3), num.int) |
---|
818 | GAtriangles[:,0] = num.take(index, triangles[:,0]) |
---|
819 | GAtriangles[:,1] = num.take(index, triangles[:,1]) |
---|
820 | GAtriangles[:,2] = num.take(index, triangles[:,2]) |
---|
821 | |
---|
822 | # Change the triangle numbering in the boundaries |
---|
823 | |
---|
824 | GAboundaries = {} |
---|
825 | for b in boundaries: |
---|
826 | GAboundaries[tri_index[b[0]], b[1]] = boundaries[b] |
---|
827 | |
---|
828 | del (index) |
---|
829 | |
---|
830 | return GAnodes, GAtriangles, GAboundaries |
---|
831 | |
---|
832 | |
---|
833 | ######################################################### |
---|
834 | # Change the communication format to that needed by the |
---|
835 | # parallel advection file. |
---|
836 | # |
---|
837 | # *) The index contains [global triangle no, |
---|
838 | # local triangle no.] |
---|
839 | # |
---|
840 | # ------------------------------------------------------- |
---|
841 | # |
---|
842 | # *) The ghost_recv and full_send dictionaries are |
---|
843 | # returned. |
---|
844 | # |
---|
845 | # *) ghost_recv dictionary is local id, global id, value |
---|
846 | # |
---|
847 | # *) full_recv dictionary is local id, global id, value |
---|
848 | # |
---|
849 | # *) The information is ordered by the global id. This |
---|
850 | # means that the communication order is predetermined and |
---|
851 | # local and global id do not need to be |
---|
852 | # compared when the information is sent/received. |
---|
853 | # |
---|
854 | ######################################################### |
---|
855 | |
---|
856 | def build_local_commun(index, ghostc, fullc, nproc): |
---|
857 | |
---|
858 | # Initialise |
---|
859 | |
---|
860 | full_send = {} |
---|
861 | ghost_recv = {} |
---|
862 | |
---|
863 | # Build the ghost_recv dictionary (sort the |
---|
864 | # information by the global numbering) |
---|
865 | |
---|
866 | ghostc = num.sort(ghostc, 0) |
---|
867 | |
---|
868 | for c in range(nproc): |
---|
869 | s = ghostc[:,0] |
---|
870 | d = num.compress(num.equal(ghostc[:,1],c), s) |
---|
871 | if len(d) > 0: |
---|
872 | ghost_recv[c] = [0, 0] |
---|
873 | ghost_recv[c][1] = d |
---|
874 | ghost_recv[c][0] = num.take(index, d) |
---|
875 | |
---|
876 | # Build a temporary copy of the full_send dictionary |
---|
877 | # (this version allows the information to be stored |
---|
878 | # by the global numbering) |
---|
879 | |
---|
880 | tmp_send = {} |
---|
881 | for global_id in fullc: |
---|
882 | for i in range(len(fullc[global_id])): |
---|
883 | neigh = fullc[global_id][i] |
---|
884 | if not tmp_send.has_key(neigh): |
---|
885 | tmp_send[neigh] = [] |
---|
886 | tmp_send[neigh].append([global_id, \ |
---|
887 | index[global_id]]) |
---|
888 | |
---|
889 | # Extract the full send information and put it in the form |
---|
890 | # required for the full_send dictionary |
---|
891 | |
---|
892 | for neigh in tmp_send: |
---|
893 | neigh_commun = num.sort(tmp_send[neigh], 0) |
---|
894 | full_send[neigh] = [0, 0] |
---|
895 | full_send[neigh][0] = neigh_commun[:,1] |
---|
896 | full_send[neigh][1] = neigh_commun[:,0] |
---|
897 | |
---|
898 | return ghost_recv, full_send |
---|
899 | |
---|
900 | |
---|
901 | ######################################################### |
---|
902 | # Convert the format of the data to that used by ANUGA |
---|
903 | # |
---|
904 | # |
---|
905 | # *) Change the nodes global ID's to an integer value, |
---|
906 | # starting from 0. The node numbering in the triangles |
---|
907 | # must also be updated to take this into account. |
---|
908 | # |
---|
909 | # *) The triangle number will also change, which affects |
---|
910 | # the boundary tag information and the communication |
---|
911 | # pattern. |
---|
912 | # |
---|
913 | # ------------------------------------------------------- |
---|
914 | # |
---|
915 | # *) The nodes, triangles, boundary edges and communication |
---|
916 | # pattern defined by the new numbering scheme are returned |
---|
917 | # |
---|
918 | ######################################################### |
---|
919 | |
---|
920 | def build_local_mesh(submesh, lower_t, upper_t, nproc): |
---|
921 | |
---|
922 | # Combine the full nodes and ghost nodes |
---|
923 | |
---|
924 | nodes = num.concatenate((submesh["full_nodes"], \ |
---|
925 | submesh["ghost_nodes"])) |
---|
926 | |
---|
927 | # Combine the full triangles and ghost triangles |
---|
928 | |
---|
929 | gtri = num.take(submesh["ghost_triangles"],(1, 2, 3),1) |
---|
930 | triangles = num.concatenate((submesh["full_triangles"], gtri)) |
---|
931 | |
---|
932 | # Combine the full boundaries and ghost boundaries |
---|
933 | |
---|
934 | boundaries = submesh["full_boundary"] |
---|
935 | for b in submesh["ghost_boundary"]: |
---|
936 | boundaries[b]=submesh["ghost_boundary"][b] |
---|
937 | |
---|
938 | # Make note of the new triangle numbers, including the ghost |
---|
939 | # triangles |
---|
940 | |
---|
941 | NGlobal = upper_t |
---|
942 | for i in range(len(submesh["ghost_triangles"])): |
---|
943 | id = submesh["ghost_triangles"][i][0] |
---|
944 | if id > NGlobal: |
---|
945 | NGlobal = id |
---|
946 | index = num.zeros(int(NGlobal)+1, num.int) |
---|
947 | index[lower_t:upper_t]=num.arange(upper_t-lower_t) |
---|
948 | for i in range(len(submesh["ghost_triangles"])): |
---|
949 | index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t |
---|
950 | |
---|
951 | # Change the node numbering (and update the numbering in the |
---|
952 | # triangles) |
---|
953 | |
---|
954 | [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index) |
---|
955 | |
---|
956 | # Extract the local quantities |
---|
957 | |
---|
958 | quantities ={} |
---|
959 | for k in submesh["full_quan"]: |
---|
960 | Nf = len(submesh["full_quan"][k]) |
---|
961 | Ng = len(submesh["ghost_quan"][k]) |
---|
962 | quantities[k] = num.zeros((Nf+Ng, 3), num.float) |
---|
963 | quantities[k][0:Nf] = submesh["full_quan"][k] |
---|
964 | quantities[k][Nf:Nf+Ng] = submesh["ghost_quan"][k] |
---|
965 | |
---|
966 | # Change the communication pattern into a form needed by |
---|
967 | # the parallel_adv |
---|
968 | |
---|
969 | gcommun = submesh["ghost_commun"] |
---|
970 | fcommun = submesh["full_commun"] |
---|
971 | [ghost_rec, full_send] = \ |
---|
972 | build_local_commun(index, gcommun, fcommun, nproc) |
---|
973 | |
---|
974 | # Clean up before exiting |
---|
975 | |
---|
976 | del(index) |
---|
977 | |
---|
978 | return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \ |
---|
979 | full_send |
---|
980 | |
---|
981 | |
---|
982 | ######################################################### |
---|
983 | # |
---|
984 | # Handle the communication between the host machine |
---|
985 | # (processor 0) and the processors. The host machine is |
---|
986 | # responsible for the doing the initial grid partitioning. |
---|
987 | # |
---|
988 | # The routines given below should be moved to the |
---|
989 | # build_submesh.py and build_local.py file to allow |
---|
990 | # overlapping of communication and computation. |
---|
991 | # This should be done after more debugging. |
---|
992 | # |
---|
993 | # |
---|
994 | # Author: Linda Stals, June 2005 |
---|
995 | # Modified: Linda Stals, Nov 2005 (optimise python code) |
---|
996 | # Steve Roberts, Aug 2009 (update to numpy) |
---|
997 | # |
---|
998 | # |
---|
999 | ######################################################### |
---|
1000 | |
---|
1001 | |
---|
1002 | ######################################################### |
---|
1003 | # |
---|
1004 | # Send the submesh to processor p. |
---|
1005 | # |
---|
1006 | # *) The order and form is strongly coupled with |
---|
1007 | # rec_submesh. |
---|
1008 | # |
---|
1009 | # ------------------------------------------------------- |
---|
1010 | # |
---|
1011 | # *) All of the information has been sent to processor p. |
---|
1012 | # |
---|
1013 | ######################################################### |
---|
1014 | |
---|
1015 | def send_submesh(submesh, triangles_per_proc, p, verbose=True): |
---|
1016 | |
---|
1017 | import pypar |
---|
1018 | |
---|
1019 | myid = pypar.rank() |
---|
1020 | |
---|
1021 | if verbose: print 'process %d sending submesh to process %d' %(myid, p) |
---|
1022 | |
---|
1023 | # build and send the tagmap for the boundary conditions |
---|
1024 | |
---|
1025 | tagmap = {} |
---|
1026 | counter = 1 |
---|
1027 | for b in submesh["full_boundary"][p]: |
---|
1028 | bkey = submesh["full_boundary"][p][b] |
---|
1029 | if not tagmap.has_key(bkey): |
---|
1030 | tagmap[bkey] = counter |
---|
1031 | counter = counter+1 |
---|
1032 | for b in submesh["ghost_boundary"][p]: |
---|
1033 | bkey = submesh["ghost_boundary"][p][b] |
---|
1034 | if not tagmap.has_key(bkey): |
---|
1035 | tagmap[bkey] = counter |
---|
1036 | counter = counter+1 |
---|
1037 | |
---|
1038 | pypar.send(tagmap, p) |
---|
1039 | |
---|
1040 | # send the quantities key information |
---|
1041 | |
---|
1042 | pypar.send(submesh["full_quan"].keys(), p) |
---|
1043 | |
---|
1044 | # send the number of triangles per processor |
---|
1045 | |
---|
1046 | pypar.send(triangles_per_proc, p) |
---|
1047 | |
---|
1048 | # compress full_commun |
---|
1049 | |
---|
1050 | flat_full_commun = [] |
---|
1051 | |
---|
1052 | for c in submesh["full_commun"][p]: |
---|
1053 | for i in range(len(submesh["full_commun"][p][c])): |
---|
1054 | flat_full_commun.append([c,submesh["full_commun"][p][c][i]]) |
---|
1055 | |
---|
1056 | # send the array sizes so memory can be allocated |
---|
1057 | |
---|
1058 | setup_array = num.zeros((9,),num.int) |
---|
1059 | setup_array[0] = len(submesh["full_nodes"][p]) |
---|
1060 | setup_array[1] = len(submesh["ghost_nodes"][p]) |
---|
1061 | setup_array[2] = len(submesh["full_triangles"][p]) |
---|
1062 | setup_array[3] = len(submesh["ghost_triangles"][p]) |
---|
1063 | setup_array[4] = len(submesh["full_boundary"][p]) |
---|
1064 | setup_array[5] = len(submesh["ghost_boundary"][p]) |
---|
1065 | setup_array[6] = len(submesh["ghost_commun"][p]) |
---|
1066 | setup_array[7] = len(flat_full_commun) |
---|
1067 | setup_array[8] = len(submesh["full_quan"]) |
---|
1068 | |
---|
1069 | pypar.send(num.array(setup_array, num.int), p) |
---|
1070 | |
---|
1071 | # send the nodes |
---|
1072 | |
---|
1073 | pypar.send(num.array(submesh["full_nodes"][p], num.float), p) |
---|
1074 | pypar.send(num.array(submesh["ghost_nodes"][p], num.float),p) |
---|
1075 | |
---|
1076 | # send the triangles |
---|
1077 | |
---|
1078 | pypar.send(num.array(submesh["full_triangles"][p], num.int), p) |
---|
1079 | pypar.send(num.array(submesh["ghost_triangles"][p], num.int), p) |
---|
1080 | |
---|
1081 | # send the boundary |
---|
1082 | |
---|
1083 | bc = [] |
---|
1084 | for b in submesh["full_boundary"][p]: |
---|
1085 | bc.append([b[0], b[1], tagmap[submesh["full_boundary"][p][b]]]) |
---|
1086 | |
---|
1087 | |
---|
1088 | pypar.send(num.array(bc, num.int), p) |
---|
1089 | |
---|
1090 | bc = [] |
---|
1091 | for b in submesh["ghost_boundary"][p]: |
---|
1092 | bc.append([b[0], b[1], tagmap[submesh["ghost_boundary"][p][b]]]) |
---|
1093 | |
---|
1094 | pypar.send(num.array(bc, num.int), p) |
---|
1095 | |
---|
1096 | # send the communication pattern |
---|
1097 | |
---|
1098 | pypar.send(submesh["ghost_commun"][p], p) |
---|
1099 | |
---|
1100 | pypar.send(num.array(flat_full_commun, num.int), p) |
---|
1101 | |
---|
1102 | # send the quantities |
---|
1103 | |
---|
1104 | for k in submesh["full_quan"]: |
---|
1105 | pypar.send(num.array(submesh["full_quan"][k][p], num.float), p) |
---|
1106 | |
---|
1107 | for k in submesh["ghost_quan"]: |
---|
1108 | pypar.send(num.array(submesh["ghost_quan"][k][p], num.float),p) |
---|
1109 | |
---|
1110 | |
---|
1111 | ######################################################### |
---|
1112 | # |
---|
1113 | # Receive the submesh from processor p. |
---|
1114 | # |
---|
1115 | # *) The order and form is strongly coupled with |
---|
1116 | # send_submesh. |
---|
1117 | # |
---|
1118 | # ------------------------------------------------------- |
---|
1119 | # |
---|
1120 | # *) All of the information has been received by the |
---|
1121 | # processor p and passed into build_local. |
---|
1122 | # |
---|
1123 | # *) The information is returned in a form needed by the |
---|
1124 | # GA datastructure. |
---|
1125 | # |
---|
1126 | ######################################################### |
---|
1127 | |
---|
1128 | def rec_submesh_flat(p, verbose=True): |
---|
1129 | |
---|
1130 | import pypar |
---|
1131 | |
---|
1132 | numproc = pypar.size() |
---|
1133 | myid = pypar.rank() |
---|
1134 | |
---|
1135 | submesh_cell = {} |
---|
1136 | |
---|
1137 | if verbose: print 'process %d receiving submesh from process %d' %(myid, p) |
---|
1138 | |
---|
1139 | # receive the tagmap for the boundary conditions |
---|
1140 | |
---|
1141 | tagmap = pypar.receive(p) |
---|
1142 | |
---|
1143 | itagmap = {} |
---|
1144 | for t in tagmap: |
---|
1145 | itagmap[tagmap[t]]=t |
---|
1146 | |
---|
1147 | # receive the quantities key information |
---|
1148 | |
---|
1149 | qkeys = pypar.receive(p) |
---|
1150 | |
---|
1151 | # receive the number of triangles per processor |
---|
1152 | |
---|
1153 | triangles_per_proc = pypar.receive(p) |
---|
1154 | |
---|
1155 | # recieve information about the array sizes |
---|
1156 | |
---|
1157 | setup_array = pypar.receive(p) |
---|
1158 | |
---|
1159 | no_full_nodes = setup_array[0] |
---|
1160 | no_ghost_nodes = setup_array[1] |
---|
1161 | no_full_triangles = setup_array[2] |
---|
1162 | no_ghost_triangles = setup_array[3] |
---|
1163 | no_full_boundary = setup_array[4] |
---|
1164 | no_ghost_boundary = setup_array[5] |
---|
1165 | no_ghost_commun = setup_array[6] |
---|
1166 | no_full_commun = setup_array[7] |
---|
1167 | no_quantities = setup_array[8] |
---|
1168 | |
---|
1169 | # receive the full nodes |
---|
1170 | |
---|
1171 | submesh_cell["full_nodes"] = pypar.receive(p) |
---|
1172 | |
---|
1173 | # receive the ghost nodes |
---|
1174 | |
---|
1175 | submesh_cell["ghost_nodes"] = pypar.receive(p) |
---|
1176 | |
---|
1177 | # receive the full triangles |
---|
1178 | |
---|
1179 | submesh_cell["full_triangles"] = pypar.receive(p) |
---|
1180 | |
---|
1181 | # receive the ghost triangles |
---|
1182 | |
---|
1183 | submesh_cell["ghost_triangles"] = pypar.receive(p) |
---|
1184 | |
---|
1185 | # receive the full boundary |
---|
1186 | |
---|
1187 | bnd_c = pypar.receive(p) |
---|
1188 | |
---|
1189 | submesh_cell["full_boundary"] = {} |
---|
1190 | for b in bnd_c: |
---|
1191 | submesh_cell["full_boundary"][b[0],b[1]]=itagmap[b[2]] |
---|
1192 | |
---|
1193 | # receive the ghost boundary |
---|
1194 | |
---|
1195 | bnd_c = pypar.receive(p) |
---|
1196 | |
---|
1197 | submesh_cell["ghost_boundary"] = {} |
---|
1198 | for b in bnd_c: |
---|
1199 | submesh_cell["ghost_boundary"][b[0],b[1]]=itagmap[b[2]] |
---|
1200 | |
---|
1201 | # receive the ghost communication pattern |
---|
1202 | |
---|
1203 | submesh_cell["ghost_commun"] = pypar.receive(p) |
---|
1204 | |
---|
1205 | # receive the full communication pattern |
---|
1206 | |
---|
1207 | full_commun = pypar.receive(p) |
---|
1208 | |
---|
1209 | submesh_cell["full_commun"] = {} |
---|
1210 | for c in full_commun: |
---|
1211 | submesh_cell["full_commun"][c[0]] = [] |
---|
1212 | for c in full_commun: |
---|
1213 | submesh_cell["full_commun"][c[0]].append(c[1]) |
---|
1214 | |
---|
1215 | # receive the quantities |
---|
1216 | |
---|
1217 | submesh_cell["full_quan"]={} |
---|
1218 | |
---|
1219 | for i in range(no_quantities): |
---|
1220 | tmp = pypar.receive(p) |
---|
1221 | submesh_cell["full_quan"][qkeys[i]]=num.zeros((no_full_triangles,3), num.float) |
---|
1222 | submesh_cell["full_quan"][qkeys[i]][:] = tmp[:] |
---|
1223 | |
---|
1224 | submesh_cell["ghost_quan"]={} |
---|
1225 | for i in range(no_quantities): |
---|
1226 | tmp = pypar.receive(p) |
---|
1227 | submesh_cell["ghost_quan"][qkeys[i]]= num.zeros((no_ghost_triangles,3), num.float) |
---|
1228 | submesh_cell["ghost_quan"][qkeys[i]][:] = tmp[:] |
---|
1229 | |
---|
1230 | return submesh_cell, triangles_per_proc,\ |
---|
1231 | no_full_nodes, no_full_triangles |
---|
1232 | |
---|
1233 | |
---|
1234 | |
---|
1235 | ######################################################### |
---|
1236 | # |
---|
1237 | # Receive the submesh from processor p. |
---|
1238 | # |
---|
1239 | # *) The order and form is strongly coupled with |
---|
1240 | # send_submesh. |
---|
1241 | # |
---|
1242 | # ------------------------------------------------------- |
---|
1243 | # |
---|
1244 | # *) All of the information has been received by the |
---|
1245 | # processor p and passed into build_local. |
---|
1246 | # |
---|
1247 | # *) The information is returned in a form needed by the |
---|
1248 | # GA datastructure. |
---|
1249 | # |
---|
1250 | ######################################################### |
---|
1251 | |
---|
1252 | def rec_submesh(p, verbose=True): |
---|
1253 | |
---|
1254 | import pypar |
---|
1255 | |
---|
1256 | numproc = pypar.size() |
---|
1257 | myid = pypar.rank() |
---|
1258 | |
---|
1259 | [submesh_cell, triangles_per_proc,\ |
---|
1260 | number_of_full_nodes, number_of_full_triangles] = rec_submesh_flat(p,verbose) |
---|
1261 | |
---|
1262 | # find the full triangles assigned to this processor |
---|
1263 | |
---|
1264 | lower_t = 0 |
---|
1265 | for i in range(myid): |
---|
1266 | lower_t = lower_t+triangles_per_proc[i] |
---|
1267 | upper_t = lower_t+triangles_per_proc[myid] |
---|
1268 | |
---|
1269 | # convert the information into a form needed by the GA |
---|
1270 | # datastructure |
---|
1271 | |
---|
1272 | [GAnodes, GAtriangles, boundary, quantities, ghost_rec, full_send] = \ |
---|
1273 | build_local_mesh(submesh_cell, lower_t, upper_t, \ |
---|
1274 | numproc) |
---|
1275 | |
---|
1276 | return GAnodes, GAtriangles, boundary, quantities,\ |
---|
1277 | ghost_rec, full_send,\ |
---|
1278 | number_of_full_nodes, number_of_full_triangles |
---|
1279 | |
---|
1280 | |
---|
1281 | ######################################################### |
---|
1282 | # |
---|
1283 | # Extract the submesh that will belong to the |
---|
1284 | # "host processor" (i.e. processor zero) |
---|
1285 | # |
---|
1286 | # *) See the documentation for build_submesh |
---|
1287 | # |
---|
1288 | # ------------------------------------------------------- |
---|
1289 | # |
---|
1290 | # *) A dictionary containing the full_triangles, |
---|
1291 | # full_nodes, full_boundary, ghost_triangles, ghost_nodes, |
---|
1292 | # ghost_boundary, ghost_commun and full_commun belonging |
---|
1293 | # to processor zero are returned. |
---|
1294 | # |
---|
1295 | ######################################################### |
---|
1296 | def extract_hostmesh(submesh, triangles_per_proc): |
---|
1297 | |
---|
1298 | |
---|
1299 | submesh_cell = {} |
---|
1300 | submesh_cell["full_nodes"] = submesh["full_nodes"][0] |
---|
1301 | submesh_cell["ghost_nodes"] = submesh["ghost_nodes"][0] |
---|
1302 | submesh_cell["full_triangles"] = submesh["full_triangles"][0] |
---|
1303 | submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0] |
---|
1304 | submesh_cell["full_boundary"] = submesh["full_boundary"][0] |
---|
1305 | submesh_cell["ghost_boundary"] = submesh["ghost_boundary"][0] |
---|
1306 | submesh_cell["ghost_commun"] = submesh["ghost_commun"][0] |
---|
1307 | submesh_cell["full_commun"] = submesh["full_commun"][0] |
---|
1308 | submesh_cell["full_quan"] ={} |
---|
1309 | submesh_cell["ghost_quan"]={} |
---|
1310 | for k in submesh["full_quan"]: |
---|
1311 | submesh_cell["full_quan"][k] = submesh["full_quan"][k][0] |
---|
1312 | submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0] |
---|
1313 | |
---|
1314 | numprocs = len(triangles_per_proc) |
---|
1315 | points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \ |
---|
1316 | build_local_mesh(submesh_cell, 0, triangles_per_proc[0], numprocs) |
---|
1317 | |
---|
1318 | |
---|
1319 | return points, vertices, boundary, quantities, ghost_recv_dict, \ |
---|
1320 | full_send_dict |
---|
1321 | |
---|
1322 | |
---|
1323 | |
---|
1324 | |
---|