source: inundation/ga/storm_surge/pyvolution/test_quantity.py @ 626

Last change on this file since 626 was 590, checked in by duncan, 20 years ago

set_function_values now handles a subset of triangles. Old and new code.
File size: 25.4 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from quantity import *
8from config import epsilon
9from Numeric import allclose, array, ones, Float
10
11       
12class TestCase(unittest.TestCase):
13    def setUp(self):
14        from domain import Domain
15       
16        a = [0.0, 0.0]
17        b = [0.0, 2.0]
18        c = [2.0,0.0]
19        d = [0.0, 4.0]
20        e = [2.0, 2.0]
21        f = [4.0,0.0]
22
23        points = [a, b, c, d, e, f]
24
25        #bac, bce, ecf, dbe
26        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
27
28        self.mesh1 = Domain(points[:3], [elements[0]])       
29        self.mesh1.check_integrity()
30       
31        self.mesh4 = Domain(points, elements)       
32        self.mesh4.check_integrity()
33       
34    def tearDown(self):
35        pass
36        #print "  Tearing down"
37
38
39    def test_creation(self):
40       
41        quantity = Quantity(self.mesh1, [[1,2,3]])
42        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
43
44        try:
45            quantity = Quantity()
46        except:
47            pass
48        else:
49            raise 'Should have raised empty quantity exception'
50
51
52        try:
53            quantity = Quantity([1,2,3])
54        except AssertionError:
55            pass
56        except:
57            raise 'Should have raised "mising mesh object" error'       
58       
59
60    def test_creation_zeros(self):
61       
62        quantity = Quantity(self.mesh1)
63        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
64
65
66        quantity = Quantity(self.mesh4)
67        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
68                                                 [0.,0.,0.], [0.,0.,0.]])       
69
70       
71    def test_interpolation(self):
72        quantity = Quantity(self.mesh1, [[1,2,3]])
73        assert allclose(quantity.centroid_values, [2.0]) #Centroid
74
75        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]]) 
76
77
78    def test_interpolation2(self):
79        quantity = Conserved_quantity(self.mesh4,
80                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
81        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
82
83
84        quantity.extrapolate_second_order()
85
86        assert allclose(quantity.vertex_values, [[2., 2., 2.],
87                                               [3.+2./3, 6.+2./3, 4.+2./3],
88                                               [7.5, 0.5, 1.],
89                                               [-5, -2.5, 7.5]])
90       
91        #print quantity.edge_values       
92        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
93                                               [5., 5., 5.],
94                                               [4.5, 4.5, 0.],
95                                               [3.0, -1.5, -1.5]])
96
97    def test_boundary_allocation(self):
98        quantity = Conserved_quantity(self.mesh4,
99                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
100
101        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
102
103
104    def test_set_values(self):
105        quantity = Quantity(self.mesh4)
106
107
108        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
109                            location = 'vertices')
110        assert allclose(quantity.vertex_values,
111                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
112        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
113        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
114                                               [5., 5., 5.],
115                                               [4.5, 4.5, 0.],
116                                               [3.0, -1.5, -1.5]])
117
118
119        #Test default
120        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
121        assert allclose(quantity.vertex_values,
122                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
123        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
124        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
125                                               [5., 5., 5.],
126                                               [4.5, 4.5, 0.],
127                                               [3.0, -1.5, -1.5]])
128
129        #Test centroids
130        quantity.set_values([1,2,3,4], location = 'centroids')
131        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
132
133        #Test edges
134        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
135                            location = 'edges')
136        assert allclose(quantity.edge_values,
137                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
138
139        #Test exceptions
140        try:
141            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
142                                location = 'bas kamel tuba')
143        except:
144            pass
145
146
147        try:
148            quantity.set_values([[1,2,3], [0,0,9]])
149        except AssertionError:
150            pass
151        except:
152            raise 'should have raised Assertionerror'
153
154
155
156    def test_set_values_const(self):
157        quantity = Quantity(self.mesh4)
158
159        quantity.set_values(1.0, location = 'vertices')
160        assert allclose(quantity.vertex_values,
161                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
162        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
163        assert allclose(quantity.edge_values, [[1, 1, 1],
164                                               [1, 1, 1],
165                                               [1, 1, 1],
166                                               [1, 1, 1]])
167
168
169        quantity.set_values(2.0, location = 'centroids')
170        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
171
172        quantity.set_values(3.0, location = 'edges')
173        assert allclose(quantity.edge_values, [[3, 3, 3],
174                                               [3, 3, 3],
175                                               [3, 3, 3],
176                                               [3, 3, 3]])       
177
178
179    def test_set_values_func(self):
180        quantity = Quantity(self.mesh4)
181
182        def f(x, y):
183            return x+y
184
185        quantity.set_values(f, location = 'vertices')
186        #print "quantity.vertex_values",quantity.vertex_values
187        assert allclose(quantity.vertex_values,
188                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])       
189        assert allclose(quantity.centroid_values,
190                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
191        assert allclose(quantity.edge_values,
192                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])               
193
194       
195        quantity.set_values(f, location = 'centroids')
196        assert allclose(quantity.centroid_values,
197                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])       
198
199
200
201
202    def test_set_vertex_values(self):
203        quantity = Quantity(self.mesh4)
204
205
206        quantity.set_vertex_values([0,1,2,3,4,5])
207
208
209        assert allclose(quantity.vertex_values,
210                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
211       
212        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid
213       
214        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
215                                               [3., 2.5, 1.5],
216                                               [3.5, 4.5, 3.],
217                                               [2.5, 3.5, 2]])
218
219
220
221
222    def test_set_vertex_values_using_general_interface(self):
223        quantity = Quantity(self.mesh4)
224
225
226        quantity.set_values([0,1,2,3,4,5])
227
228
229        assert allclose(quantity.vertex_values,
230                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
231       
232        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3]) #Centroid
233       
234        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
235                                               [3., 2.5, 1.5],
236                                               [3.5, 4.5, 3.],
237                                               [2.5, 3.5, 2]])
238
239
240
241
242    def test_gradient(self):
243        quantity = Conserved_quantity(self.mesh4)
244
245        #Set up for a gradient of (3,0) at mid triangle
246        quantity.set_values([2.0, 4.0, 8.0, 2.0],
247                            location = 'centroids')
248
249
250
251        #Gradients
252        a, b = quantity.compute_gradients()
253
254        #gradient bewteen t0 and t1 is undefined as det==0
255        assert a[0] == 0.0
256        assert b[0] == 0.0
257        #The others are OK
258        for i in range(1,4):
259            assert a[i] == 3.0
260            assert b[i] == 0.0
261
262
263        quantity.extrapolate_second_order()
264
265        #print quantity.vertex_values
266        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
267                                                 [0., 6.,  6.],
268                                                 [6., 6., 12.],
269                                                 [0., 0.,  6.]])
270
271               
272    def test_second_order_extrapolation2(self):
273        quantity = Conserved_quantity(self.mesh4)       
274
275        #Set up for a gradient of (3,1), f(x) = 3x+y
276        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
277                            location = 'centroids')
278       
279        #Gradients
280        a, b = quantity.compute_gradients()
281
282        #gradient bewteen t0 and t1 is undefined as det==0
283        assert a[0] == 0.0
284        assert b[0] == 0.0
285        #The others are OK
286        for i in range(1,4):
287            assert allclose(a[i], 3.0)
288            assert allclose(b[i], 1.0)
289
290        quantity.extrapolate_second_order()
291       
292        #print quantity.vertex_values
293        assert allclose(quantity.vertex_values[1,0], 2.0)
294        assert allclose(quantity.vertex_values[1,1], 6.0)       
295        assert allclose(quantity.vertex_values[1,2], 8.0)
296
297
298
299    def test_first_order_extrapolator(self):
300        quantity = Conserved_quantity(self.mesh4)               
301
302        #Test centroids
303        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
304        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
305
306        #Extrapolate
307        quantity.extrapolate_first_order()
308
309        #Check vertices but not edge values
310        assert allclose(quantity.vertex_values,
311                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
312
313
314    def test_second_order_extrapolator(self):
315        quantity = Conserved_quantity(self.mesh4)                       
316
317        #Set up for a gradient of (3,0) at mid triangle
318        quantity.set_values([2.0, 4.0, 8.0, 2.0],
319                            location = 'centroids')
320       
321
322
323        quantity.extrapolate_second_order()
324        quantity.limit()
325
326
327        #Assert that central triangle is limited by neighbours
328        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
329        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
330       
331        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
332        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
333       
334        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
335        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
336
337       
338        #Assert that quantities are conserved
339        from Numeric import sum
340        for k in range(quantity.centroid_values.shape[0]):
341            assert allclose (quantity.centroid_values[k],
342                             sum(quantity.vertex_values[k,:])/3)
343       
344       
345
346       
347
348    def test_limiter(self):
349        quantity = Conserved_quantity(self.mesh4)
350
351        #Create a deliberate overshoot (e.g. from gradient computation)
352        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
353
354
355        #Limit
356        quantity.limit()
357
358        #Assert that central triangle is limited by neighbours
359        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
360        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
361       
362        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
363        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
364       
365        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
366        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
367
368
369       
370        #Assert that quantities are conserved
371        from Numeric import sum
372        for k in range(quantity.centroid_values.shape[0]):
373            assert allclose (quantity.centroid_values[k],
374                             sum(quantity.vertex_values[k,:])/3)
375       
376
377
378    def test_distribute_first_order(self):
379        quantity = Conserved_quantity(self.mesh4)       
380
381        #Test centroids
382        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
383        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
384
385
386        #Extrapolate
387        quantity.extrapolate_first_order()
388
389        #Interpolate
390        quantity.interpolate_from_vertices_to_edges()       
391
392        assert allclose(quantity.vertex_values,
393                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
394        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
395                                               [3,3,3], [4, 4, 4]])
396
397
398
399    def test_update_explicit(self):
400        quantity = Conserved_quantity(self.mesh4)
401
402        #Test centroids
403        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
404        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
405
406        #Set explicit_update
407        quantity.explicit_update = array( [1.,1.,1.,1.] )
408
409        #Update with given timestep
410        quantity.update(0.1)
411
412        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
413        assert allclose( quantity.centroid_values, x)
414
415    def test_update_semi_implicit(self):
416        quantity = Conserved_quantity(self.mesh4)
417
418        #Test centroids
419        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
420        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
421
422        #Set semi implicit update
423        quantity.semi_implicit_update = array([1.,1.,1.,1.])
424
425        #Update with given timestep
426        timestep = 0.1
427        quantity.update(timestep)
428
429        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
430        denom = ones(4, Float)-timestep*sem
431       
432        x = array([1, 2, 3, 4])/denom
433        assert allclose( quantity.centroid_values, x)
434
435
436    def test_both_updates(self):
437        quantity = Conserved_quantity(self.mesh4)
438
439        #Test centroids
440        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
441        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
442
443        #Set explicit_update
444        quantity.explicit_update = array( [4.,3.,2.,1.] )
445       
446        #Set semi implicit update
447        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
448
449        #Update with given timestep
450        timestep = 0.1       
451        quantity.update(0.1)
452
453        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
454        denom = ones(4, Float)-timestep*sem
455
456        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
457        x /= denom
458        assert allclose( quantity.centroid_values, x)       
459
460       
461
462
463    #Test smoothing   
464    def test_smoothing(self):
465
466        from mesh_factory import rectangular
467        from shallow_water import Domain, Transmissive_boundary
468        from Numeric import zeros, Float
469        from util import mean
470
471        #Create basic mesh
472        points, vertices, boundary = rectangular(2, 2)
473
474        #Create shallow water domain
475        domain = Domain(points, vertices, boundary)
476        domain.default_order=2
477        domain.reduction = mean
478
479       
480        #Set some field values
481        domain.set_quantity('elevation', lambda x,y: x)       
482        domain.set_quantity('friction', 0.03)
483
484
485        ######################
486        # Boundary conditions
487        B = Transmissive_boundary(domain)
488        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
489       
490
491        ######################
492        #Initial condition - with jumps
493
494        bed = domain.quantities['elevation'].vertex_values
495        level = zeros(bed.shape, Float)
496
497        h = 0.03
498        for i in range(level.shape[0]):
499            if i % 2 == 0:           
500                level[i,:] = bed[i,:] + h
501            else:
502                level[i,:] = bed[i,:]
503               
504        domain.set_quantity('level', level)
505
506        level = domain.quantities['level']
507
508        #Get smoothed level
509        A, V = level.get_vertex_values(xy=False, smooth=True)
510        Q = level.vertex_values
511
512       
513        assert A.shape[0] == 9 
514        assert V.shape[0] == 8 
515        assert V.shape[1] == 3                 
516       
517        #First four points
518        assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 
519        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
520        assert allclose(A[2], Q[3,0])
521        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
522
523        #Center point
524        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
525                               Q[5,0] + Q[6,2] + Q[7,1])/6)
526                                 
527
528        #Check V
529        assert allclose(V[0,:], [3,4,0])
530        assert allclose(V[1,:], [1,0,4])
531        assert allclose(V[2,:], [4,5,1])                       
532        assert allclose(V[3,:], [2,1,5])
533        assert allclose(V[4,:], [6,7,3])       
534        assert allclose(V[5,:], [4,3,7])
535        assert allclose(V[6,:], [7,8,4])
536        assert allclose(V[7,:], [5,4,8])                       
537
538        #Get smoothed level with XY
539        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=True)
540
541        assert allclose(A, A1)
542        assert allclose(V, V1)       
543
544        #Check XY
545        assert allclose(X[4], 0.5)
546        assert allclose(Y[4], 0.5)
547
548        assert allclose(X[7], 1.0)
549        assert allclose(Y[7], 0.5)                 
550
551
552
553
554    def test_vertex_values_no_smoothing(self):
555
556        from mesh_factory import rectangular
557        from shallow_water import Domain, Transmissive_boundary
558        from Numeric import zeros, Float
559        from util import mean
560
561       
562        #Create basic mesh
563        points, vertices, boundary = rectangular(2, 2)
564
565        #Create shallow water domain
566        domain = Domain(points, vertices, boundary)
567        domain.default_order=2
568        domain.reduction = mean
569
570       
571        #Set some field values
572        domain.set_quantity('elevation', lambda x,y: x)       
573        domain.set_quantity('friction', 0.03)
574
575
576        ######################
577        #Initial condition - with jumps
578
579        bed = domain.quantities['elevation'].vertex_values
580        level = zeros(bed.shape, Float)
581
582        h = 0.03
583        for i in range(level.shape[0]):
584            if i % 2 == 0:           
585                level[i,:] = bed[i,:] + h
586            else:
587                level[i,:] = bed[i,:]
588               
589        domain.set_quantity('level', level)
590
591        #Get level
592        level = domain.quantities['level']       
593        A, V = level.get_vertex_values(xy=False, smooth=False)
594        Q = level.vertex_values.flat
595
596        for k in range(8):
597            assert allclose(A[k], Q[k])
598
599           
600        for k in range(8):
601            assert V[k, 0] == 3*k
602            assert V[k, 1] == 3*k+1
603            assert V[k, 2] == 3*k+2           
604           
605
606
607        X, Y, A1, V1 = level.get_vertex_values(xy=True, smooth=False)
608
609       
610        assert allclose(A, A1)
611        assert allclose(V, V1)       
612
613        #Check XY
614        assert allclose(X[1], 0.5)
615        assert allclose(Y[1], 0.5)
616        assert allclose(X[4], 0.0)
617        assert allclose(Y[4], 0.0)
618        assert allclose(X[12], 1.0)
619        assert allclose(Y[12], 0.0)               
620       
621       
622 
623    def set_array_values_by_index(self):
624
625        from mesh_factory import rectangular
626        from shallow_water import Domain
627        from Numeric import zeros, Float
628       
629        #Create basic mesh
630        points, vertices, boundary = rectangular(1, 1)
631
632        #Create shallow water domain
633        domain = Domain(points, vertices, boundary)
634        #print "domain.number_of_elements ",domain.number_of_elements
635        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
636        value = [7]
637        indexes = [1]
638        quantity.set_array_values_by_index(value,
639                                           location = 'centroids',
640                                           indexes = indexes)
641        #print "quantity.centroid_values",quantity.centroid_values
642       
643        assert allclose(quantity.centroid_values, [1,7])
644
645        quantity.set_array_values([15,20,25], indexes = indexes)
646        assert allclose(quantity.centroid_values, [1,20])
647
648        quantity.set_array_values([15,20,25], indexes = indexes)
649        assert allclose(quantity.centroid_values, [1,20])
650       
651    def test_setting_some_vertex_values(self):
652        """
653        set values based on triangle lists.
654        """
655        from mesh_factory import rectangular
656        from shallow_water import Domain
657        from Numeric import zeros, Float
658       
659        #Create basic mesh
660        points, vertices, boundary = rectangular(1, 3)
661
662        #Create shallow water domain
663        domain = Domain(points, vertices, boundary)
664        #print "domain.number_of_elements ",domain.number_of_elements
665        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
666                                    [4,4,4],[5,5,5],[6,6,6]])
667        value = [7]
668        indexes = [1]
669        quantity.set_values(value,
670                                  location = 'centroids',
671                                  indexes = indexes)
672        #print "quantity.centroid_values",quantity.centroid_values
673        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
674       
675        value = [[15,20,25]]
676        quantity.set_values(value, indexes = indexes)
677        #print "1 quantity.vertex_values",quantity.vertex_values
678        assert allclose(quantity.vertex_values[1], value[0])
679
680        # FIXMEDSG this seems wrong -DSG
681        # It's right.  The values is an array, indexed by the element id's
682        # given in indexes.  It sets values per triangle.
683        values = [10,100,50]
684        quantity.set_values(values, indexes = [0,1,5]) # , location = 'centroids'
685        #print "2 quantity.vertex_values",quantity.vertex_values
686        assert allclose(quantity.vertex_values[0], [10,10,10])
687        assert allclose(quantity.vertex_values[5], [50,50,50])
688        #quantity.interpolate()
689        #print "quantity.centroid_values",quantity.centroid_values
690        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
691
692        values = [[31,30,29],[400,400,400],[1000,999,998]]
693        quantity.set_values(values, indexes = [3,3,5])
694        quantity.interpolate()
695        assert allclose(quantity.centroid_values, [10,100,3,400,5,999])
696
697        values = [[1,1,1],[2,2,2],[3,3,3],
698                                    [4,4,4],[5,5,5],[6,6,6]]
699        quantity.set_values(values)
700       
701        # testing the standard set values by vertex
702        # indexed by vertex_id in general_mesh.coordinates
703        values = [0,1,2,3,4,5,6,7]
704       
705        quantity.set_values(values)
706        #print "1 quantity.vertex_values",quantity.vertex_values
707        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
708                                                [ 1.,  0.,  5.],
709                                                [ 5.,  6.,  1.],
710                                                [ 2.,  1.,  6.],
711                                                [ 6.,  7.,  2.],
712                                                [ 3.,  2.,  7.]])
713       
714    def test_getting_some_vertex_values(self):
715        """
716        get values based on triangle lists.
717        """
718        from mesh_factory import rectangular
719        from shallow_water import Domain
720        from Numeric import zeros, Float
721       
722        #Create basic mesh
723        points, vertices, boundary = rectangular(1, 3)
724
725        #print "points",points
726        #print "vertices",vertices
727        #print "boundary",boundary
728
729        #Create shallow water domain
730        domain = Domain(points, vertices, boundary)
731        #print "domain.number_of_elements ",domain.number_of_elements
732        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
733                                    [4,4,4],[5,5,5],[6,6,6]])
734        value = [7]
735        indexes = [1]
736        quantity.set_values(value,
737                                  location = 'centroids',
738                                  indexes = indexes)
739        #print "quantity.centroid_values",quantity.centroid_values
740        #print "quantity.get_values(location = 'centroids') ",\
741        #      quantity.get_values(location = 'centroids')
742        assert allclose(quantity.centroid_values,
743                        quantity.get_values(location = 'centroids'))
744
745       
746        value = [[15,20,25]]
747        quantity.set_values(value, indexes = indexes)
748        #print "1 quantity.vertex_values",quantity.vertex_values
749        assert allclose(quantity.vertex_values, quantity.get_values())
750
751        assert allclose(quantity.edge_values,
752                        quantity.get_values(location = 'edges'))
753
754        # get a subset of elements
755        subset = quantity.get_values(location='centroids', indexes=[0,5])
756        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
757        assert allclose(subset, answer)
758
759
760        subset = quantity.get_values(location='edges', indexes=[0,5])
761        answer = [quantity.edge_values[0],quantity.edge_values[5]]
762        #print "subset",subset
763        #print "answer",answer
764        assert allclose(subset, answer)
765
766        subset = quantity.get_values( indexes=[1,5])
767        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
768        #print "subset",subset
769        #print "answer",answer
770        assert allclose(subset, answer)
771
772       
773       
774#-------------------------------------------------------------
775if __name__ == "__main__":
776    suite = unittest.makeSuite(TestCase,'test')
777    runner = unittest.TextTestRunner()
778    runner.run(suite)
779
780
Note: See TracBrowser for help on using the repository browser.