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

Last change on this file since 3568 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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.pyvolution.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.