source: trunk/anuga_work/anuga_cuda/src/utilities/sort_domain.py

Last change on this file was 9017, checked in by steve, 11 years ago

Adding in Zhe (John) Weng's anuga_cuda code as obtained from googlecode ​https://code.google.com/p/anuga-cuda
File size: 28.9 KB
Line 
1#!/usr/bin/env python
2
3def swap_domain(domain, i, j, k):
4    """Swap information of two mesh entities.
5   
6    This is used to ensure the calculation sequence in parallel.
7    """
8
9    neighbours = domain.neighbours[i]
10   
11    s_neighbours = domain.surrogate_neighbours[i]
12    normals = domain.normals[i]
13    edgelengths = domain.edgelengths[i]
14
15    vertex_coordinates = domain.vertex_coordinates
16    edge_coordinates = domain.edge_coordinates
17
18
19    stage_edge_values = domain.quantities['stage'].edge_values[i]
20    xmom_edge_values = domain.quantities['xmomentum'].edge_values[i]
21    ymom_edge_values = domain.quantities['ymomentum'].edge_values[i]
22    bed_edge_values = domain.quantities['elevation'].edge_values[i]
23
24    stage_vertex_values = domain.quantities['stage'].vertex_values[i]
25    xmom_vertex_values = domain.quantities['xmomentum'].vertex_values[i]
26    ymom_vertex_values = domain.quantities['ymomentum'].vertex_values[i]
27    bed_vertex_values = domain.quantities['elevation'].vertex_values[i]
28   
29
30    neighbours[j], neighbours[k] = neighbours[k], neighbours[j]
31
32    s_neighbours[j], s_neighbours[k] = s_neighbours[k], s_neighbours[j]
33
34    normals[j*2], normals[k*2] = normals[k*2], normals[j*2]
35
36    normals[j*2+1], normals[k*2+1] = normals[k*2+1], normals[j*2+1]
37
38    edgelengths[j], edgelengths[k] = edgelengths[k], edgelengths[j]
39
40    vertex_coordinates[i*3+j][0], vertex_coordinates[i*3+k][0] = vertex_coordinates[i*3+k][0], vertex_coordinates[i*3+j][0] 
41    vertex_coordinates[i*3+j][1], vertex_coordinates[i*3+k][1] = vertex_coordinates[i*3+k][1], vertex_coordinates[i*3+j][1] 
42
43    edge_coordinates[i*3+j][0], edge_coordinates[i*3+k][0] = edge_coordinates[i*3+k][0], edge_coordinates[i*3+j][0] 
44    edge_coordinates[i*3+j][1], edge_coordinates[i*3+k][1] = edge_coordinates[i*3+k][1], edge_coordinates[i*3+j][1] 
45
46    stage_edge_values[j], stage_edge_values[k] = stage_edge_values[k], stage_edge_values[j]
47    xmom_edge_values[j], xmom_edge_values[k] = xmom_edge_values[k], xmom_edge_values[j]
48    ymom_edge_values[j], ymom_edge_values[k] = ymom_edge_values[k], ymom_edge_values[j]
49    bed_edge_values[j], bed_edge_values[k] = bed_edge_values[k], bed_edge_values[j]
50
51    stage_vertex_values[j], stage_vertex_values[k] = stage_vertex_values[k], stage_vertex_values[j]
52    xmom_vertex_values[j], xmom_vertex_values[k] = xmom_vertex_values[k], xmom_vertex_values[j]
53    ymom_vertex_values[j], ymom_vertex_values[k] = ymom_vertex_values[k], ymom_vertex_values[j]
54    bed_vertex_values[j], bed_vertex_values[k] = bed_vertex_values[k], bed_vertex_values[j]
55
56
57
58def sort_domain(domain):
59    """Sort the neighbour order of mesh entities
60
61    This is used to ensure the calculation sequence in paralle as same as
62    that in sequential
63    """
64
65    for i in range(domain.number_of_elements):
66       
67        neighbours = domain.neighbours[i]
68
69        if neighbours[2]>=0 and neighbours[2] < i and \
70                (neighbours[1]<0 or neighbours[2]<neighbours[1]):
71            swap_domain(domain, i, 2, 1)
72
73        if neighbours[1]>=0 and neighbours[1]<i and \
74                (neighbours[0]<0 or neighbours[1]<neighbours[0]):
75            swap_domain(domain, i, 0, 1)
76
77        if neighbours[2]>=0 and neighbours[2]<i and \
78               (neighbours[1]<0 or neighbours[2]<neighbours[1]):
79            swap_domain(domain, i, 2, 1)
80   
81   
82    for i in range(domain.number_of_elements):
83        n1, n2, n3 = domain.neighbours[i]
84       
85        for index in range(3):
86            if domain.neighbours[n1][index] == i:
87                domain.neighbour_edges[n1][index] = 0
88            if domain.neighbours[n2][index] == i:
89                domain.neighbour_edges[n2][index] = 1
90            if domain.neighbours[n3][index] == i:
91                domain.neighbour_edges[n3][index] = 2               
92                   
93
94
95def sort_domain_check(domain2):
96    """Check the sorted domain information."""
97
98    from anuga_cuda.merimbula_data.generate_domain import domain_create
99
100    domain1 = domain_create()
101
102    for k in range(domain1.number_of_elements):
103        b = [0,1,2]
104        spe_bubble_sort(b, domain1.neighbours[k], k)
105        for i in range(3):
106            if domain1.neighbours[k][b[i]] != domain2.neighbours[k][i]:
107                print "###### tri: %ld, edge: %d - %d " % (k, b[i], i)
108                continue
109           
110            if domain1.surrogate_neighbours[k][ b[i] ] != domain2.surrogate_neighbours[k][i]:
111                print "tri: %ld, edge: %d - %d @ surrogate_neighbours " % (k, b[i], i)
112               
113            if domain1.normals[k][ b[i]*2] != domain2.normals[k][i*2] or\
114                     domain1.normals[k][ b[i]*2+1] != domain2.normals[k][i*2+1]:
115                print "tri: %ld, edge: %ld - %d @ normals " % (k, b[i], i)
116
117            if domain1.edgelengths[k][ b[i] ] != domain2.edgelengths[k][i]:
118                print "tri: %ld, edge: %ld - %d @ edgelengths " % (k, b[i], i)
119
120
121            if domain1.vertex_coordinates[k*3+b[i]][0] != domain2.vertex_coordinates[k*3+i][0] or\
122                     domain1.vertex_coordinates[k*3+b[i]][1] != domain2.vertex_coordinates[k*3+i][1]:
123                print "tri: %ld, edge: %ld - %d @ vertex_coordinates " % (k, b[i], i)
124
125            if domain1.edge_coordinates[k*3+b[i]][0] != domain2.edge_coordinates[k*3+i][0] or\
126                     domain1.edge_coordinates[k*3+b[i]][1] != domain2.edge_coordinates[k*3+i][1]:
127                print "tri: %ld, edge: %ld - %d @ edge_coordinates " % (k, b[i], i)
128
129
130            if domain1.quantities['stage'].edge_values[k][ b[i] ] != \
131                    domain2.quantities['stage'].edge_values[k][i]:
132                print "tri: %ld, edge: %ld - %d @ stage_edge_values " % (k, b[i], i)
133   
134            if domain1.quantities['xmomentum'].edge_values[k][ b[i] ] != \
135                    domain2.quantities['xmomentum'].edge_values[k][i]:
136                print "tri: %ld, edge: %ld - %d @ xmom_edge_values " % (k, b[i], i)
137
138            if domain1.quantities['ymomentum'].edge_values[k][ b[i] ] != \
139                    domain2.quantities['ymomentum'].edge_values[k][i]:
140                print "tri: %ld, edge: %ld - %d @ ymom_edge_values " % (k, b[i], i)
141
142            if domain1.quantities['elevation'].edge_values[k][ b[i] ] != \
143                    domain2.quantities['elevation'].edge_values[k][i]:
144                print "tri: %ld, edge: %ld - %d @ bed_edge_values " % (k, b[i], i)
145
146
147            if domain1.quantities['stage'].vertex_values[k][ b[i] ] != \
148                    domain2.quantities['stage'].vertex_values[k][i]:
149                print "tri: %ld, vertex: %ld - %d @ stage_vertex_values " % (k, b[i], i)
150   
151            if domain1.quantities['xmomentum'].vertex_values[k][ b[i] ] != \
152                    domain2.quantities['xmomentum'].vertex_values[k][i]:
153                print "tri: %ld, vertex: %ld - %d @ xmom_vertex_values " % (k, b[i], i)
154
155            if domain1.quantities['ymomentum'].vertex_values[k][ b[i] ] != \
156                    domain2.quantities['ymomentum'].vertex_values[k][i]:
157                print "tri: %ld, vertex: %ld - %d @ ymom_vertex_values " % (k, b[i], i)
158
159            if domain1.quantities['elevation'].vertex_values[k][ b[i] ] != \
160                    domain2.quantities['elevation'].vertex_values[k][i]:
161                print "tri: %ld, vertex: %ld - %d @ bed_vertex_values " % (k, b[i], i)
162
163
164
165def rearrange_domain(domain, sort_flag=True, spare_domain=None):
166    """Rearrange the mesh information stored in memory.
167   
168    This will help to better meet the requirements of coalesced memory
169    access.
170    """
171
172    import copy 
173
174    if spare_domain == None:
175        new_dom = copy.deepcopy(domain)
176       
177        if sort_flag:
178            sort_domain(new_dom)
179            temp_dom = copy.deepcopy(new_dom)
180        else:
181            temp_dom = domain
182    else:
183        new_dom = domain
184        temp_dom = spare_domain
185
186    N = new_dom.number_of_elements
187
188    neighbours = new_dom.neighbours
189    neighbour_edges = new_dom.neighbour_edges
190    surrogate_neighbours = new_dom.surrogate_neighbours
191    normals = new_dom.normals
192    edgelengths = new_dom.edgelengths
193
194    vertex_coordinates = new_dom.vertex_coordinates
195    edge_coordinates = new_dom.edge_coordinates
196    centroid_coordinates = new_dom.centroid_coordinates
197
198
199    stage_edge_values = new_dom.quantities['stage'].edge_values
200    xmom_edge_values = new_dom.quantities['xmomentum'].edge_values
201    ymom_edge_values = new_dom.quantities['ymomentum'].edge_values
202    bed_edge_values = new_dom.quantities['elevation'].edge_values
203
204    stage_vertex_values = new_dom.quantities['stage'].vertex_values
205    xmom_vertex_values = new_dom.quantities['xmomentum'].vertex_values
206    ymom_vertex_values = new_dom.quantities['ymomentum'].vertex_values
207    bed_vertex_values = new_dom.quantities['elevation'].vertex_values
208
209
210    for i in range(N):
211        neighbours[i/3][i%3] = temp_dom.neighbours[i][0]
212        neighbours[(i+N)/3][(i+N)%3] = temp_dom.neighbours[i][1]
213        neighbours[(i+2*N)/3][(i+2*N)%3] = temp_dom.neighbours[i][2]
214
215        neighbour_edges[i/3][i%3] = temp_dom.neighbour_edges[i][0]
216        neighbour_edges[(i+N)/3][(i+N)%3] = temp_dom.neighbour_edges[i][1]
217        neighbour_edges[(i+2*N)/3][(i+2*N)%3] = temp_dom.neighbour_edges[i][2]
218
219        surrogate_neighbours[i/3][i%3] = \
220                                                                                temp_dom.surrogate_neighbours[i][0]
221        surrogate_neighbours[(i+N)/3][(i+N)%3] = \
222                                                                                temp_dom.surrogate_neighbours[i][1]
223        surrogate_neighbours[(i+2*N)/3][(i+2*N)%3] = \
224                                                                                temp_dom.surrogate_neighbours[i][2]
225
226
227        normals[i/6][i%6] = temp_dom.normals[i][0]
228        normals[(i+N)/6][(i+N)%6] = temp_dom.normals[i][1]
229        normals[(i+2*N)/6][(i+2*N)%6] = temp_dom.normals[i][2]
230        normals[(i+3*N)/6][(i+3*N)%6] = temp_dom.normals[i][3]
231        normals[(i+4*N)/6][(i+4*N)%6] = temp_dom.normals[i][4]
232        normals[(i+5*N)/6][(i+5*N)%6] = temp_dom.normals[i][5]
233
234
235        edgelengths[i/3][i%3] = temp_dom.edgelengths[i][0]
236        edgelengths[(i+N)/3][(i+N)%3] = temp_dom.edgelengths[i][1]
237        edgelengths[(i+2*N)/3][(i+2*N)%3] = temp_dom.edgelengths[i][2]
238
239       
240        #vertex_coordinates[i/2][i%2] = \
241                #                                                       temp_dom.vertex_coordinates[i*3][0]
242        #vertex_coordinates[(i+N)/2][(i+N)%2] = \
243                #                                                       temp_dom.vertex_coordinates[i*3][1]
244        #vertex_coordinates[(i+N*2)/2][(i+N*2)%2] = \
245                #                                                       temp_dom.vertex_coordinates[i*3+1][0]
246        #vertex_coordinates[(i+N*3)/2][(i+N*3)%2] = \
247                #                                                       temp_dom.vertex_coordinates[i*3+1][1]
248        #vertex_coordinates[(i+N*4)/2][(i+N*4)%2] = \
249                #                                                       temp_dom.vertex_coordinates[i*3+2][0]
250        #vertex_coordinates[(i+N*5)/2][(i+N*5)%2] = \
251        #                                                       temp_dom.vertex_coordinates[i*3+2][1]
252
253        vertex_coordinates.flat[i] = \
254                                                                        temp_dom.vertex_coordinates[i*3][0]
255        vertex_coordinates.flat[i+N] = \
256                                                                        temp_dom.vertex_coordinates[i*3][1]
257        vertex_coordinates.flat[i+N*2] = \
258                                                                        temp_dom.vertex_coordinates[i*3+1][0]
259        vertex_coordinates.flat[i+N*3] = \
260                                                                        temp_dom.vertex_coordinates[i*3+1][1]
261        vertex_coordinates.flat[i+N*4] = \
262                                                                        temp_dom.vertex_coordinates[i*3+2][0]
263        vertex_coordinates.flat[i+N*5] = \
264                                                                temp_dom.vertex_coordinates[i*3+2][1]
265
266       
267        edge_coordinates[i/2][i%2] = temp_dom.edge_coordinates[i*3][0]
268        edge_coordinates[(i+N)/2][(i+N)%2] = temp_dom.edge_coordinates[i*3][1]
269        edge_coordinates[(i+N*2)/2][(i+N*2)%2] = temp_dom.edge_coordinates[i*3+1][0]
270        edge_coordinates[(i+N*3)/2][(i+N*3)%2] = temp_dom.edge_coordinates[i*3+1][1]
271        edge_coordinates[(i+N*4)/2][(i+N*4)%2] = temp_dom.edge_coordinates[i*3+2][0]
272        edge_coordinates[(i+N*5)/2][(i+N*5)%2] = temp_dom.edge_coordinates[i*3+2][1]
273
274
275        centroid_coordinates[i/2][i%2] = temp_dom.centroid_coordinates[i][0]
276        centroid_coordinates[(i+N)/2][(i+N)%2] = temp_dom.centroid_coordinates[i][1]
277
278        stage_edge_values[i/3][i%3] = \
279                                                        temp_dom.quantities['stage'].edge_values[i][0]
280        stage_edge_values[(i+N)/3][(i+N)%3] = \
281                                                        temp_dom.quantities['stage'].edge_values[i][1]
282        stage_edge_values[(i+2*N)/3][(i+2*N)%3] = \
283                                                        temp_dom.quantities['stage'].edge_values[i][2]
284
285
286        xmom_edge_values[i/3][i%3] = \
287                                                temp_dom.quantities['xmomentum'].edge_values[i][0]
288        xmom_edge_values[(i+N)/3][(i+N)%3] = \
289                                                temp_dom.quantities['xmomentum'].edge_values[i][1]
290        xmom_edge_values[(i+2*N)/3][(i+2*N)%3] = \
291                        temp_dom.quantities['xmomentum'].edge_values[i][2]
292
293
294        ymom_edge_values[i/3][i%3] = \
295                                                        temp_dom.quantities['ymomentum'].edge_values[i][0]
296        ymom_edge_values[(i+N)/3][(i+N)%3] = \
297                                                        temp_dom.quantities['ymomentum'].edge_values[i][1]
298        ymom_edge_values[(i+2*N)/3][(i+2*N)%3] = \
299                                                        temp_dom.quantities['ymomentum'].edge_values[i][2]
300
301
302        bed_edge_values[i/3][i%3] = \
303                                                        temp_dom.quantities['elevation'].edge_values[i][0]
304        bed_edge_values[(i+N)/3][(i+N)%3] = \
305                                                        temp_dom.quantities['elevation'].edge_values[i][1]
306        bed_edge_values[(i+2*N)/3][(i+2*N)%3] = \
307                                                        temp_dom.quantities['elevation'].edge_values[i][2]
308
309
310
311        stage_vertex_values[i/3][i%3] = \
312                                                        temp_dom.quantities['stage'].vertex_values[i][0]
313        stage_vertex_values[(i+N)/3][(i+N)%3] = \
314                                                        temp_dom.quantities['stage'].vertex_values[i][1]
315        stage_vertex_values[(i+2*N)/3][(i+2*N)%3] = \
316                                                        temp_dom.quantities['stage'].vertex_values[i][2]
317
318
319        xmom_vertex_values[i/3][i%3] = \
320                                                        temp_dom.quantities['xmomentum'].vertex_values[i][0]
321        xmom_vertex_values[(i+N)/3][(i+N)%3] = \
322                                                        temp_dom.quantities['xmomentum'].vertex_values[i][1]
323        xmom_vertex_values[(i+2*N)/3][(i+2*N)%3] = \
324                                                        temp_dom.quantities['xmomentum'].vertex_values[i][2]
325
326
327        ymom_vertex_values[i/3][i%3] = \
328                                                        temp_dom.quantities['ymomentum'].vertex_values[i][0]
329        ymom_vertex_values[(i+N)/3][(i+N)%3] = \
330                                                        temp_dom.quantities['ymomentum'].vertex_values[i][1]
331        ymom_vertex_values[(i+2*N)/3][(i+2*N)%3] = \
332                                                        temp_dom.quantities['ymomentum'].vertex_values[i][2]
333
334
335        bed_vertex_values[i/3][i%3] = \
336                                                        temp_dom.quantities['elevation'].vertex_values[i][0]
337        bed_vertex_values[(i+N)/3][(i+N)%3] = \
338                                                        temp_dom.quantities['elevation'].vertex_values[i][1]
339        bed_vertex_values[(i+2*N)/3][(i+2*N)%3] = \
340                                                        temp_dom.quantities['elevation'].vertex_values[i][2]
341
342    return new_dom
343
344
345
346def check_rearranged_array(a, b, configure=3):
347    """Check rearranged array."""
348
349    import numpy 
350    counter = 0
351    N = a.shape[0]
352    if configure == 2:
353        """ N * 2 """
354        for i in range(N):
355            if a[i][0] != b[i/2][i%2] or \
356                        a[i][1] != b[(i+N)/2][(i+N)%2] :
357                counter +=1
358
359    if configure == 3:
360        """ N * 3 """
361        for i in range(N):
362            if not ( numpy.allclose( a[i][0], b[i/3][i%3] ) and \
363                        numpy.allclose( a[i][1], b[(i+N)/3][(i+N)%3] ) and \
364                numpy.allclose( a[i][2], b[(i+N*2)/3][(i+N*2)%3] ) ):
365                if counter < 5:
366                    print i, a[i], b[i/3][i%3], b[(i+N)/3][(i+N)%3], \
367                            b[(i+N*2)/3][(i+N*2)%3]
368
369                counter += 1
370
371    elif configure == 6:
372        """ N * 6 """
373        for i in range(N):
374            if a[i][0] != b[i/6][i%6] or \
375                a[i][1] != b[(i+N)/6][(i+N)%6] or \
376                    a[i][2] != b[(i+N*2)/6][(i+N*2)%6] or \
377                    a[i][3] != b[(i+N*3)/6][(i+N*3)%6] or \
378                    a[i][4] != b[(i+N*4)/6][(i+N*4)%6] or \
379                    a[i][5] != b[(i+N*5)/6][(i+N*5)%6]:
380
381                        counter +=1
382
383    elif configure == 32:
384        """ N*3 * 2"""
385        N = N/3
386        for i in range(N):
387            if a[i*3][0] != b[i/2][i%2] or \
388                        a[i*3][1] != b[(i+N)/2][(i+N)%2] or \
389                        a[i*3 + 1][0] != b[(i+N*2)/2][(i+N*2)%2] or \
390                        a[i*3 + 1][1] != b[(i+N*3)/2][(i+N*3)%2] or \
391                        a[i*3 + 2][0] != b[(i+N*4)/2][(i+N*4)%2] or \
392                        a[i*3 + 2][1] != b[(i+N*5)/2][(i+N*5)%2] :
393               
394                counter +=1
395   
396    if counter == 0:
397        # This is used to tell apart the False and 0
398        counter = -1
399    return counter
400
401
402   
403def rearrange_domain_check(domain1, domain2):
404    """Check rearranged mesh information"""
405
406    neighbours = domain1.neighbours
407    neighbour_edges = domain1.neighbour_edges
408    surrogate_neighbours = domain1.surrogate_neighbours
409    normals = domain1.normals
410    edgelengths = domain1.edgelengths
411
412    vertex_coordinates = domain1.vertex_coordinates
413    edge_coordinates = domain1.edge_coordinates
414    centroid_coordinates = domain1.centroid_coordinates
415
416
417    stage_edge_values = domain1.quantities['stage'].edge_values
418    xmom_edge_values = domain1.quantities['xmomentum'].edge_values
419    ymom_edge_values = domain1.quantities['ymomentum'].edge_values
420    bed_edge_values = domain1.quantities['elevation'].edge_values
421
422    stage_vertex_values = domain1.quantities['stage'].vertex_values
423    xmom_vertex_values = domain1.quantities['xmomentum'].vertex_values
424    ymom_vertex_values = domain1.quantities['ymomentum'].vertex_values
425    bed_vertex_values = domain1.quantities['elevation'].vertex_values
426
427
428    vc = domain2.vertex_coordinates
429    xc = domain2.edge_coordinates
430    cc = domain2.centroid_coordinates
431
432
433    se = domain2.quantities['stage'].edge_values
434    xe = domain2.quantities['xmomentum'].edge_values
435    ye = domain2.quantities['ymomentum'].edge_values
436    be = domain2.quantities['elevation'].edge_values
437
438    sv = domain2.quantities['stage'].vertex_values
439    xv = domain2.quantities['xmomentum'].vertex_values
440    yv = domain2.quantities['ymomentum'].vertex_values
441    bv = domain2.quantities['elevation'].vertex_values
442
443
444
445    N = domain1.number_of_elements
446
447    cnt_nb = check_rearranged_array(
448            neighbours, domain2.neighbours, 3)
449    cnt_nbedge = check_rearranged_array(
450            neighbour_edges, domain2.neighbour_edges, 3)
451    cnt_surrnb = check_rearranged_array(
452            surrogate_neighbours, domain2.surrogate_neighbours, 3)
453    cnt_norm = check_rearranged_array(
454            normals, domain2.normals, 6)
455    cnt_el = check_rearranged_array(
456            edgelengths, domain2.edgelengths, 3)
457    cnt_vc = check_rearranged_array(
458            vertex_coordinates, domain2.vertex_coordinates, 32)
459    cnt_ec = check_rearranged_array(
460            edge_coordinates, domain2.edge_coordinates, 32)
461   
462    edgeC = edge_coordinates
463    ec_2 = domain2.edge_coordinates.flat
464    for i in range(N):
465        if edgeC[i*3][0] != ec_2[i] or edgeC[i*3][1] != ec_2[i+N] or \
466            edgeC[i*3+1][0]!=ec_2[i+N*2] or edgeC[i*3+1][1]!=ec_2[i+N*3] or\
467            edgeC[i*3+2][0]!=ec_2[i+N*4] or edgeC[i*3+2][1]!=ec_2[i+N*5]:
468            print i
469
470    cnt_cc = check_rearranged_array(
471            centroid_coordinates, domain2.centroid_coordinates, 2)
472    cnt_se = check_rearranged_array( stage_edge_values,
473            domain2.quantities['stage'].edge_values, 3)
474    cnt_xe = check_rearranged_array( xmom_edge_values,
475            domain2.quantities['xmomentum'].edge_values, 3)
476    cnt_ye = check_rearranged_array( ymom_edge_values,
477            domain2.quantities['ymomentum'].edge_values, 3)
478    cnt_be = check_rearranged_array( bed_edge_values,
479            domain2.quantities['elevation'].edge_values, 3)
480    cnt_sv = check_rearranged_array( stage_vertex_values,
481            domain2.quantities['stage'].vertex_values, 3)
482    cnt_xv = check_rearranged_array( xmom_vertex_values,
483            domain2.quantities['xmomentum'].vertex_values, 3)
484    cnt_yv = check_rearranged_array( ymom_vertex_values,
485            domain2.quantities['ymomentum'].vertex_values, 3)
486    cnt_bv = check_rearranged_array( bed_vertex_values,
487            domain2.quantities["elevation"].vertex_values, 3)
488
489    #for i in range(N):
490    #    if neighbours[i][0] != domain2.neighbours[i/3][i%3] or \
491        #               neighbours[i][1] != domain2.neighbours[(i+N)/3][(i+N)%3] or \
492        #               neighbours[i][2] != domain2.neighbours[(i+N*2)/3][(i+N*2)%3]:
493    #        if cnt_nb < 5:
494    #            print i, neighbours[i], domain2.neighbours[i/3][i%3], \
495    #                    domain2.neighbours[(i+N)/3][(i+N)%3], \
496    #                    domain2.neighbours[(i+N*2)/3][(i+N*2)%3]
497    #        cnt_nb+=1
498
499    #    if neighbour_edges[i][0] != domain2.neighbour_edges[i/3][i%3] or \
500        #               neighbour_edges[i][1] != \
501    #            domain2.neighbour_edges[(i+N)/3][(i+N)%3]or\
502        #               neighbour_edges[i][2] != \
503    #            domain2.neighbour_edges[(i+N*2)/3][(i+N*2)%3]:
504        #               cnt_nbedge+=1
505    #   
506    #    if surrogate_neighbours[i][0] != \
507    #            domain2.surrogate_neighbours[i/3][i%3]or\
508        #               surrogate_neighbours[i][1] !=\
509    #            domain2.surrogate_neighbours[(i+N)/3][(i+N)%3]or\
510        #               surrogate_neighbours[i][2] !=\
511    #            domain2.surrogate_neighbours[(i+N*2)/3][(i+N*2)%3]:
512        #               cnt_surrnb +=1
513
514    #    if normals[i][0] != domain2.normals[i/6][i%6] or \
515        #                   normals[i][1] != domain2.normals[(i+N)/6][(i+N)%6] or \
516        #                   normals[i][2] != domain2.normals[(i+N*2)/6][(i+N*2)%6] or \
517        #                   normals[i][3] != domain2.normals[(i+N*3)/6][(i+N*3)%6] or \
518        #                   normals[i][4] != domain2.normals[(i+N*4)/6][(i+N*4)%6] or \
519        #                   normals[i][5] != domain2.normals[(i+N*5)/6][(i+N*5)%6]:
520        #               cnt_norm +=1
521
522    #    if edgelengths[i][0] != domain2.edgelengths[i/3][i%3] or \
523        #               edgelengths[i][1] != domain2.edgelengths[(i+N)/3][(i+N)%3] or \
524        #               edgelengths[i][2] != domain2.edgelengths[(i+N*2)/3][(i+N*2)%3]:
525    #        if cnt_nb < 5:
526    #            print i, edgelengths[i], domain2.edgelengths[i/3][i%3], \
527    #                    domain2.edgelengths[(i+N)/3][(i+N)%3], \
528    #                    domain2.edgelengths[(i+N*2)/3][(i+N*2)%3]
529    #        cnt_el +=1
530
531
532    #    if vertex_coordinates[i*3][0] != \
533    #            domain2.vertex_coordinates[i/2][i%2] or \
534        #               vertex_coordinates[i*3][1] != \
535    #            domain2.vertex_coordinates[(i+N)/2][(i+N)%2] or \
536        #               vertex_coordinates[i*3 + 1][0] != \
537    #            domain2.vertex_coordinates[(i+N*2)/2][(i+N*2)%2] or \
538        #               vertex_coordinates[i*3 + 1][1] != \
539    #            domain2.vertex_coordinates[(i+N*3)/2][(i+N*3)%2] or \
540        #               vertex_coordinates[i*3 + 2][0] != \
541    #            domain2.vertex_coordinates[(i+N*4)/2][(i+N*4)%2] or \
542        #               vertex_coordinates[i*3 + 2][1] != \
543    #            domain2.vertex_coordinates[(i+N*5)/2][(i+N*5)%2] :
544    #        cnt_vc +=1
545
546
547    #    if edge_coordinates[i*3][0] != \
548    #            domain2.edge_coordinates[i/2][i%2] or \
549        #               edge_coordinates[i*3][1] != \
550    #            domain2.edge_coordinates[(i+N)/2][(i+N)%2] or \
551        #               edge_coordinates[i*3 + 1][0] != \
552    #            domain2.edge_coordinates[(i+N*2)/2][(i+N*2)%2] or \
553        #               edge_coordinates[i*3 + 1][1] != \
554    #            domain2.edge_coordinates[(i+N*3)/2][(i+N*3)%2] or \
555        #               edge_coordinates[i*3 + 2][0] != \
556    #            domain2.edge_coordinates[(i+N*4)/2][(i+N*4)%2] or \
557        #               edge_coordinates[i*3 + 2][1] != \
558    #            domain2.edge_coordinates[(i+N*5)/2][(i+N*5)%2] :
559    #        cnt_ec +=1
560
561
562    #    if centroid_coordinates[i][0] !=\
563    #            domain2.centroid_coordinates[i/2][i%2] or \
564        #               centroid_coordinates[i][1] != \
565    #            domain2.centroid_coordinates[(i+N)/2][(i+N)%2] :
566    #        cnt_cc +=1
567
568
569    #    if stage_edge_values[i][0] != \
570    #            domain2.quantities["stage"].edge_values[i/3][i%3] or \
571        #           stage_edge_values[i][1] != \
572    #            domain2.quantities["stage"].edge_values[(i+N)/3][(i+N)%3] \
573    #        or stage_edge_values[i][2] != \
574    #        domain2.quantities["stage"].edge_values[(i+N*2)/3][(i+N*2)%3]:
575    #        if cnt_nb < 5:
576    #            print i, stage_edge_values[i], \
577    #                domain2.quantities["stage"].edge_values[i/3][i%3], \
578    #                domain2.quantities["stage"].edge_values[(i+N)/3][(i+N)%3],\
579    #                domain2.quantities["stage"].edge_values[(i+N*2)/3][(i+N*2)%3]
580    #        cnt_se+=1
581
582    #    if xmom_edge_values[i][0] != \
583    #            domain2.quantities["xmomentum"].edge_values[i/3][i%3] or \
584        #               xmom_edge_values[i][1] != \
585    #            domain2.quantities["xmomentum"].edge_values[(i+N)/3][(i+N)%3] or \
586        #               xmom_edge_values[i][2] != \
587    #            domain2.quantities["xmomentum"].edge_values[(i+N*2)/3][(i+N*2)%3]:
588    #        if cnt_nb < 5:
589    #            print i, xmom_edge_values[i], domain2.quantities["xmomentum"].edge_values[i/3][i%3], domain2.quantities["xmomentum"].edge_values[(i+N)/3][(i+N)%3], domain2.quantities["xmomentum"].edge_values[(i+N*2)/3][(i+N*2)%3]
590    #        cnt_xe+=1
591
592    #    if ymom_edge_values[i][0] != \
593    #        domain2.quantities["ymomentum"].edge_values[i/3][i%3] or \
594        #               ymom_edge_values[i][1] != \
595    #            domain2.quantities["ymomentum"].edge_values[(i+N)/3][(i+N)%3] or \
596        #               ymom_edge_values[i][2] != \
597    #            domain2.quantities["ymomentum"].edge_values[(i+N*2)/3][(i+N*2)%3]:
598    #        if cnt_nb < 5:
599    #            print i, ymom_edge_values[i], domain2.quantities["ymomentum"].edge_values[i/3][i%3], domain2.quantities["ymomentum"].edge_values[(i+N)/3][(i+N)%3], domain2.quantities["ymomentum"].edge_values[(i+N*2)/3][(i+N*2)%3]
600    #        cnt_ye+=1
601
602    #    if bed_edge_values[i][0] != \
603    #        domain2.quantities["elevation"].edge_values[i/3][i%3] or \
604        #               bed_edge_values[i][1] != \
605    #            domain2.quantities["elevation"].edge_values[(i+N)/3][(i+N)%3] or \
606        #               bed_edge_values[i][2] != \
607    #            domain2.quantities["elevation"].edge_values[(i+N*2)/3][(i+N*2)%3]:
608    #        if cnt_nb < 5:
609    #            print i, bed_edge_values[i], domain2.quantities["elevation"].edge_values[i/3][i%3], domain2.quantities["elevation"].edge_values[(i+N)/3][(i+N)%3], domain2.quantities["elevation"].edge_values[(i+N*2)/3][(i+N*2)%3]
610    #        cnt_be+=1
611
612    #    # Vertex values
613
614    #    if stage_vertex_values[i][0] != \
615    #        domain2.quantities["stage"].vertex_values[i/3][i%3] or \
616        #               stage_vertex_values[i][1] != \
617    #            domain2.quantities["stage"].vertex_values[(i+N)/3][(i+N)%3] or \
618        #               stage_vertex_values[i][2] != \
619    #            domain2.quantities["stage"].vertex_values[(i+N*2)/3][(i+N*2)%3]:
620    #        if cnt_nb < 5:
621    #            print i, stage_vertex_values[i], domain2.quantities["stage"].vertex_values[i/3][i%3], domain2.quantities["stage"].vertex_values[(i+N)/3][(i+N)%3], domain2.quantities["stage"].vertex_values[(i+N*2)/3][(i+N*2)%3]
622    #        cnt_sv +=1
623
624    #    if xmom_vertex_values[i][0] != \
625    #        domain2.quantities["xmomentum"].vertex_values[i/3][i%3] or \
626        #               xmom_vertex_values[i][1] != \
627    #            domain2.quantities["xmomentum"].vertex_values[(i+N)/3][(i+N)%3] or \
628        #               xmom_vertex_values[i][2] != \
629    #            domain2.quantities["xmomentum"].vertex_values[(i+N*2)/3][(i+N*2)%3]:
630    #        if cnt_nb < 5:
631    #            print i, xmom_vertex_values[i], domain2.quantities["xmomentum"].vertex_values[i/3][i%3], domain2.quantities["xmomentum"].vertex_values[(i+N)/3][(i+N)%3], domain2.quantities["xmomentum"].vertex_values[(i+N*2)/3][(i+N*2)%3]
632    #        cnt_xv +=1
633
634    #    if ymom_vertex_values[i][0] != \
635    #        domain2.quantities["ymomentum"].vertex_values[i/3][i%3] or \
636        #               ymom_vertex_values[i][1] != \
637    #            domain2.quantities["ymomentum"].vertex_values[(i+N)/3][(i+N)%3] or \
638        #               ymom_vertex_values[i][2] != \
639    #            domain2.quantities["ymomentum"].vertex_values[(i+N*2)/3][(i+N*2)%3]:
640    #        if cnt_nb < 5:
641    #            print i, ymom_vertex_values[i], domain2.quantities["ymomentum"].vertex_values[i/3][i%3], domain2.quantities["ymomentum"].vertex_values[(i+N)/3][(i+N)%3], domain2.quantities["ymomentum"].vertex_values[(i+N*2)/3][(i+N*2)%3]
642    #        cnt_yv +=1
643
644    #    if bed_vertex_values[i][0] != \
645    #        domain2.quantities["elevation"].vertex_values[i/3][i%3] or \
646        #               bed_vertex_values[i][1] != \
647    #            domain2.quantities["elevation"].vertex_values[(i+N)/3][(i+N)%3] or \
648        #               bed_vertex_values[i][2] != \
649    #            domain2.quantities["elevation"].vertex_values[(i+N*2)/3][(i+N*2)%3]:
650    #        if cnt_nb < 5:
651    #            print i, bed_vertex_values[i], domain2.quantities["elevation"].vertex_values[i/3][i%3], domain2.quantities["elevation"].vertex_values[(i+N)/3][(i+N)%3], domain2.quantities["elevation"].vertex_values[(i+N*2)/3][(i+N*2)%3]
652    #        cnt_bv +=1
653
654    if cnt_nb or cnt_nbedge or cnt_surrnb or cnt_norm:
655        print "nb=%d, nb_edge=%d surrnb=%d norm=%d" % \
656            (cnt_nb, cnt_nbedge, cnt_surrnb, cnt_norm)
657    if cnt_el or cnt_vc or cnt_ec or cnt_cc:
658        print "elen=%d, vc=%d ec=%d cc=%d" % \
659            (cnt_el, cnt_vc, cnt_ec, cnt_cc)
660    if cnt_se or cnt_xe or cnt_ye or cnt_be:
661        print "se=%d, xe=%d ye=%d be=%d" % \
662            (cnt_se, cnt_xe, cnt_ye, cnt_be)
663    if cnt_sv or cnt_xv or cnt_yv or cnt_bv:
664        print "sv=%d, xv=%d yv=%d bv=%d" % \
665            (cnt_sv, cnt_xv, cnt_yv, cnt_bv)
666
667
668
669
670if __name__ == "__main__":
671    from anuga_cuda import generate_merimbula_domain
672
673    from anuga_cuda.compute_fluxes.compute_fluxes import spe_bubble_sort
674
675    domain1 = generate_merimbula_domain()
676    domain2 = generate_merimbula_domain()
677    #sort_domain(domain2)
678    #sort_domain_check(domain2)
679       
680    #sort_domain(domain1)
681    print " **** rearrange_domain **** "
682    domain2=rearrange_domain(domain2, False)
683    print " **** check domain ****"
684    rearrange_domain_check(domain1, domain2)
Note: See TracBrowser for help on using the repository browser.