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

Last change on this file since 7400 was 7400, checked in by steve, 15 years ago

Commit a working copy of numpy version of build_commun

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