source: inundation/ga/storm_surge/parallel/parallel_meshes.py @ 1472

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