source: inundation/parallel/build_submesh.py @ 2130

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

Modified the parallel code to agree with the python style files

File size: 15.2 KB
Line 
1#########################################################
2#
3# Subdivide the GA domain. This module is primarily
4# responsible for building the ghost layer and
5# communication pattern
6#
7#
8#  Author: Linda Stals, June 2005
9#  Modified: Linda Stals, Nov 2005 (optimise python code)
10#
11#
12#########################################################
13
14import sys
15
16from Numeric import zeros, Float, Int, concatenate, \
17     reshape, arrayrange, take, nonzero
18
19from mesh import Mesh
20
21
22#########################################################
23#
24# Subdivide the triangles into non-overlapping domains.
25#
26#  *)  The subdivision is controlled by triangles_per_proc.
27# The first triangles_per_proc[0] triangles are assigned
28# to the first processor, the second triangles_per_proc[1]
29# are assigned to the second processor etc.
30#
31#  *) nodes, triangles and boundary contains all of the
32# nodes, triangles and boundary tag information for the
33# whole domain. The triangles should be orientated in the
34# correct way and the nodes number consecutively from 0.
35#
36# -------------------------------------------------------
37#
38#  *) A dictionary containing the full_nodes, full_triangles
39# and full_boundary information for each processor is
40# returned. The node information consists of
41# [global_id, x_coord, y_coord].
42#
43#########################################################
44
45def submesh_full(nodes, triangles, boundary, triangles_per_proc):
46
47    # Initialise
48
49    tlower = 0
50    nproc = len(triangles_per_proc)
51    nnodes = len(nodes)
52    node_list = []
53    triangle_list = []
54    boundary_list = []
55    submesh = {}
56    node_range = reshape(arrayrange(nnodes),(nnodes,1))
57    tsubnodes = concatenate((node_range, nodes), 1)
58
59    # Loop over processors
60
61    for p in range(nproc):
62
63        # Find triangles on processor p
64
65        tupper = triangles_per_proc[p]+tlower
66        subtriangles = triangles[tlower:tupper]
67        triangle_list.append(subtriangles)
68
69        # Find the boundary edges on processor p
70
71        subboundary = {}
72        for k in boundary:
73            if (k[0] >=tlower and k[0] < tupper):
74                subboundary[k]=boundary[k]
75        boundary_list.append(subboundary)
76
77        # Find nodes in processor p
78
79        nodemap = zeros(nnodes, 'i')
80        for t in subtriangles:
81            nodemap[t[0]]=1
82            nodemap[t[1]]=1
83            nodemap[t[2]]=1
84
85        node_list.append(take(tsubnodes,nonzero(nodemap)))
86
87        # Move to the next processor
88
89        tlower = tupper
90
91    # Put the results in a dictionary
92
93    submesh["full_nodes"] = node_list
94    submesh["full_triangles"] = triangle_list
95    submesh["full_boundary"] = boundary_list
96
97    # Clean up before exiting
98
99    del (nodemap)
100
101    return submesh
102
103
104#########################################################
105#
106# Build the ghost layer of triangles
107#
108#  *) Given the triangle subpartion for the processor
109# build a ghost layer of triangles. The ghost layer
110# consists of two layers of neighbouring triangles.
111#
112#  *) The vertices in the ghost triangles must also
113# be added to the node list for the current processor
114#
115#  *) The boundary edges for the ghost triangles are
116# ignored.
117#
118# -------------------------------------------------------
119#
120#  *) The extra triangles and nodes are returned.
121#
122#  *)  The node information consists of
123# [global_id, x_coord, y_coord].
124#
125#  *) The triangle information consists of
126# [triangle number, t], where t = [v1, v2, v3].
127#
128#########################################################
129
130def ghost_layer(submesh, mesh, p, tupper, tlower):
131
132    ncoord = len(mesh.coordinates)
133    ntriangles = len(mesh.triangles)
134
135    # Find the first layer of boundary triangles
136
137    trianglemap = zeros(ntriangles, 'i')
138    for t in range(tlower, tupper):
139        n = mesh.neighbours[t, 0]
140        if n > 0:
141            if n < tlower or n >= tupper:
142                trianglemap[n] = 1
143        n = mesh.neighbours[t, 1]
144        if n > 0:
145            if n < tlower or n >= tupper:
146                trianglemap[n] = 1
147        n = mesh.neighbours[t, 2]
148        if n > 0:
149            if n < tlower or n >= tupper:
150                trianglemap[n] = 1
151
152    # Find the second layer of boundary triangles
153
154    for t in range(len(trianglemap)):
155        if trianglemap[t]==1:
156            n = mesh.neighbours[t, 0]
157            if n > 0:
158                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
159                    trianglemap[n] = 1
160            n = mesh.neighbours[t, 1]
161            if n > 0:
162                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
163                    trianglemap[n] = 1
164            n = mesh.neighbours[t, 2]
165            if n > 0:
166                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
167                    trianglemap[n] = 1
168
169    # Build the triangle list and make note of the vertices
170
171    nodemap = zeros(ncoord, 'i')
172    fullnodes = submesh["full_nodes"][p]
173
174    subtriangles = []
175    for i in range(len(trianglemap)):
176        if trianglemap[i] == 1:
177            t = list(mesh.triangles[i])
178            nodemap[t[0]] = 1
179            nodemap[t[1]] = 1
180            nodemap[t[2]] = 1
181
182    trilist = reshape(arrayrange(ntriangles),(ntriangles,1))
183    tsubtriangles = concatenate((trilist, mesh.triangles), 1)
184    subtriangles = take(tsubtriangles, nonzero(trianglemap))
185
186    # Keep a record of the triangle vertices, if they are not already there
187
188    subnodes = []
189    for n in fullnodes:
190        nodemap[int(n[0])] = 0
191
192    nodelist = reshape(arrayrange(ncoord),(ncoord,1))
193    tsubnodes = concatenate((nodelist, mesh.coordinates), 1)
194    subnodes = take(tsubnodes, nonzero(nodemap))
195
196    # Clean up before exiting
197
198    del (nodelist)
199    del (trilist)
200    del (tsubnodes)
201    del (nodemap)
202    del (trianglemap)
203
204    # Return the triangles and vertices sitting on the boundary layer
205
206    return subnodes, subtriangles
207
208
209#########################################################
210#
211# The ghost triangles on the current processor will need
212# to get updated information from the neighbouring
213# processor containing the corresponding full triangles.
214#
215#  *) The tri_per_proc is used to determine which
216# processor contains the full node copy.
217#
218# -------------------------------------------------------
219#
220#  *) The ghost communication pattern consists of
221# [global node number, neighbour processor number].
222#
223#########################################################
224
225def ghost_commun_pattern(subtri, p, tri_per_proc):
226
227    # Loop over the ghost triangles
228
229    ghost_commun = zeros((len(subtri), 2), Int)
230
231    for i in range(len(subtri)):
232        global_no = subtri[i][0]
233
234        # Find which processor contains the full triangle
235
236        nproc = len(tri_per_proc)
237        neigh = nproc-1
238        sum = 0
239        for q in range(nproc-1):
240            if (global_no < sum+tri_per_proc[q]):
241                neigh = q
242                break
243            sum = sum+tri_per_proc[q]
244
245        # Keep a copy of the neighbour processor number
246
247        ghost_commun[i] = [global_no, neigh]
248
249    return ghost_commun
250
251#########################################################
252#
253# The full triangles in this processor must communicate
254# updated information to neighbouring processor that
255# contain ghost triangles
256#
257#  *) The ghost communication pattern for all of the
258# processor must be built before calling this processor.
259#
260#  *) The full communication pattern is found by looping
261# through the ghost communication pattern for all of the
262# processors. Recall that this information is stored in
263# the form [global node number, neighbour processor number].
264# The full communication for the neighbour processor is
265# then updated.
266#
267# -------------------------------------------------------
268#
269#  *) The full communication pattern consists of
270# [global id, [p1, p2, ...]], where p1, p2 etc contain
271# a ghost node copy of the triangle global id.
272#
273#########################################################
274
275def full_commun_pattern(submesh, tri_per_proc):
276    tlower = 0
277    nproc = len(tri_per_proc)
278    full_commun = []
279
280    # Loop over the processor
281
282    for p in range(nproc):
283
284        # Loop over the full triangles in the current processor
285        # and build an empty dictionary
286
287        fcommun = {}
288        tupper = tri_per_proc[p]+tlower
289        for i in range(tlower, tupper):
290            fcommun[i] = []
291        full_commun.append(fcommun)
292        tlower = tupper
293
294    # Loop over the processor again
295
296    for p in range(nproc):
297
298        # Loop over the ghost triangles in the current processor,
299        # find which processor contains the corresponding full copy
300        # and note that the processor must send updates to this
301        # processor
302
303        for g in submesh["ghost_commun"][p]:
304            neigh = g[1]
305            full_commun[neigh][g[0]].append(p)
306
307    return full_commun
308
309
310#########################################################
311#
312# Given the non-overlapping grid partition, an extra layer
313# of triangles are included to help with the computations.
314# The triangles in this extra layer are not updated by
315# the processor, their updated values must be sent by the
316# processor containing the original, full, copy of the
317# triangle. The communication pattern that controls these
318# updates must also be built.
319#
320#  *) See the documentation for ghost_layer,
321# ghost_commun_pattern and full_commun_pattern
322#
323# -------------------------------------------------------
324#
325#  *) The additional information is added to the submesh
326# dictionary. See the documentation for ghost_layer,
327# ghost_commun_pattern and full_commun_pattern
328#
329#  *) The ghost_triangles, ghost_nodes, ghost_commun and
330# full_commun is added to submesh
331#########################################################
332
333def submesh_ghost(submesh, mesh, triangles_per_proc):
334
335    nproc = len(triangles_per_proc)
336    tlower = 0
337    ghost_triangles = []
338    ghost_nodes = []
339    ghost_commun = []
340
341    # Loop over the processors
342
343    for p in range(nproc):
344
345        # Find the full triangles in this processor
346
347        tupper = triangles_per_proc[p]+tlower
348
349        # Build the ghost boundary layer
350
351        [subnodes, subtri] = \
352                   ghost_layer(submesh, mesh, p, tupper, tlower)
353        ghost_triangles.append(subtri)
354        ghost_nodes.append(subnodes)
355
356        # Build the communication pattern for the ghost nodes
357
358        gcommun = \
359                ghost_commun_pattern(subtri, p, triangles_per_proc)
360        ghost_commun.append(gcommun)
361
362        # Move to the next processor
363
364        tlower = tupper
365
366    # Record the ghost layer and communication pattern
367
368    submesh["ghost_nodes"] = ghost_nodes
369    submesh["ghost_triangles"] = ghost_triangles
370    submesh["ghost_commun"] = ghost_commun
371
372    # Build the communication pattern for the full triangles
373
374    full_commun = full_commun_pattern(submesh, triangles_per_proc)
375    submesh["full_commun"] = full_commun
376
377    # Return the submesh
378
379    return submesh
380
381
382#########################################################
383#
384# Certain quantities may be assigned to the triangles,
385# these quantities must be subdivided in the same way
386# as the triangles
387#
388#  *) The quantities are ordered in the same way as the
389# triangles
390#
391# -------------------------------------------------------
392#
393#  *) The quantites attached to the full triangles are
394# stored in full_quan
395#
396#  *) The quantities attached to the ghost triangles are
397# stored in ghost_quan
398#########################################################
399
400def submesh_quantities(submesh, quantities, triangles_per_proc):
401
402    nproc = len(triangles_per_proc)
403
404    lower = 0
405
406    # Build an empty dictionary to hold the quantites
407
408    submesh["full_quan"] = {}
409    submesh["ghost_quan"] = {}
410    for k in quantities:
411        submesh["full_quan"][k] = []
412        submesh["ghost_quan"][k] = []
413
414    # Loop trough the subdomains
415
416    for p in range(nproc):
417        upper =   lower+triangles_per_proc[p]
418
419        # Find the global ID of the ghost triangles
420
421        global_id = []
422        M = len(submesh["ghost_triangles"][p])
423        for j in range(M):
424            global_id.append(submesh["ghost_triangles"][p][j][0])
425
426        # Use the global ID to extract the quantites information from
427        # the full domain
428
429        for k in quantities:
430            submesh["full_quan"][k].append(quantities[k][lower:upper])
431            submesh["ghost_quan"][k].append(zeros( (M,3) , Float))
432            for j in range(M):
433                submesh["ghost_quan"][k][p][j] = \
434                                               quantities[k][global_id[j]]
435
436        lower = upper
437
438    return submesh
439
440#########################################################
441#
442# Build the grid partition on the host.
443#
444#  *) See the documentation for submesh_ghost and
445# submesh_full
446#
447# -------------------------------------------------------
448#
449#  *) A dictionary containing the full_triangles,
450# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
451# ghost_commun and full_commun is returned.
452#
453#########################################################
454
455def build_submesh(nodes, triangles, edges, quantities,
456                  triangles_per_proc):
457
458    # Temporarily build the mesh to find the neighbouring
459    # triangles
460
461    mesh = Mesh(nodes, triangles)
462
463    # Subdivide into non-overlapping partitions
464
465    submeshf = submesh_full(nodes, triangles, edges, \
466                            triangles_per_proc)
467
468    # Add any extra ghost boundary layer information
469
470    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
471
472    # Order the quantities information to be the same as the triangle
473    # information
474
475    submesh = submesh_quantities(submeshg, quantities, \
476                                 triangles_per_proc)
477
478    return submesh
479
480#########################################################
481#
482# Extract the submesh that will belong to the
483# "host processor" (i.e. processor zero)
484#
485#  *) See the documentation for build_submesh
486#
487# -------------------------------------------------------
488#
489#  *) A dictionary containing the full_triangles,
490# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
491# ghost_commun and full_commun belonging to processor zero
492# are returned.
493#
494#########################################################
495def extract_hostmesh(submesh):
496
497    submesh_cell = {}
498    submesh_cell["full_nodes"] = submesh["full_nodes"][0]
499    submesh_cell["ghost_nodes"] = submesh["ghost_nodes"][0]
500    submesh_cell["full_triangles"] = submesh["full_triangles"][0]
501    submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0]
502    submesh_cell["full_boundary"] = submesh["full_boundary"][0]
503    submesh_cell["ghost_commun"] = submesh["ghost_commun"][0]
504    submesh_cell["full_commun"] = submesh["full_commun"][0]
505    submesh_cell["full_quan"] ={}
506    submesh_cell["ghost_quan"]={}
507    for k in submesh["full_quan"]:
508        submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
509        submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
510
511    return submesh_cell
512
Note: See TracBrowser for help on using the repository browser.