source: anuga_core/source/anuga_parallel/distribute_mesh.py @ 7461

Last change on this file since 7461 was 7461, checked in by steve, 11 years ago
File size: 38.4 KB
Line 
1#########################################################
2#
3#
4#  Read in a data file and subdivide the triangle list
5#
6#
7#  The final routine, pmesh_divide_metis, does automatic
8# grid partitioning. Once testing has finished on this
9# routine the others should be removed.
10#
11#  Authors: Linda Stals and Matthew Hardy, June 2005
12#  Modified: Linda Stals, Nov 2005
13#            Jack Kelly, Nov 2005
14#            Steve Roberts, Aug 2009 (updating to numpy)
15#
16#
17#########################################################
18
19
20import sys
21from os import sep
22from sys import path
23from math import floor
24
25import numpy as num
26
27from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
28
29#########################################################
30#
31# If the triangles list is reordered, the quantities
32# assigned to the triangles must also be reorded.
33#
34# *) quantities contain the quantites in the old ordering
35# *) proc_sum[i] contains the number of triangles in
36# processor i
37# *) tri_index is a map from the old triangle ordering to
38# the new ordering, where the new number for triangle
39# i is proc_sum[tri_index[i][0]]+tri_index[i][1]
40#
41# -------------------------------------------------------
42#
43# *) The quantaties are returned in the new ordering
44#
45#########################################################
46
47def reorder(quantities, tri_index, proc_sum):
48
49    # Find the number triangles
50
51    N = len(tri_index)
52
53    # Temporary storage area
54
55    index = num.zeros(N, num.int)
56    q_reord = {}
57
58    # Find the new ordering of the triangles
59
60    for i in range(N):
61        bin = tri_index[i][0]
62        bin_off_set = tri_index[i][1]
63        index[i] = proc_sum[bin]+bin_off_set
64
65    # Reorder each quantity according to the new ordering
66
67    for k in quantities:
68        q_reord[k] = num.zeros((N, 3), num.float)
69        for i in range(N):
70            q_reord[k][index[i]]=quantities[k].vertex_values[i]
71    del index
72
73    return q_reord
74
75
76#########################################################
77#
78# Divide the mesh using a call to metis, through pymetis.
79#
80# -------------------------------------------------------
81#
82# *)  The nodes, triangles, boundary, and quantities are
83# returned. triangles_per_proc defines the subdivision.
84# The first triangles_per_proc[0] triangles are assigned
85# to processor 0, the next triangles_per_proc[1] are
86# assigned to processor 1 etc. The boundary and quantites
87# are ordered the same way as the triangles
88#
89#########################################################
90
91#path.append('..' + sep + 'pymetis')
92
93try:
94    from pymetis.metis import partMeshNodal
95except ImportError:
96    print "***************************************************"
97    print "         Metis is probably not compiled."
98    print "         Read \anuga_core\source\pymetis\README"
99    print "***************************************************"
100    raise ImportError
101
102def pmesh_divide_metis(domain, n_procs):
103   
104    # Initialise the lists
105    # List, indexed by processor of # triangles.
106   
107    triangles_per_proc = []
108   
109    # List of lists, indexed by processor of vertex numbers
110   
111    tri_list = []
112   
113    # List indexed by processor of cumulative total of triangles allocated.
114   
115    proc_sum = []
116    for i in range(n_procs):
117        tri_list.append([])
118        triangles_per_proc.append(0)
119        proc_sum.append([])
120
121    # Prepare variables for the metis call
122   
123    n_tri = len(domain.triangles)
124    if n_procs != 1: #Because metis chokes on it...
125        n_vert = domain.get_number_of_nodes()
126        t_list = domain.triangles.copy()
127        t_list = num.reshape(t_list, (-1,))
128   
129        # The 1 here is for triangular mesh elements.
130        edgecut, epart, npart = partMeshNodal(n_tri, n_vert, t_list, 1, n_procs)
131        # print edgecut
132        # print npart
133        # print epart
134        del edgecut
135        del npart
136
137        # Sometimes (usu. on x86_64), partMeshNodal returnes an array of zero
138        # dimensional arrays. Correct this.
139        if type(epart[0]) == num.ndarray:
140            epart_new = num.zeros(len(epart), num.int)
141            for i in range(len(epart)):
142                epart_new[i] = epart[i][0]
143            epart = epart_new
144            del epart_new
145        # Assign triangles to processes, according to what metis told us.
146       
147        # tri_index maps triangle number -> processor, new triangle number
148        # (local to the processor)
149       
150        tri_index = {}
151        triangles = []       
152        for i in range(n_tri):
153            triangles_per_proc[epart[i]] = triangles_per_proc[epart[i]] + 1
154            tri_list[epart[i]].append(domain.triangles[i])
155            tri_index[i] = ([epart[i], len(tri_list[epart[i]]) - 1])
156       
157        # Order the triangle list so that all of the triangles belonging
158        # to processor i are listed before those belonging to processor
159        # i+1
160
161        for i in range(n_procs):
162            for t in tri_list[i]:
163                triangles.append(t)
164           
165        # The boundary labels have to changed in accoradance with the
166        # new triangle ordering, proc_sum and tri_index help with this
167
168        proc_sum[0] = 0
169        for i in range(n_procs - 1):
170            proc_sum[i+1]=proc_sum[i]+triangles_per_proc[i]
171
172        # Relabel the boundary elements to fit in with the new triangle
173        # ordering
174
175        boundary = {}
176        for b in domain.boundary:
177            t =  tri_index[b[0]]
178            boundary[proc_sum[t[0]]+t[1], b[1]] = domain.boundary[b]
179
180        quantities = reorder(domain.quantities, tri_index, proc_sum)
181    else:
182        boundary = domain.boundary.copy()
183        triangles_per_proc[0] = n_tri
184        triangles = domain.triangles.copy()
185       
186        # This is essentially the same as a chunk of code from reorder.
187       
188        quantities = {}
189        for k in domain.quantities:
190            quantities[k] = num.zeros((n_tri, 3), num.float)
191            for i in range(n_tri):
192                quantities[k][i] = domain.quantities[k].vertex_values[i]
193       
194    # Extract the node list
195   
196    nodes = domain.get_nodes().copy()
197   
198    # Convert the triangle datastructure to be an array type,
199    # this helps with the communication
200
201    ttriangles = num.zeros((len(triangles), 3), num.int)
202    for i in range(len(triangles)):
203        ttriangles[i] = triangles[i]
204   
205    return nodes, ttriangles, boundary, triangles_per_proc, quantities
206
207
208#########################################################
209#
210# Subdivide the domain. This module is primarily
211# responsible for building the ghost layer and
212# communication pattern
213#
214#
215#  Author: Linda Stals, June 2005
216#  Modified: Linda Stals, Nov 2005 (optimise python code)
217#            Steve Roberts, Aug 2009 (convert to numpy)
218#
219#
220#########################################################
221
222
223
224#########################################################
225#
226# Subdivide the triangles into non-overlapping domains.
227#
228#  *)  The subdivision is controlled by triangles_per_proc.
229# The first triangles_per_proc[0] triangles are assigned
230# to the first processor, the second triangles_per_proc[1]
231# are assigned to the second processor etc.
232#
233#  *) nodes, triangles and boundary contains all of the
234# nodes, triangles and boundary tag information for the
235# whole domain. The triangles should be orientated in the
236# correct way and the nodes number consecutively from 0.
237#
238# -------------------------------------------------------
239#
240#  *) A dictionary containing the full_nodes, full_triangles
241# and full_boundary information for each processor is
242# returned. The node information consists of
243# [global_id, x_coord, y_coord].
244#
245#########################################################
246
247def submesh_full(nodes, triangles, boundary, triangles_per_proc):
248
249    # Initialise
250
251    tlower = 0
252    nproc = len(triangles_per_proc)
253    nnodes = len(nodes)
254    node_list = []
255    triangle_list = []
256    boundary_list = []
257    submesh = {}
258    node_range = num.reshape(num.arange(nnodes),(nnodes,1))
259
260    #print node_range
261    tsubnodes = num.concatenate((node_range, nodes), 1)
262
263
264    # Loop over processors
265
266    for p in range(nproc):
267
268        # Find triangles on processor p
269
270        tupper = triangles_per_proc[p]+tlower
271        subtriangles = triangles[tlower:tupper]
272        triangle_list.append(subtriangles)
273
274        # Find the boundary edges on processor p
275
276        subboundary = {}
277        for k in boundary:
278            if (k[0] >=tlower and k[0] < tupper):
279                subboundary[k]=boundary[k]
280        boundary_list.append(subboundary)
281
282        # Find nodes in processor p
283
284        nodemap = num.zeros(nnodes, 'i')
285        for t in subtriangles:
286            nodemap[t[0]]=1
287            nodemap[t[1]]=1
288            nodemap[t[2]]=1
289
290       
291        node_list.append(tsubnodes.take(num.flatnonzero(nodemap),axis=0))
292
293        # Move to the next processor
294
295        tlower = tupper
296
297    # Put the results in a dictionary
298
299    submesh["full_nodes"] = node_list
300    submesh["full_triangles"] = triangle_list
301    submesh["full_boundary"] = boundary_list
302
303    # Clean up before exiting
304
305    del (nodemap)
306
307    return submesh
308
309
310#########################################################
311#
312# Build the ghost layer of triangles
313#
314#  *) Given the triangle subpartion for the processor
315# build a ghost layer of triangles. The ghost layer
316# consists of two layers of neighbouring triangles.
317#
318#  *) The vertices in the ghost triangles must also
319# be added to the node list for the current processor
320#
321#
322# -------------------------------------------------------
323#
324#  *) The extra triangles and nodes are returned.
325#
326#  *)  The node information consists of
327# [global_id, x_coord, y_coord].
328#
329#  *) The triangle information consists of
330# [triangle number, t], where t = [v1, v2, v3].
331#
332#########################################################
333
334def ghost_layer(submesh, mesh, p, tupper, tlower):
335
336    ncoord = mesh.number_of_nodes
337    ntriangles = mesh.number_of_triangles
338
339    # Find the first layer of boundary triangles
340
341    trianglemap = num.zeros(ntriangles, 'i')
342    for t in range(tlower, tupper):
343       
344        n = mesh.neighbours[t, 0]
345
346        if n >= 0:
347            if n < tlower or n >= tupper:
348                trianglemap[n] = 1
349        n = mesh.neighbours[t, 1]
350        if n >= 0:
351            if n < tlower or n >= tupper:
352                trianglemap[n] = 1
353        n = mesh.neighbours[t, 2]
354        if n >= 0:
355            if n < tlower or n >= tupper:
356                trianglemap[n] = 1
357
358    # Find the second layer of boundary triangles
359
360    for t in range(len(trianglemap)):
361        if trianglemap[t]==1:
362            n = mesh.neighbours[t, 0]
363            if n >= 0:
364                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
365                    trianglemap[n] = 2
366            n = mesh.neighbours[t, 1]
367            if n >= 0:
368                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
369                    trianglemap[n] = 2
370            n = mesh.neighbours[t, 2]
371            if n >= 0:
372                if (n < tlower or n >= tupper) and trianglemap[n] == 0:
373                    trianglemap[n] = 2
374
375    # Build the triangle list and make note of the vertices
376
377    nodemap = num.zeros(ncoord, 'i')
378    fullnodes = submesh["full_nodes"][p]
379
380    subtriangles = []
381    for i in range(len(trianglemap)):
382        if trianglemap[i] != 0:
383            t = list(mesh.triangles[i])
384            nodemap[t[0]] = 1
385            nodemap[t[1]] = 1
386            nodemap[t[2]] = 1
387
388    trilist = num.reshape(num.arange(ntriangles),(ntriangles,1))
389    tsubtriangles = num.concatenate((trilist, mesh.triangles), 1)
390    subtriangles = tsubtriangles.take(num.flatnonzero(trianglemap),axis=0)
391
392   
393    # Keep a record of the triangle vertices, if they are not already there
394
395    subnodes = []
396    for n in fullnodes:
397        nodemap[int(n[0])] = 0
398
399    nodelist = num.reshape(num.arange(ncoord),(ncoord,1))
400    tsubnodes = num.concatenate((nodelist, mesh.get_nodes()), 1)
401    subnodes = tsubnodes.take(num.flatnonzero(nodemap),axis=0)
402
403    # Clean up before exiting
404
405    del (nodelist)
406    del (trilist)
407    del (tsubnodes)
408    del (nodemap)
409    del (trianglemap)
410
411    # Return the triangles and vertices sitting on the boundary layer
412
413    return subnodes, subtriangles
414
415#########################################################
416#
417# Find the edges of the ghost trianlges that do not
418# have a neighbour in the current cell. These are
419# treated as a special type of boundary edge.
420#
421#  *) Given the ghost triangles in a particular
422# triangle, use the mesh to find its neigbours. If
423# the neighbour is not in the processor set it to
424# be a boundary edge
425#
426#  *) The vertices in the ghost triangles must also
427# be added to the node list for the current processor
428#
429#  *) The boundary edges for the ghost triangles are
430# ignored.
431#
432# -------------------------------------------------------
433#
434#  *) The type assigned to the ghost boundary edges is 'ghost'
435#
436#  *)  The boundary information is returned as a directorier
437# with the key = (triangle id, edge no) and the values
438# assigned to the key is 'ghost'
439#
440#
441#########################################################
442def is_in_processor(ghost_list, tlower, tupper, n):
443
444    return num.equal(ghost_list,n).any() or (tlower <= n and tupper > n)
445
446
447def ghost_bnd_layer(ghosttri, tlower, tupper, mesh, p):
448
449    ghost_list = []
450    subboundary = {}
451
452
453    for t in ghosttri:
454        ghost_list.append(t[0])
455   
456    for t in ghosttri:
457
458        n = mesh.neighbours[t[0], 0]
459        if not is_in_processor(ghost_list, tlower, tupper, n):
460            subboundary[t[0], 0] = 'ghost'
461
462        n = mesh.neighbours[t[0], 1]
463        if not is_in_processor(ghost_list, tlower, tupper, n):
464            subboundary[t[0], 1] = 'ghost'
465
466        n = mesh.neighbours[t[0], 2]
467        if not is_in_processor(ghost_list, tlower, tupper, n):
468            subboundary[t[0], 2] = 'ghost'
469           
470    return subboundary
471
472#########################################################
473#
474# The ghost triangles on the current processor will need
475# to get updated information from the neighbouring
476# processor containing the corresponding full triangles.
477#
478#  *) The tri_per_proc is used to determine which
479# processor contains the full node copy.
480#
481# -------------------------------------------------------
482#
483#  *) The ghost communication pattern consists of
484# [global node number, neighbour processor number].
485#
486#########################################################
487
488def ghost_commun_pattern(subtri, p, tri_per_proc):
489
490    # Loop over the ghost triangles
491
492    ghost_commun = num.zeros((len(subtri), 2), num.int)
493
494    for i in range(len(subtri)):
495        global_no = subtri[i][0]
496
497        # Find which processor contains the full triangle
498
499        nproc = len(tri_per_proc)
500        neigh = nproc-1
501        sum = 0
502        for q in range(nproc-1):
503            if (global_no < sum+tri_per_proc[q]):
504                neigh = q
505                break
506            sum = sum+tri_per_proc[q]
507
508        # Keep a copy of the neighbour processor number
509
510        ghost_commun[i] = [global_no, neigh]
511
512    return ghost_commun
513
514#########################################################
515#
516# The full triangles in this processor must communicate
517# updated information to neighbouring processor that
518# contain ghost triangles
519#
520#  *) The ghost communication pattern for all of the
521# processor must be built before calling this processor.
522#
523#  *) The full communication pattern is found by looping
524# through the ghost communication pattern for all of the
525# processors. Recall that this information is stored in
526# the form [global node number, neighbour processor number].
527# The full communication for the neighbour processor is
528# then updated.
529#
530# -------------------------------------------------------
531#
532#  *) The full communication pattern consists of
533# [global id, [p1, p2, ...]], where p1, p2 etc contain
534# a ghost node copy of the triangle global id.
535#
536#########################################################
537
538def full_commun_pattern(submesh, tri_per_proc):
539    tlower = 0
540    nproc = len(tri_per_proc)
541    full_commun = []
542
543    # Loop over the processor
544
545    for p in range(nproc):
546
547        # Loop over the full triangles in the current processor
548        # and build an empty dictionary
549
550        fcommun = {}
551        tupper = tri_per_proc[p]+tlower
552        for i in range(tlower, tupper):
553            fcommun[i] = []
554        full_commun.append(fcommun)
555        tlower = tupper
556
557    # Loop over the processor again
558
559    for p in range(nproc):
560
561        # Loop over the ghost triangles in the current processor,
562        # find which processor contains the corresponding full copy
563        # and note that the processor must send updates to this
564        # processor
565
566        for g in submesh["ghost_commun"][p]:
567            neigh = g[1]
568            full_commun[neigh][g[0]].append(p)
569
570    return full_commun
571
572
573#########################################################
574#
575# Given the non-overlapping grid partition, an extra layer
576# of triangles are included to help with the computations.
577# The triangles in this extra layer are not updated by
578# the processor, their updated values must be sent by the
579# processor containing the original, full, copy of the
580# triangle. The communication pattern that controls these
581# updates must also be built.
582#
583#  *) Assumes that full triangles, nodes etc have already
584# been found and stored in submesh
585#
586#  *) See the documentation for ghost_layer,
587# ghost_commun_pattern and full_commun_pattern
588#
589# -------------------------------------------------------
590#
591#  *) The additional information is added to the submesh
592# dictionary. See the documentation for ghost_layer,
593# ghost_commun_pattern and full_commun_pattern
594#
595#  *) The ghost_triangles, ghost_nodes, ghost_boundary,
596# ghost_commun and full_commun is added to submesh
597#########################################################
598
599def submesh_ghost(submesh, mesh, triangles_per_proc):
600
601    nproc = len(triangles_per_proc)
602    tlower = 0
603    ghost_triangles = []
604    ghost_nodes = []
605    ghost_commun = []
606    ghost_bnd = []
607
608    # Loop over the processors
609
610    for p in range(nproc):
611
612        # Find the full triangles in this processor
613
614        tupper = triangles_per_proc[p]+tlower
615
616        # Build the ghost boundary layer
617
618        [subnodes, subtri] = \
619                   ghost_layer(submesh, mesh, p, tupper, tlower)
620        ghost_triangles.append(subtri)
621        ghost_nodes.append(subnodes)
622
623
624        # Find the boundary layer formed by the ghost triangles
625       
626        subbnd = ghost_bnd_layer(subtri, tlower, tupper, mesh, p)
627        ghost_bnd.append(subbnd)
628       
629        # Build the communication pattern for the ghost nodes
630
631        gcommun = \
632                ghost_commun_pattern(subtri, p, triangles_per_proc)
633        ghost_commun.append(gcommun)
634
635        # Move to the next processor
636
637        tlower = tupper
638
639
640    # Record the ghost layer and communication pattern
641
642    submesh["ghost_nodes"] = ghost_nodes
643    submesh["ghost_triangles"] = ghost_triangles
644    submesh["ghost_commun"] = ghost_commun
645    submesh["ghost_boundary"] = ghost_bnd
646   
647    # Build the communication pattern for the full triangles
648
649    full_commun = full_commun_pattern(submesh, triangles_per_proc)
650    submesh["full_commun"] = full_commun
651
652    # Return the submesh
653
654    return submesh
655
656
657#########################################################
658#
659# Certain quantities may be assigned to the triangles,
660# these quantities must be subdivided in the same way
661# as the triangles
662#
663#  *) The quantities are ordered in the same way as the
664# triangles
665#
666# -------------------------------------------------------
667#
668#  *) The quantites attached to the full triangles are
669# stored in full_quan
670#
671#  *) The quantities attached to the ghost triangles are
672# stored in ghost_quan
673#########################################################
674
675def submesh_quantities(submesh, quantities, triangles_per_proc):
676
677    nproc = len(triangles_per_proc)
678
679    lower = 0
680
681    # Build an empty dictionary to hold the quantites
682
683    submesh["full_quan"] = {}
684    submesh["ghost_quan"] = {}
685    for k in quantities:
686        submesh["full_quan"][k] = []
687        submesh["ghost_quan"][k] = []
688
689    # Loop trough the subdomains
690
691    for p in range(nproc):
692        upper =   lower+triangles_per_proc[p]
693
694        # Find the global ID of the ghost triangles
695
696        global_id = []
697        M = len(submesh["ghost_triangles"][p])
698        for j in range(M):
699            global_id.append(submesh["ghost_triangles"][p][j][0])
700
701        # Use the global ID to extract the quantites information from
702        # the full domain
703
704        for k in quantities:
705            submesh["full_quan"][k].append(quantities[k][lower:upper])
706            submesh["ghost_quan"][k].append(num.zeros( (M,3) , num.float))
707            for j in range(M):
708                submesh["ghost_quan"][k][p][j] = \
709                                               quantities[k][global_id[j]]
710
711        lower = upper
712
713    return submesh
714
715#########################################################
716#
717# Build the grid partition on the host.
718#
719#  *) See the documentation for submesh_ghost and
720# submesh_full
721#
722# -------------------------------------------------------
723#
724#  *) A dictionary containing the full_triangles,
725# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
726# ghost_boundary, ghost_commun and full_commun and true boundary polygon is returned.
727#
728#########################################################
729
730def build_submesh(nodes, triangles, edges, quantities,
731                  triangles_per_proc):
732
733    # Temporarily build the mesh to find the neighbouring
734    # triangles and true boundary polygon
735
736    mesh = Mesh(nodes, triangles)
737    boundary_polygon = mesh.get_boundary_polygon()
738   
739
740    # Subdivide into non-overlapping partitions
741
742    submeshf = submesh_full(nodes, triangles, edges, \
743                            triangles_per_proc)
744   
745    # Add any extra ghost boundary layer information
746
747    submeshg = submesh_ghost(submeshf, mesh, triangles_per_proc)
748
749    # Order the quantities information to be the same as the triangle
750    # information
751
752    submesh = submesh_quantities(submeshg, quantities, \
753                                 triangles_per_proc)
754
755    submesh["boundary_polygon"] = boundary_polygon
756    return submesh
757
758#########################################################
759#
760#  Given the subdivision of the grid assigned to the
761# current processor convert it into a form that is
762# appropriate for the GA datastructure.
763#
764#  The main function of these modules is to change the
765# node numbering. The GA datastructure assumes they
766# are numbered consecutively from 0.
767#
768#  The module also changes the communication pattern
769# datastructure into a form needed by parallel_advection
770#
771#  Authors: Linda Stals and Matthew Hardy, June 2005
772#  Modified: Linda Stals, Nov 2005 (optimise python code)
773#            Steve Roberts, Aug 2009 (updating to numpy)
774#
775#
776#########################################################
777
778 
779#########################################################
780# Convert the format of the data to that used by ANUGA
781#
782#
783# *) Change the nodes global ID's to an integer value,
784#starting from 0.
785#
786# *) The triangles and boundary edges must also be
787# updated accordingly.
788#
789# -------------------------------------------------------
790#
791# *) The nodes, triangles and boundary edges defined by
792# the new numbering scheme are returned
793#
794#########################################################
795
796def build_local_GA(nodes, triangles, boundaries, tri_index):
797
798    Nnodes =len(nodes)
799    Ntriangles = len(triangles)
800   
801    # Extract the nodes (using the local ID)
802   
803    GAnodes = num.take(nodes, (1, 2), 1)
804
805    # Build a global ID to local ID mapping
806
807    NGlobal = 0
808    for i in range(Nnodes):
809        if nodes[i][0] > NGlobal:
810            NGlobal = nodes[i][0]
811    index = num.zeros(int(NGlobal)+1, num.int)
812    num.put(index, num.take(nodes, (0,), 1).astype(num.int), \
813        num.arange(Nnodes))
814       
815    # Change the global IDs in the triangles to the local IDs
816
817    GAtriangles = num.zeros((Ntriangles, 3), num.int)
818    GAtriangles[:,0] = num.take(index, triangles[:,0])
819    GAtriangles[:,1] = num.take(index, triangles[:,1])
820    GAtriangles[:,2] = num.take(index, triangles[:,2])
821
822    # Change the triangle numbering in the boundaries
823
824    GAboundaries = {}
825    for b in boundaries:
826        GAboundaries[tri_index[b[0]], b[1]] = boundaries[b]
827       
828    del (index)
829   
830    return GAnodes, GAtriangles, GAboundaries
831
832
833#########################################################
834# Change the communication format to that needed by the
835# parallel advection file.
836#
837# *) The index contains [global triangle no,
838# local triangle no.]
839#
840# -------------------------------------------------------
841#
842# *) The ghost_recv and full_send dictionaries are
843# returned.
844#
845# *) ghost_recv dictionary is local id, global id, value
846#
847# *) full_recv dictionary is local id, global id, value
848#
849# *) The information is ordered by the global id. This
850# means that the communication order is predetermined and
851# local and global id do not need to be
852# compared when the information is sent/received.
853#
854#########################################################
855
856def build_local_commun(index, ghostc, fullc, nproc):
857
858    # Initialise
859
860    full_send = {}
861    ghost_recv = {}
862
863    # Build the ghost_recv dictionary (sort the
864    # information by the global numbering)
865   
866    ghostc = num.sort(ghostc, 0)
867   
868    for c in range(nproc):
869        s = ghostc[:,0]
870        d = num.compress(num.equal(ghostc[:,1],c), s)
871        if len(d) > 0:
872            ghost_recv[c] = [0, 0]
873            ghost_recv[c][1] = d
874            ghost_recv[c][0] = num.take(index, d)
875           
876    # Build a temporary copy of the full_send dictionary
877    # (this version allows the information to be stored
878    # by the global numbering)
879
880    tmp_send = {}
881    for global_id in fullc:
882        for i in range(len(fullc[global_id])):
883            neigh = fullc[global_id][i]
884            if not tmp_send.has_key(neigh):
885                tmp_send[neigh] = []
886            tmp_send[neigh].append([global_id, \
887                                    index[global_id]])
888
889    # Extract the full send information and put it in the form
890    # required for the full_send dictionary
891
892    for neigh in tmp_send:
893        neigh_commun = num.sort(tmp_send[neigh], 0)
894        full_send[neigh] = [0, 0]
895        full_send[neigh][0] = neigh_commun[:,1]
896        full_send[neigh][1] = neigh_commun[:,0]
897
898    return ghost_recv, full_send
899
900
901#########################################################
902# Convert the format of the data to that used by ANUGA
903#
904#
905# *) Change the nodes global ID's to an integer value,
906# starting from 0. The node numbering in the triangles
907# must also be updated to take this into account.
908#
909# *) The triangle number will also change, which affects
910# the boundary tag information and the communication
911# pattern.
912#
913# -------------------------------------------------------
914#
915# *) The nodes, triangles, boundary edges and communication
916# pattern defined by the new numbering scheme are returned
917#
918#########################################################
919
920def build_local_mesh(submesh, lower_t, upper_t, nproc):
921
922    # Combine the full nodes and ghost nodes
923
924    nodes = num.concatenate((submesh["full_nodes"], \
925                         submesh["ghost_nodes"]))
926   
927    # Combine the full triangles and ghost triangles
928
929    gtri =  num.take(submesh["ghost_triangles"],(1, 2, 3),1)
930    triangles = num.concatenate((submesh["full_triangles"], gtri))
931
932    # Combine the full boundaries and ghost boundaries
933
934    boundaries = submesh["full_boundary"]
935    for b in submesh["ghost_boundary"]:
936        boundaries[b]=submesh["ghost_boundary"][b]
937
938    # Make note of the new triangle numbers, including the ghost
939    # triangles
940
941    NGlobal = upper_t
942    for i in range(len(submesh["ghost_triangles"])):
943        id = submesh["ghost_triangles"][i][0]
944        if id > NGlobal:
945            NGlobal = id
946    index = num.zeros(int(NGlobal)+1, num.int)
947    index[lower_t:upper_t]=num.arange(upper_t-lower_t)
948    for i in range(len(submesh["ghost_triangles"])):
949        index[submesh["ghost_triangles"][i][0]] = i+upper_t-lower_t
950   
951    # Change the node numbering (and update the numbering in the
952    # triangles)
953
954    [GAnodes, GAtriangles, GAboundary] = build_local_GA(nodes, triangles, boundaries, index)
955
956    # Extract the local quantities
957   
958    quantities ={}
959    for k in submesh["full_quan"]:
960        Nf = len(submesh["full_quan"][k])
961        Ng = len(submesh["ghost_quan"][k])
962        quantities[k] = num.zeros((Nf+Ng, 3), num.float)
963        quantities[k][0:Nf] = submesh["full_quan"][k] 
964        quantities[k][Nf:Nf+Ng] = submesh["ghost_quan"][k]
965                             
966    # Change the communication pattern into a form needed by
967    # the parallel_adv
968
969    gcommun = submesh["ghost_commun"]
970    fcommun = submesh["full_commun"]
971    [ghost_rec, full_send] = \
972                build_local_commun(index, gcommun, fcommun, nproc)
973
974    # Clean up before exiting
975
976    del(index)
977
978    return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \
979           full_send
980
981
982#########################################################
983#
984# Handle the communication between the host machine
985# (processor 0) and the processors. The host machine is
986# responsible for the doing the initial grid partitioning.
987#
988# The routines given below should be moved to the
989# build_submesh.py and build_local.py file to allow
990# overlapping of  communication and computation.
991# This should be done after more debugging.
992#
993#
994#  Author: Linda Stals, June 2005
995#  Modified: Linda Stals, Nov 2005 (optimise python code)
996#            Steve Roberts, Aug 2009 (update to numpy)
997#
998#
999#########################################################
1000
1001
1002#########################################################
1003#
1004# Send the submesh to processor p.
1005#
1006# *) The order and form is strongly coupled with
1007# rec_submesh.
1008#
1009# -------------------------------------------------------
1010#
1011# *) All of the information has been sent to processor p.
1012#
1013#########################################################
1014
1015def send_submesh(submesh, triangles_per_proc, p, verbose=True):
1016
1017    import pypar
1018   
1019    myid = pypar.rank()
1020   
1021    if verbose: print 'process %d sending submesh to process %d' %(myid, p)
1022   
1023    # build and send the tagmap for the boundary conditions
1024   
1025    tagmap = {}
1026    counter = 1
1027    for b in submesh["full_boundary"][p]:
1028         bkey = submesh["full_boundary"][p][b]
1029         if not tagmap.has_key(bkey):
1030             tagmap[bkey] = counter
1031             counter = counter+1
1032    for b in submesh["ghost_boundary"][p]:
1033         bkey = submesh["ghost_boundary"][p][b]
1034         if not tagmap.has_key(bkey):
1035             tagmap[bkey] = counter
1036             counter = counter+1
1037
1038    pypar.send(tagmap, p)
1039
1040    # send the quantities key information
1041
1042    pypar.send(submesh["full_quan"].keys(), p)
1043
1044    # send the number of triangles per processor
1045
1046    pypar.send(triangles_per_proc, p)
1047
1048    # compress full_commun
1049
1050    flat_full_commun = []
1051
1052    for c in submesh["full_commun"][p]:
1053        for i in range(len(submesh["full_commun"][p][c])):
1054            flat_full_commun.append([c,submesh["full_commun"][p][c][i]])
1055
1056    # send the array sizes so memory can be allocated
1057
1058    setup_array = num.zeros((9,),num.int)
1059    setup_array[0] = len(submesh["full_nodes"][p])
1060    setup_array[1] = len(submesh["ghost_nodes"][p])
1061    setup_array[2] = len(submesh["full_triangles"][p])
1062    setup_array[3] = len(submesh["ghost_triangles"][p])
1063    setup_array[4] = len(submesh["full_boundary"][p])
1064    setup_array[5] = len(submesh["ghost_boundary"][p])
1065    setup_array[6] = len(submesh["ghost_commun"][p])
1066    setup_array[7] = len(flat_full_commun)
1067    setup_array[8] = len(submesh["full_quan"])
1068   
1069    pypar.send(num.array(setup_array, num.int), p)
1070   
1071    # send the nodes
1072
1073    pypar.send(num.array(submesh["full_nodes"][p], num.float), p)
1074    pypar.send(num.array(submesh["ghost_nodes"][p], num.float),p)
1075
1076    # send the triangles
1077
1078    pypar.send(num.array(submesh["full_triangles"][p],  num.int), p)
1079    pypar.send(num.array(submesh["ghost_triangles"][p], num.int), p)
1080
1081    # send the boundary
1082
1083    bc = []
1084    for b in submesh["full_boundary"][p]:
1085        bc.append([b[0], b[1], tagmap[submesh["full_boundary"][p][b]]])
1086
1087
1088    pypar.send(num.array(bc, num.int), p)
1089
1090    bc = []
1091    for b in submesh["ghost_boundary"][p]:
1092        bc.append([b[0], b[1], tagmap[submesh["ghost_boundary"][p][b]]])
1093
1094    pypar.send(num.array(bc, num.int), p)
1095
1096    # send the communication pattern
1097
1098    pypar.send(submesh["ghost_commun"][p], p)
1099
1100    pypar.send(num.array(flat_full_commun, num.int), p)
1101
1102    # send the quantities
1103   
1104    for k in submesh["full_quan"]:
1105        pypar.send(num.array(submesh["full_quan"][k][p], num.float), p)
1106       
1107    for k in submesh["ghost_quan"]:
1108        pypar.send(num.array(submesh["ghost_quan"][k][p], num.float),p)
1109       
1110
1111#########################################################
1112#
1113# Receive the submesh from processor p.
1114#
1115# *) The order and form is strongly coupled with
1116# send_submesh.
1117#
1118# -------------------------------------------------------
1119#
1120# *) All of the information has been received by the
1121# processor p and passed into build_local.
1122#
1123# *) The information is returned in a form needed by the
1124# GA datastructure.
1125#
1126#########################################################
1127
1128def rec_submesh_flat(p, verbose=True):
1129
1130    import pypar
1131   
1132    numproc = pypar.size()
1133    myid = pypar.rank()
1134
1135    submesh_cell = {}
1136   
1137    if verbose: print 'process %d receiving submesh from process %d' %(myid, p)
1138
1139    # receive the tagmap for the boundary conditions
1140   
1141    tagmap = pypar.receive(p)
1142
1143    itagmap = {}
1144    for t in tagmap:
1145        itagmap[tagmap[t]]=t
1146
1147    # receive the quantities key information
1148
1149    qkeys = pypar.receive(p)
1150
1151    # receive the number of triangles per processor
1152
1153    triangles_per_proc = pypar.receive(p)
1154
1155    # recieve information about the array sizes
1156
1157    setup_array = pypar.receive(p)
1158
1159    no_full_nodes      = setup_array[0]
1160    no_ghost_nodes     = setup_array[1]
1161    no_full_triangles  = setup_array[2]
1162    no_ghost_triangles = setup_array[3]
1163    no_full_boundary   = setup_array[4]
1164    no_ghost_boundary  = setup_array[5]
1165    no_ghost_commun    = setup_array[6]
1166    no_full_commun     = setup_array[7]
1167    no_quantities      = setup_array[8]
1168   
1169    # receive the full nodes
1170
1171    submesh_cell["full_nodes"] = pypar.receive(p)
1172
1173    # receive the ghost nodes
1174
1175    submesh_cell["ghost_nodes"] = pypar.receive(p)
1176   
1177    # receive the full triangles
1178
1179    submesh_cell["full_triangles"] = pypar.receive(p)
1180   
1181    # receive the ghost triangles
1182
1183    submesh_cell["ghost_triangles"] = pypar.receive(p)
1184
1185    # receive the full boundary
1186
1187    bnd_c = pypar.receive(p)
1188
1189    submesh_cell["full_boundary"] = {}
1190    for b in bnd_c:
1191        submesh_cell["full_boundary"][b[0],b[1]]=itagmap[b[2]]
1192
1193    # receive the ghost boundary
1194
1195    bnd_c = pypar.receive(p)
1196
1197    submesh_cell["ghost_boundary"] = {}
1198    for b in bnd_c:
1199        submesh_cell["ghost_boundary"][b[0],b[1]]=itagmap[b[2]]
1200
1201    # receive the ghost communication pattern
1202
1203    submesh_cell["ghost_commun"] = pypar.receive(p)
1204   
1205    # receive the full communication pattern
1206
1207    full_commun = pypar.receive(p)
1208
1209    submesh_cell["full_commun"] = {}
1210    for c in full_commun:
1211        submesh_cell["full_commun"][c[0]] = []
1212    for c in full_commun:
1213        submesh_cell["full_commun"][c[0]].append(c[1])
1214
1215    # receive the quantities
1216
1217    submesh_cell["full_quan"]={}
1218   
1219    for i in range(no_quantities):
1220        tmp = pypar.receive(p)
1221        submesh_cell["full_quan"][qkeys[i]]=num.zeros((no_full_triangles,3), num.float)
1222        submesh_cell["full_quan"][qkeys[i]][:] = tmp[:]
1223
1224    submesh_cell["ghost_quan"]={}
1225    for i in range(no_quantities):
1226        tmp = pypar.receive(p)
1227        submesh_cell["ghost_quan"][qkeys[i]]= num.zeros((no_ghost_triangles,3), num.float)
1228        submesh_cell["ghost_quan"][qkeys[i]][:] = tmp[:]
1229   
1230    return submesh_cell, triangles_per_proc,\
1231           no_full_nodes, no_full_triangles
1232
1233
1234
1235#########################################################
1236#
1237# Receive the submesh from processor p.
1238#
1239# *) The order and form is strongly coupled with
1240# send_submesh.
1241#
1242# -------------------------------------------------------
1243#
1244# *) All of the information has been received by the
1245# processor p and passed into build_local.
1246#
1247# *) The information is returned in a form needed by the
1248# GA datastructure.
1249#
1250#########################################################
1251
1252def rec_submesh(p, verbose=True):
1253
1254    import pypar
1255   
1256    numproc = pypar.size()
1257    myid = pypar.rank()
1258
1259    [submesh_cell, triangles_per_proc,\
1260     number_of_full_nodes, number_of_full_triangles] = rec_submesh_flat(p,verbose)
1261   
1262    # find the full triangles assigned to this processor
1263
1264    lower_t = 0
1265    for i in range(myid):
1266        lower_t = lower_t+triangles_per_proc[i]
1267    upper_t = lower_t+triangles_per_proc[myid]
1268
1269    # convert the information into a form needed by the GA
1270    # datastructure
1271
1272    [GAnodes, GAtriangles, boundary, quantities, ghost_rec, full_send] = \
1273              build_local_mesh(submesh_cell, lower_t, upper_t, \
1274                               numproc)
1275   
1276    return GAnodes, GAtriangles, boundary, quantities,\
1277           ghost_rec, full_send,\
1278           number_of_full_nodes, number_of_full_triangles
1279
1280
1281#########################################################
1282#
1283# Extract the submesh that will belong to the
1284# "host processor" (i.e. processor zero)
1285#
1286#  *) See the documentation for build_submesh
1287#
1288# -------------------------------------------------------
1289#
1290#  *) A dictionary containing the full_triangles,
1291# full_nodes, full_boundary, ghost_triangles, ghost_nodes,
1292# ghost_boundary, ghost_commun and full_commun belonging
1293# to processor zero are returned.
1294#
1295#########################################################
1296def extract_hostmesh(submesh, triangles_per_proc):
1297
1298   
1299    submesh_cell = {}
1300    submesh_cell["full_nodes"] = submesh["full_nodes"][0]
1301    submesh_cell["ghost_nodes"] = submesh["ghost_nodes"][0]
1302    submesh_cell["full_triangles"] = submesh["full_triangles"][0]
1303    submesh_cell["ghost_triangles"] = submesh["ghost_triangles"][0]
1304    submesh_cell["full_boundary"] = submesh["full_boundary"][0]
1305    submesh_cell["ghost_boundary"] = submesh["ghost_boundary"][0]
1306    submesh_cell["ghost_commun"] = submesh["ghost_commun"][0]
1307    submesh_cell["full_commun"] = submesh["full_commun"][0]
1308    submesh_cell["full_quan"] ={}
1309    submesh_cell["ghost_quan"]={}
1310    for k in submesh["full_quan"]:
1311        submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
1312        submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
1313
1314    numprocs = len(triangles_per_proc)
1315    points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict = \
1316            build_local_mesh(submesh_cell, 0, triangles_per_proc[0], numprocs)
1317
1318
1319    return  points, vertices, boundary, quantities, ghost_recv_dict, \
1320           full_send_dict
1321
1322
1323
1324
Note: See TracBrowser for help on using the repository browser.