source: inundation/parallel/build_submesh.py @ 2378

Last change on this file since 2378 was 2130, checked in by linda, 19 years ago

Modified the parallel code to agree with the python style files

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