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

Last change on this file since 3818 was 3818, checked in by ole, 17 years ago

Work on parallel test and attempt to get TRUE bounding polygon for each submesh

File size: 16.1 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 and true boundary polygon 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 and true boundary polygon
523
524    mesh = Mesh(nodes, triangles)
525    boundary_polygon = mesh.get_boundary_polygon()
526   
527
528    # Subdivide into non-overlapping partitions
529
530    submeshf = submesh_full(nodes, triangles, edges, \
531                            triangles_per_proc)
532   
533    # Add any extra ghost boundary layer information
534
535    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
536
537    # Order the quantities information to be the same as the triangle
538    # information
539
540    submesh = submesh_quantities(submeshg, quantities, \
541                                 triangles_per_proc)
542
543    submesh["boundary_polygon"] = boundary_polygon
544    return submesh
545
Note: See TracBrowser for help on using the repository browser.