source: anuga_core/source/anuga_parallel/parallel_meshes.py @ 7447

Last change on this file since 7447 was 7447, checked in by steve, 14 years ago

Concentrating code

File size: 10.3 KB
Line 
1"""parallel-meshes -
22D triangular domains for parallel finite-volume computations of
3the advection equation, with extra structures to define the
4sending and receiving communications define in dictionaries
5full_send_dict and ghost_recv_dict
6
7
8Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
9Geoscience Australia, 2005
10
11Modified by Linda Stals, March 2006, to include ghost boundaries
12
13"""
14
15
16import sys
17#from Numeric import array, zeros, Float, Int, ones, sum
18
19import numpy as num
20import pypar
21
22from anuga.config import epsilon
23
24
25
26
27
28def parallel_rectangle(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
29
30
31    """Setup a rectangular grid of triangles
32    with m+1 by n+1 grid points
33    and side lengths len1, len2. If side lengths are omitted
34    the mesh defaults to the unit square, divided between all the
35    processors
36
37    len1: x direction (left to right)
38    len2: y direction (bottom to top)
39
40    """
41
42    processor = pypar.rank()
43    numproc   = pypar.size()
44
45    #print 'numproc',numproc
46    #print 'processor ',processor
47
48    m_low, m_high = pypar.balance(m_g, numproc, processor)
49   
50    n = n_g
51    m_low  = m_low-1
52    m_high = m_high+1
53
54    #print 'm_low, m_high', m_low, m_high
55    m = m_high - m_low
56
57    delta1 = float(len1_g)/m_g
58    delta2 = float(len2_g)/n_g
59
60    len1 = len1_g*float(m)/float(m_g)
61    len2 = len2_g
62    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
63
64    #Calculate number of points
65    Np = (m+1)*(n+1)
66
67    class VIndex:
68
69        def __init__(self, n,m):
70            self.n = n
71            self.m = m
72
73        def __call__(self, i,j):
74            return j+i*(self.n+1)
75
76    class EIndex:
77
78        def __init__(self, n,m):
79            self.n = n
80            self.m = m
81
82        def __call__(self, i,j):
83            return 2*(j+i*self.n)
84
85
86    I = VIndex(n,m)
87    E = EIndex(n,m)
88
89    points = num.zeros( (Np,2), num.float)
90
91    for i in range(m+1):
92        for j in range(n+1):
93
94            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
95
96    #Construct 2 triangles per rectangular element and assign tags to boundary
97    #Calculate number of triangles
98    Nt = 2*m*n
99
100
101    elements = num.zeros( (Nt,3), num.int)
102    boundary = {}
103    Idgl = []
104    Idfl = []
105    Idgr = []
106    Idfr = []
107
108    full_send_dict = {}
109    ghost_recv_dict = {}
110    nt = -1
111    for i in range(m):
112        for j in range(n):
113
114            i1 = I(i,j+1)
115            i2 = I(i,j)
116            i3 = I(i+1,j+1)
117            i4 = I(i+1,j)
118
119            #Lower Element
120            nt = E(i,j)
121            if i == 0:
122                Idgl.append(nt)
123
124            if i == 1:
125                Idfl.append(nt)
126
127            if i == m-2:
128                Idfr.append(nt)
129
130            if i == m-1:
131                Idgr.append(nt)
132
133            if i == m-1:
134                if processor == numproc-1:
135                    boundary[nt, 2] = 'right'
136                else:
137                    boundary[nt, 2] = 'ghost'
138       
139            if j == 0:
140                boundary[nt, 1] = 'bottom'
141            elements[nt,:] = [i4,i3,i2]
142
143            #Upper Element
144            nt = E(i,j)+1
145            if i == 0:
146                Idgl.append(nt)
147
148            if i == 1:
149                Idfl.append(nt)
150
151            if i == m-2:
152                Idfr.append(nt)
153
154            if i == m-1:
155                Idgr.append(nt)
156
157            if i == 0:
158                if processor == 0:
159                    boundary[nt, 2] = 'left'
160                else:
161                    boundary[nt, 2] = 'ghost'
162            if j == n-1:
163                boundary[nt, 1] = 'top'
164            elements[nt,:] = [i1,i2,i3]
165
166    if numproc==1:
167        Idfl.extend(Idfr)
168        Idgr.extend(Idgl)
169
170        #print Idfl
171        #print Idgr
172       
173        Idfl = num.array(Idfl,num.int)
174        Idgr = num.array(Idgr,num.int)
175
176        #print Idfl
177        #print Idgr
178       
179        full_send_dict[processor]  = [Idfl, Idfl]
180        ghost_recv_dict[processor] = [Idgr, Idgr]
181
182        #print  full_send_dict[processor]
183        #print ghost_recv_dict[processor]
184    elif numproc == 2:
185        Idfl.extend(Idfr)
186        Idgr.extend(Idgl)
187        Idfl = num.array(Idfl,num.int)
188        Idgr = num.array(Idgr,num.int)
189        full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl]
190        ghost_recv_dict[(processor-1)%numproc] = [Idgr, Idgr]
191    else:
192        Idfl = num.array(Idfl,num.int)
193        Idgl = num.array(Idgl,num.int)
194
195        Idfr = num.array(Idfr,num.int)
196        Idgr = num.array(Idgr,num.int)
197
198        full_send_dict[(processor-1)%numproc]  = [Idfl, Idfl]
199        ghost_recv_dict[(processor-1)%numproc] = [Idgl, Idgl]
200        full_send_dict[(processor+1)%numproc]  = [Idfr, Idfr]
201        ghost_recv_dict[(processor+1)%numproc] = [Idgr, Idgr]
202
203
204
205    #print full_send_dict
206    #print ghost_recv_dict       
207   
208    return  points, elements, boundary, full_send_dict, ghost_recv_dict
209
210
211
212def rectangular_periodic(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
213
214
215    """Setup a rectangular grid of triangles
216    with m+1 by n+1 grid points
217    and side lengths len1, len2. If side lengths are omitted
218    the mesh defaults to the unit square.
219
220    len1: x direction (left to right)
221    len2: y direction (bottom to top)
222
223    Return to lists: points and elements suitable for creating a Mesh or
224    FVMesh object, e.g. Mesh(points, elements)
225    """
226
227    delta1 = float(len1)/m
228    delta2 = float(len2)/n
229
230    #Calculate number of points
231    Np = (m+1)*(n+1)
232
233    class VIndex:
234
235        def __init__(self, n,m):
236            self.n = n
237            self.m = m
238
239        def __call__(self, i,j):
240            return j+i*(self.n+1)
241
242    class EIndex:
243
244        def __init__(self, n,m):
245            self.n = n
246            self.m = m
247
248        def __call__(self, i,j):
249            return 2*(j+i*self.n)
250
251
252    I = VIndex(n,m)
253    E = EIndex(n,m)
254
255    points = num.zeros( (Np,2), num.float)
256
257    for i in range(m+1):
258        for j in range(n+1):
259
260            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
261
262    #Construct 2 triangles per rectangular element and assign tags to boundary
263    #Calculate number of triangles
264    Nt = 2*m*n
265
266
267    elements = num.zeros( (Nt,3), num.int)
268    boundary = {}
269    ghosts = {}
270    nt = -1
271    for i in range(m):
272        for j in range(n):
273
274            i1 = I(i,j+1)
275            i2 = I(i,j)
276            i3 = I(i+1,j+1)
277            i4 = I(i+1,j)
278
279            #Lower Element
280            nt = E(i,j)
281            if i == m-1:
282                ghosts[nt] = E(1,j)
283            if i == 0:
284                ghosts[nt] = E(m-2,j)
285
286            if j == n-1:
287                ghosts[nt] = E(i,1)
288
289            if j == 0:
290                ghosts[nt] = E(i,n-2)
291
292            if i == m-1:
293                if processor == numproc-1:
294                    boundary[nt, 2] = 'right'
295                else:
296                    boundary[nt, 2] = 'ghost'
297            if j == 0:
298                boundary[nt, 1] = 'bottom'
299            elements[nt,:] = [i4,i3,i2]
300
301            #Upper Element
302            nt = E(i,j)+1
303            if i == m-1:
304                ghosts[nt] = E(1,j)+1
305            if i == 0:
306                ghosts[nt] = E(m-2,j)+1
307
308            if j == n-1:
309                ghosts[nt] = E(i,1)+1
310
311            if j == 0:
312                ghosts[nt] = E(i,n-2)+1
313
314            if i == 0:
315                if processor == 0:
316                    boundary[nt, 2] = 'left'
317                else:
318                    boundary[nt, 2] = 'ghost'
319            if j == n-1:
320                boundary[nt, 1] = 'top'
321            elements[nt,:] = [i1,i2,i3]
322
323    #bottom left
324    nt = E(0,0)
325    nf = E(m-2,n-2)
326    ghosts[nt]   = nf
327    ghosts[nt+1] = nf+1
328
329    #bottom right
330    nt = E(m-1,0)
331    nf = E(1,n-2)
332    ghosts[nt]   = nf
333    ghosts[nt+1] = nf+1
334
335    #top left
336    nt = E(0,n-1)
337    nf = E(m-2,1)
338    ghosts[nt]   = nf
339    ghosts[nt+1] = nf+1
340
341    #top right
342    nt = E(m-1,n-1)
343    nf = E(1,1)
344    ghosts[nt]   = nf
345    ghosts[nt+1] = nf+1
346
347    return points, elements, boundary, ghosts
348
349def rectangular_periodic_lr(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
350
351
352    """Setup a rectangular grid of triangles
353    with m+1 by n+1 grid points
354    and side lengths len1, len2. If side lengths are omitted
355    the mesh defaults to the unit square.
356
357    len1: x direction (left to right)
358    len2: y direction (bottom to top)
359
360    Return to lists: points and elements suitable for creating a Mesh or
361    Domain object, e.g. Mesh(points, elements)
362    """
363
364    delta1 = float(len1)/m
365    delta2 = float(len2)/n
366
367    #Calculate number of points
368    Np = (m+1)*(n+1)
369
370    class VIndex:
371
372        def __init__(self, n,m):
373            self.n = n
374            self.m = m
375
376        def __call__(self, i,j):
377            return j+i*(self.n+1)
378
379    class EIndex:
380
381        def __init__(self, n,m):
382            self.n = n
383            self.m = m
384
385        def __call__(self, i,j):
386            return 2*(j+i*self.n)
387
388
389    I = VIndex(n,m)
390    E = EIndex(n,m)
391
392    points = num.zeros( (Np,2), num.float)
393
394    for i in range(m+1):
395        for j in range(n+1):
396
397            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
398
399    #Construct 2 triangles per rectangular element and assign tags to boundary
400    #Calculate number of triangles
401    Nt = 2*m*n
402
403
404    elements = num.zeros( (Nt,3), num.int)
405    boundary = {}
406    ghosts = {}
407    nt = -1
408    for i in range(m):
409        for j in range(n):
410
411            i1 = I(i,j+1)
412            i2 = I(i,j)
413            i3 = I(i+1,j+1)
414            i4 = I(i+1,j)
415
416            #Lower Element
417            nt = E(i,j)
418            if i == m-1:
419                ghosts[nt] = E(1,j)
420            if i == 0:
421                ghosts[nt] = E(m-2,j)
422
423            if i == m-1:
424                if processor == numproc-1:
425                    boundary[nt, 2] = 'right'
426                else:
427                    boundary[nt, 2] = 'ghost'
428            if j == 0:
429                boundary[nt, 1] = 'bottom'
430            elements[nt,:] = [i4,i3,i2]
431
432            #Upper Element
433            nt = E(i,j)+1
434            if i == m-1:
435                ghosts[nt] = E(1,j)+1
436            if i == 0:
437                ghosts[nt] = E(m-2,j)+1
438
439            if i == 0:
440                if processor == 0:
441                    boundary[nt, 2] = 'left'
442                else:
443                    boundary[nt, 2] = 'ghost'
444            if j == n-1:
445                boundary[nt, 1] = 'top'
446            elements[nt,:] = [i1,i2,i3]
447
448
449    return points, elements, boundary, ghosts
Note: See TracBrowser for help on using the repository browser.