source: anuga_core/source/anuga_parallel/build_submesh.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: 15.9 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 anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
20
21
22
23#########################################################
24#
25# Subdivide the triangles into non-overlapping domains.
26#
27#  *)  The subdivision is controlled by triangles_per_proc.
28# The first triangles_per_proc[0] triangles are assigned
29# to the first processor, the second triangles_per_proc[1]
30# are assigned to the second processor etc.
31#
32#  *) nodes, triangles and boundary contains all of the
33# nodes, triangles and boundary tag information for the
34# whole domain. The triangles should be orientated in the
35# correct way and the nodes number consecutively from 0.
36#
37# -------------------------------------------------------
38#
39#  *) A dictionary containing the full_nodes, full_triangles
40# and full_boundary information for each processor is
41# returned. The node information consists of
42# [global_id, x_coord, y_coord].
43#
44#########################################################
45
46def submesh_full(nodes, triangles, boundary, triangles_per_proc):
47
48    # Initialise
49
50    tlower = 0
51    nproc = len(triangles_per_proc)
52    nnodes = len(nodes)
53    node_list = []
54    triangle_list = []
55    boundary_list = []
56    submesh = {}
57    node_range = reshape(arrayrange(nnodes),(nnodes,1))
58    tsubnodes = concatenate((node_range, nodes), 1)
59
60    # Loop over processors
61
62    for p in range(nproc):
63
64        # Find triangles on processor p
65
66        tupper = triangles_per_proc[p]+tlower
67        subtriangles = triangles[tlower:tupper]
68        triangle_list.append(subtriangles)
69
70        # Find the boundary edges on processor p
71
72        subboundary = {}
73        for k in boundary:
74            if (k[0] >=tlower and k[0] < tupper):
75                subboundary[k]=boundary[k]
76        boundary_list.append(subboundary)
77
78        # Find nodes in processor p
79
80        nodemap = zeros(nnodes, 'i')
81        for t in subtriangles:
82            nodemap[t[0]]=1
83            nodemap[t[1]]=1
84            nodemap[t[2]]=1
85
86        node_list.append(take(tsubnodes,nonzero(nodemap)))
87
88        # Move to the next processor
89
90        tlower = tupper
91
92    # Put the results in a dictionary
93
94    submesh["full_nodes"] = node_list
95    submesh["full_triangles"] = triangle_list
96    submesh["full_boundary"] = boundary_list
97
98    # Clean up before exiting
99
100    del (nodemap)
101
102    return submesh
103
104
105#########################################################
106#
107# Build the ghost layer of triangles
108#
109#  *) Given the triangle subpartion for the processor
110# build a ghost layer of triangles. The ghost layer
111# consists of two layers of neighbouring triangles.
112#
113#  *) The vertices in the ghost triangles must also
114# be added to the node list for the current processor
115#
116#
117# -------------------------------------------------------
118#
119#  *) The extra triangles and nodes are returned.
120#
121#  *)  The node information consists of
122# [global_id, x_coord, y_coord].
123#
124#  *) The triangle information consists of
125# [triangle number, t], where t = [v1, v2, v3].
126#
127#########################################################
128
129def ghost_layer(submesh, mesh, p, tupper, tlower):
130
131    ncoord = len(mesh.coordinates)
132    ntriangles = len(mesh.triangles)
133
134    # Find the first layer of boundary triangles
135
136    trianglemap = zeros(ntriangles, 'i')
137    for t in range(tlower, tupper):
138       
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] = 2
160            n = mesh.neighbours[t, 1]
161            if n >= 0:
162                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
163                    trianglemap[n] = 2
164            n = mesh.neighbours[t, 2]
165            if n >= 0:
166                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
167                    trianglemap[n] = 2
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] != 0:
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# Find the edges of the ghost trianlges that do not
211# have a neighbour in the current cell. These are
212# treated as a special type of boundary edge.
213#
214#  *) Given the ghost triangles in a particular
215# triangle, use the mesh to find its neigbours. If
216# the neighbour is not in the processor set it to
217# be a boundary edge
218#
219#  *) The vertices in the ghost triangles must also
220# be added to the node list for the current processor
221#
222#  *) The boundary edges for the ghost triangles are
223# ignored.
224#
225# -------------------------------------------------------
226#
227#  *) The type assigned to the ghost boundary edges is 'ghost'
228#
229#  *)  The boundary information is returned as a directorier
230# with the key = (triangle id, edge no) and the values
231# assigned to the key is 'ghost'
232#
233#
234#########################################################
235def is_in_processor(ghost_list, tlower, tupper, n):
236    return (n in ghost_list) or (tlower <= n and tupper > n)
237
238def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
239
240    ghost_list = []
241    subboundary = {}
242       
243    for t in ghosttri:
244        ghost_list.append(t[0])
245   
246    for t in ghosttri:
247        n = mesh.neighbours[t[0], 0]
248        if not is_in_processor(ghost_list, tlower, tupper, n):
249            subboundary[t[0], 0] = 'ghost'
250
251        n = mesh.neighbours[t[0], 1]
252        if not is_in_processor(ghost_list, tlower, tupper, n):
253            subboundary[t[0], 1] = 'ghost'
254
255        n = mesh.neighbours[t[0], 2]
256        if not is_in_processor(ghost_list, tlower, tupper, n):
257            subboundary[t[0], 2] = 'ghost'
258           
259    return subboundary
260
261#########################################################
262#
263# The ghost triangles on the current processor will need
264# to get updated information from the neighbouring
265# processor containing the corresponding full triangles.
266#
267#  *) The tri_per_proc is used to determine which
268# processor contains the full node copy.
269#
270# -------------------------------------------------------
271#
272#  *) The ghost communication pattern consists of
273# [global node number, neighbour processor number].
274#
275#########################################################
276
277def ghost_commun_pattern(subtri, p, tri_per_proc):
278
279    # Loop over the ghost triangles
280
281    ghost_commun = zeros((len(subtri), 2), Int)
282
283    for i in range(len(subtri)):
284        global_no = subtri[i][0]
285
286        # Find which processor contains the full triangle
287
288        nproc = len(tri_per_proc)
289        neigh = nproc-1
290        sum = 0
291        for q in range(nproc-1):
292            if (global_no < sum+tri_per_proc[q]):
293                neigh = q
294                break
295            sum = sum+tri_per_proc[q]
296
297        # Keep a copy of the neighbour processor number
298
299        ghost_commun[i] = [global_no, neigh]
300
301    return ghost_commun
302
303#########################################################
304#
305# The full triangles in this processor must communicate
306# updated information to neighbouring processor that
307# contain ghost triangles
308#
309#  *) The ghost communication pattern for all of the
310# processor must be built before calling this processor.
311#
312#  *) The full communication pattern is found by looping
313# through the ghost communication pattern for all of the
314# processors. Recall that this information is stored in
315# the form [global node number, neighbour processor number].
316# The full communication for the neighbour processor is
317# then updated.
318#
319# -------------------------------------------------------
320#
321#  *) The full communication pattern consists of
322# [global id, [p1, p2, ...]], where p1, p2 etc contain
323# a ghost node copy of the triangle global id.
324#
325#########################################################
326
327def full_commun_pattern(submesh, tri_per_proc):
328    tlower = 0
329    nproc = len(tri_per_proc)
330    full_commun = []
331
332    # Loop over the processor
333
334    for p in range(nproc):
335
336        # Loop over the full triangles in the current processor
337        # and build an empty dictionary
338
339        fcommun = {}
340        tupper = tri_per_proc[p]+tlower
341        for i in range(tlower, tupper):
342            fcommun[i] = []
343        full_commun.append(fcommun)
344        tlower = tupper
345
346    # Loop over the processor again
347
348    for p in range(nproc):
349
350        # Loop over the ghost triangles in the current processor,
351        # find which processor contains the corresponding full copy
352        # and note that the processor must send updates to this
353        # processor
354
355        for g in submesh["ghost_commun"][p]:
356            neigh = g[1]
357            full_commun[neigh][g[0]].append(p)
358
359    return full_commun
360
361
362#########################################################
363#
364# Given the non-overlapping grid partition, an extra layer
365# of triangles are included to help with the computations.
366# The triangles in this extra layer are not updated by
367# the processor, their updated values must be sent by the
368# processor containing the original, full, copy of the
369# triangle. The communication pattern that controls these
370# updates must also be built.
371#
372#  *) Assumes that full triangles, nodes etc have already
373# been found and stored in submesh
374#
375#  *) See the documentation for ghost_layer,
376# ghost_commun_pattern and full_commun_pattern
377#
378# -------------------------------------------------------
379#
380#  *) The additional information is added to the submesh
381# dictionary. See the documentation for ghost_layer,
382# ghost_commun_pattern and full_commun_pattern
383#
384#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
385# ghost_commun and full_commun is added to submesh
386#########################################################
387
388def submesh_ghost(submesh, mesh, triangles_per_proc):
389
390    nproc = len(triangles_per_proc)
391    tlower = 0
392    ghost_triangles = []
393    ghost_nodes = []
394    ghost_commun = []
395    ghost_bnd = []
396
397    # Loop over the processors
398
399    for p in range(nproc):
400
401        # Find the full triangles in this processor
402
403        tupper = triangles_per_proc[p]+tlower
404
405        # Build the ghost boundary layer
406
407        [subnodes, subtri] = \
408                   ghost_layer(submesh, mesh, p, tupper, tlower)
409        ghost_triangles.append(subtri)
410        ghost_nodes.append(subnodes)
411
412        # Find the boundary layer formed by the ghost triangles
413       
414        subbnd = ghost_bnd_layer(subtri, tlower, tupper, mesh, p)
415        ghost_bnd.append(subbnd)
416       
417        # Build the communication pattern for the ghost nodes
418
419        gcommun = \
420                ghost_commun_pattern(subtri, p, triangles_per_proc)
421        ghost_commun.append(gcommun)
422
423        # Move to the next processor
424
425        tlower = tupper
426
427
428    # Record the ghost layer and communication pattern
429
430    submesh["ghost_nodes"] = ghost_nodes
431    submesh["ghost_triangles"] = ghost_triangles
432    submesh["ghost_commun"] = ghost_commun
433    submesh["ghost_boundary"] = ghost_bnd
434   
435    # Build the communication pattern for the full triangles
436
437    full_commun = full_commun_pattern(submesh, triangles_per_proc)
438    submesh["full_commun"] = full_commun
439
440    # Return the submesh
441
442    return submesh
443
444
445#########################################################
446#
447# Certain quantities may be assigned to the triangles,
448# these quantities must be subdivided in the same way
449# as the triangles
450#
451#  *) The quantities are ordered in the same way as the
452# triangles
453#
454# -------------------------------------------------------
455#
456#  *) The quantites attached to the full triangles are
457# stored in full_quan
458#
459#  *) The quantities attached to the ghost triangles are
460# stored in ghost_quan
461#########################################################
462
463def submesh_quantities(submesh, quantities, triangles_per_proc):
464
465    nproc = len(triangles_per_proc)
466
467    lower = 0
468
469    # Build an empty dictionary to hold the quantites
470
471    submesh["full_quan"] = {}
472    submesh["ghost_quan"] = {}
473    for k in quantities:
474        submesh["full_quan"][k] = []
475        submesh["ghost_quan"][k] = []
476
477    # Loop trough the subdomains
478
479    for p in range(nproc):
480        upper =   lower+triangles_per_proc[p]
481
482        # Find the global ID of the ghost triangles
483
484        global_id = []
485        M = len(submesh["ghost_triangles"][p])
486        for j in range(M):
487            global_id.append(submesh["ghost_triangles"][p][j][0])
488
489        # Use the global ID to extract the quantites information from
490        # the full domain
491
492        for k in quantities:
493            submesh["full_quan"][k].append(quantities[k][lower:upper])
494            submesh["ghost_quan"][k].append(zeros( (M,3) , Float))
495            for j in range(M):
496                submesh["ghost_quan"][k][p][j] = \
497                                               quantities[k][global_id[j]]
498
499        lower = upper
500
501    return submesh
502
503#########################################################
504#
505# Build the grid partition on the host.
506#
507#  *) See the documentation for submesh_ghost and
508# submesh_full
509#
510# -------------------------------------------------------
511#
512#  *) A dictionary containing the full_triangles,
513# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
514# ghost_boundary, ghost_commun and full_commun is returned.
515#
516#########################################################
517
518def build_submesh(nodes, triangles, edges, quantities,
519                  triangles_per_proc):
520
521    # Temporarily build the mesh to find the neighbouring
522    # triangles
523
524    mesh = Mesh(nodes, triangles)
525
526    # Subdivide into non-overlapping partitions
527
528    submeshf = submesh_full(nodes, triangles, edges, \
529                            triangles_per_proc)
530   
531    # Add any extra ghost boundary layer information
532
533    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
534
535    # Order the quantities information to be the same as the triangle
536    # information
537
538    submesh = submesh_quantities(submeshg, quantities, \
539                                 triangles_per_proc)
540
541    return submesh
542
Note: See TracBrowser for help on using the repository browser.