source: anuga_core/source/anuga_parallel/build_submesh.py @ 3579

Last change on this file since 3579 was 3579, checked in by ole, 18 years ago

Removed all references to pyvolution in parallel code

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        n = mesh.neighbours[t, 0]
139        if n > 0:
140            if n < tlower or n >= tupper:
141                trianglemap[n] = 1
142        n = mesh.neighbours[t, 1]
143        if n > 0:
144            if n < tlower or n >= tupper:
145                trianglemap[n] = 1
146        n = mesh.neighbours[t, 2]
147        if n > 0:
148            if n < tlower or n >= tupper:
149                trianglemap[n] = 1
150
151    # Find the second layer of boundary triangles
152
153    for t in range(len(trianglemap)):
154        if trianglemap[t]==1:
155            n = mesh.neighbours[t, 0]
156            if n > 0:
157                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
158                    trianglemap[n] = 1
159            n = mesh.neighbours[t, 1]
160            if n > 0:
161                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
162                    trianglemap[n] = 1
163            n = mesh.neighbours[t, 2]
164            if n > 0:
165                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
166                    trianglemap[n] = 1
167
168    # Build the triangle list and make note of the vertices
169
170    nodemap = zeros(ncoord, 'i')
171    fullnodes = submesh["full_nodes"][p]
172
173    subtriangles = []
174    for i in range(len(trianglemap)):
175        if trianglemap[i] == 1:
176            t = list(mesh.triangles[i])
177            nodemap[t[0]] = 1
178            nodemap[t[1]] = 1
179            nodemap[t[2]] = 1
180
181    trilist = reshape(arrayrange(ntriangles),(ntriangles,1))
182    tsubtriangles = concatenate((trilist, mesh.triangles), 1)
183    subtriangles = take(tsubtriangles, nonzero(trianglemap))
184
185    # Keep a record of the triangle vertices, if they are not already there
186
187    subnodes = []
188    for n in fullnodes:
189        nodemap[int(n[0])] = 0
190
191    nodelist = reshape(arrayrange(ncoord),(ncoord,1))
192    tsubnodes = concatenate((nodelist, mesh.coordinates), 1)
193    subnodes = take(tsubnodes, nonzero(nodemap))
194
195    # Clean up before exiting
196
197    del (nodelist)
198    del (trilist)
199    del (tsubnodes)
200    del (nodemap)
201    del (trianglemap)
202
203    # Return the triangles and vertices sitting on the boundary layer
204
205    return subnodes, subtriangles
206
207#########################################################
208#
209# Find the edges of the ghost trianlges that do not
210# have a neighbour in the current cell. These are
211# treated as a special type of boundary edge.
212#
213#  *) Given the ghost triangles in a particular
214# triangle, use the mesh to find its neigbours. If
215# the neighbour is not in the processor set it to
216# be a boundary edge
217#
218#  *) The vertices in the ghost triangles must also
219# be added to the node list for the current processor
220#
221#  *) The boundary edges for the ghost triangles are
222# ignored.
223#
224# -------------------------------------------------------
225#
226#  *) The type assigned to the ghost boundary edges is 'ghost'
227#
228#  *)  The boundary information is returned as a directorier
229# with the key = (triangle id, edge no) and the values
230# assigned to the key is 'ghost'
231#
232#
233#########################################################
234def is_in_processor(ghost_list, tlower, tupper, n):
235    return (n in ghost_list) or (tlower <= n and tupper >= n)
236
237def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
238
239    ghost_list = []
240    subboundary = {}
241       
242    for t in ghosttri:
243        ghost_list.append(t[0])
244
245    for t in ghosttri:
246        n = mesh.neighbours[t[0], 0]
247        if not is_in_processor(ghost_list, tlower, tupper, n):
248            subboundary[t[0], 0] = 'ghost'
249
250        n = mesh.neighbours[t[0], 1]
251        if not is_in_processor(ghost_list, tlower, tupper, n):
252            subboundary[t[0], 1] = 'ghost'
253
254        n = mesh.neighbours[t[0], 2]
255        if not is_in_processor(ghost_list, tlower, tupper, n):
256            subboundary[t[0], 2] = 'ghost'
257   
258    return subboundary
259
260#########################################################
261#
262# The ghost triangles on the current processor will need
263# to get updated information from the neighbouring
264# processor containing the corresponding full triangles.
265#
266#  *) The tri_per_proc is used to determine which
267# processor contains the full node copy.
268#
269# -------------------------------------------------------
270#
271#  *) The ghost communication pattern consists of
272# [global node number, neighbour processor number].
273#
274#########################################################
275
276def ghost_commun_pattern(subtri, p, tri_per_proc):
277
278    # Loop over the ghost triangles
279
280    ghost_commun = zeros((len(subtri), 2), Int)
281
282    for i in range(len(subtri)):
283        global_no = subtri[i][0]
284
285        # Find which processor contains the full triangle
286
287        nproc = len(tri_per_proc)
288        neigh = nproc-1
289        sum = 0
290        for q in range(nproc-1):
291            if (global_no < sum+tri_per_proc[q]):
292                neigh = q
293                break
294            sum = sum+tri_per_proc[q]
295
296        # Keep a copy of the neighbour processor number
297
298        ghost_commun[i] = [global_no, neigh]
299
300    return ghost_commun
301
302#########################################################
303#
304# The full triangles in this processor must communicate
305# updated information to neighbouring processor that
306# contain ghost triangles
307#
308#  *) The ghost communication pattern for all of the
309# processor must be built before calling this processor.
310#
311#  *) The full communication pattern is found by looping
312# through the ghost communication pattern for all of the
313# processors. Recall that this information is stored in
314# the form [global node number, neighbour processor number].
315# The full communication for the neighbour processor is
316# then updated.
317#
318# -------------------------------------------------------
319#
320#  *) The full communication pattern consists of
321# [global id, [p1, p2, ...]], where p1, p2 etc contain
322# a ghost node copy of the triangle global id.
323#
324#########################################################
325
326def full_commun_pattern(submesh, tri_per_proc):
327    tlower = 0
328    nproc = len(tri_per_proc)
329    full_commun = []
330
331    # Loop over the processor
332
333    for p in range(nproc):
334
335        # Loop over the full triangles in the current processor
336        # and build an empty dictionary
337
338        fcommun = {}
339        tupper = tri_per_proc[p]+tlower
340        for i in range(tlower, tupper):
341            fcommun[i] = []
342        full_commun.append(fcommun)
343        tlower = tupper
344
345    # Loop over the processor again
346
347    for p in range(nproc):
348
349        # Loop over the ghost triangles in the current processor,
350        # find which processor contains the corresponding full copy
351        # and note that the processor must send updates to this
352        # processor
353
354        for g in submesh["ghost_commun"][p]:
355            neigh = g[1]
356            full_commun[neigh][g[0]].append(p)
357
358    return full_commun
359
360
361#########################################################
362#
363# Given the non-overlapping grid partition, an extra layer
364# of triangles are included to help with the computations.
365# The triangles in this extra layer are not updated by
366# the processor, their updated values must be sent by the
367# processor containing the original, full, copy of the
368# triangle. The communication pattern that controls these
369# updates must also be built.
370#
371#  *) Assumes that full triangles, nodes etc have already
372# been found and stored in submesh
373#
374#  *) See the documentation for ghost_layer,
375# ghost_commun_pattern and full_commun_pattern
376#
377# -------------------------------------------------------
378#
379#  *) The additional information is added to the submesh
380# dictionary. See the documentation for ghost_layer,
381# ghost_commun_pattern and full_commun_pattern
382#
383#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
384# ghost_commun and full_commun is added to submesh
385#########################################################
386
387def submesh_ghost(submesh, mesh, triangles_per_proc):
388
389    nproc = len(triangles_per_proc)
390    tlower = 0
391    ghost_triangles = []
392    ghost_nodes = []
393    ghost_commun = []
394    ghost_bnd = []
395
396    # Loop over the processors
397
398    for p in range(nproc):
399
400        # Find the full triangles in this processor
401
402        tupper = triangles_per_proc[p]+tlower
403
404        # Build the ghost boundary layer
405
406        [subnodes, subtri] = \
407                   ghost_layer(submesh, mesh, p, tupper, tlower)
408        ghost_triangles.append(subtri)
409        ghost_nodes.append(subnodes)
410
411        # Find the boundary layer formed by the ghost triangles
412       
413        subbnd = ghost_bnd_layer(subtri, tupper, tlower, mesh, p)
414        ghost_bnd.append(subbnd)
415       
416        # Build the communication pattern for the ghost nodes
417
418        gcommun = \
419                ghost_commun_pattern(subtri, p, triangles_per_proc)
420        ghost_commun.append(gcommun)
421
422        # Move to the next processor
423
424        tlower = tupper
425
426    # Record the ghost layer and communication pattern
427
428    submesh["ghost_nodes"] = ghost_nodes
429    submesh["ghost_triangles"] = ghost_triangles
430    submesh["ghost_commun"] = ghost_commun
431    submesh["ghost_boundary"] = ghost_bnd
432   
433    # Build the communication pattern for the full triangles
434
435    full_commun = full_commun_pattern(submesh, triangles_per_proc)
436    submesh["full_commun"] = full_commun
437
438    # Return the submesh
439
440    return submesh
441
442
443#########################################################
444#
445# Certain quantities may be assigned to the triangles,
446# these quantities must be subdivided in the same way
447# as the triangles
448#
449#  *) The quantities are ordered in the same way as the
450# triangles
451#
452# -------------------------------------------------------
453#
454#  *) The quantites attached to the full triangles are
455# stored in full_quan
456#
457#  *) The quantities attached to the ghost triangles are
458# stored in ghost_quan
459#########################################################
460
461def submesh_quantities(submesh, quantities, triangles_per_proc):
462
463    nproc = len(triangles_per_proc)
464
465    lower = 0
466
467    # Build an empty dictionary to hold the quantites
468
469    submesh["full_quan"] = {}
470    submesh["ghost_quan"] = {}
471    for k in quantities:
472        submesh["full_quan"][k] = []
473        submesh["ghost_quan"][k] = []
474
475    # Loop trough the subdomains
476
477    for p in range(nproc):
478        upper =   lower+triangles_per_proc[p]
479
480        # Find the global ID of the ghost triangles
481
482        global_id = []
483        M = len(submesh["ghost_triangles"][p])
484        for j in range(M):
485            global_id.append(submesh["ghost_triangles"][p][j][0])
486
487        # Use the global ID to extract the quantites information from
488        # the full domain
489
490        for k in quantities:
491            submesh["full_quan"][k].append(quantities[k][lower:upper])
492            submesh["ghost_quan"][k].append(zeros( (M,3) , Float))
493            for j in range(M):
494                submesh["ghost_quan"][k][p][j] = \
495                                               quantities[k][global_id[j]]
496
497        lower = upper
498
499    return submesh
500
501#########################################################
502#
503# Build the grid partition on the host.
504#
505#  *) See the documentation for submesh_ghost and
506# submesh_full
507#
508# -------------------------------------------------------
509#
510#  *) A dictionary containing the full_triangles,
511# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
512# ghost_boundary, ghost_commun and full_commun is returned.
513#
514#########################################################
515
516def build_submesh(nodes, triangles, edges, quantities,
517                  triangles_per_proc):
518
519    # Temporarily build the mesh to find the neighbouring
520    # triangles
521
522    mesh = Mesh(nodes, triangles)
523
524    # Subdivide into non-overlapping partitions
525
526    submeshf = submesh_full(nodes, triangles, edges, \
527                            triangles_per_proc)
528
529    # Add any extra ghost boundary layer information
530
531    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
532
533    # Order the quantities information to be the same as the triangle
534    # information
535
536    submesh = submesh_quantities(submeshg, quantities, \
537                                 triangles_per_proc)
538
539    return submesh
540
Note: See TracBrowser for help on using the repository browser.