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

Last change on this file since 826 was 773, checked in by ole, 20 years ago

Changed quantity name 'level' to 'stage'
File size: 30.3 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    def test_set_vertex_values(self):
201        quantity = Quantity(self.mesh4)
202        quantity.set_vertex_values([0,1,2,3,4,5])
203
204        assert allclose(quantity.vertex_values,
205                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
206        assert allclose(quantity.centroid_values,
207                        [1., 7./3, 11./3, 8./3]) #Centroid
208        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
209                                               [3., 2.5, 1.5],
210                                               [3.5, 4.5, 3.],
211                                               [2.5, 3.5, 2]])
212
213
214    def test_set_vertex_values_subset(self):
215        quantity = Quantity(self.mesh4)
216        quantity.set_vertex_values([0,1,2,3,4,5])
217        quantity.set_vertex_values([0,20,30,50], indexes = [0,2,3,5])
218
219        assert allclose(quantity.vertex_values,
220                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
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        #Centroid
233        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
234       
235        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
236                                               [3., 2.5, 1.5],
237                                               [3.5, 4.5, 3.],
238                                               [2.5, 3.5, 2]])
239
240
241
242
243    def test_gradient(self):
244        quantity = Conserved_quantity(self.mesh4)
245
246        #Set up for a gradient of (3,0) at mid triangle
247        quantity.set_values([2.0, 4.0, 8.0, 2.0],
248                            location = 'centroids')
249
250
251
252        #Gradients
253        a, b = quantity.compute_gradients()
254
255        #gradient bewteen t0 and t1 is undefined as det==0
256        assert a[0] == 0.0
257        assert b[0] == 0.0
258        #The others are OK
259        for i in range(1,4):
260            assert a[i] == 3.0
261            assert b[i] == 0.0
262
263
264        quantity.extrapolate_second_order()
265
266        #print quantity.vertex_values
267        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
268                                                 [0., 6.,  6.],
269                                                 [6., 6., 12.],
270                                                 [0., 0.,  6.]])
271
272               
273    def test_second_order_extrapolation2(self):
274        quantity = Conserved_quantity(self.mesh4)       
275
276        #Set up for a gradient of (3,1), f(x) = 3x+y
277        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
278                            location = 'centroids')
279       
280        #Gradients
281        a, b = quantity.compute_gradients()
282
283        #gradient bewteen t0 and t1 is undefined as det==0
284        assert a[0] == 0.0
285        assert b[0] == 0.0
286        #The others are OK
287        for i in range(1,4):
288            assert allclose(a[i], 3.0)
289            assert allclose(b[i], 1.0)
290
291        quantity.extrapolate_second_order()
292       
293        #print quantity.vertex_values
294        assert allclose(quantity.vertex_values[1,0], 2.0)
295        assert allclose(quantity.vertex_values[1,1], 6.0)       
296        assert allclose(quantity.vertex_values[1,2], 8.0)
297
298
299
300    def test_first_order_extrapolator(self):
301        quantity = Conserved_quantity(self.mesh4)               
302
303        #Test centroids
304        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
305        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
306
307        #Extrapolate
308        quantity.extrapolate_first_order()
309
310        #Check vertices but not edge values
311        assert allclose(quantity.vertex_values,
312                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
313
314
315    def test_second_order_extrapolator(self):
316        quantity = Conserved_quantity(self.mesh4)                       
317
318        #Set up for a gradient of (3,0) at mid triangle
319        quantity.set_values([2.0, 4.0, 8.0, 2.0],
320                            location = 'centroids')
321       
322
323
324        quantity.extrapolate_second_order()
325        quantity.limit()
326
327
328        #Assert that central triangle is limited by neighbours
329        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
330        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
331       
332        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
333        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
334       
335        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
336        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
337
338       
339        #Assert that quantities are conserved
340        from Numeric import sum
341        for k in range(quantity.centroid_values.shape[0]):
342            assert allclose (quantity.centroid_values[k],
343                             sum(quantity.vertex_values[k,:])/3)
344       
345       
346
347       
348
349    def test_limiter(self):
350        quantity = Conserved_quantity(self.mesh4)
351
352        #Create a deliberate overshoot (e.g. from gradient computation)
353        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
354
355
356        #Limit
357        quantity.limit()
358
359        #Assert that central triangle is limited by neighbours
360        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
361        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
362       
363        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
364        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
365       
366        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
367        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
368
369
370       
371        #Assert that quantities are conserved
372        from Numeric import sum
373        for k in range(quantity.centroid_values.shape[0]):
374            assert allclose (quantity.centroid_values[k],
375                             sum(quantity.vertex_values[k,:])/3)
376       
377
378    def test_limiter2(self):
379        """Taken from test_shallow_water
380        """
381        quantity = Conserved_quantity(self.mesh4)
382
383        #Test centroids
384        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
385        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
386
387
388        #Extrapolate
389        quantity.extrapolate_second_order()
390
391        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
392
393        #Limit
394        quantity.limit()
395
396       
397        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
398       
399
400        #Assert that quantities are conserved
401        from Numeric import sum
402        for k in range(quantity.centroid_values.shape[0]):
403            assert allclose (quantity.centroid_values[k],
404                             sum(quantity.vertex_values[k,:])/3)
405
406
407       
408
409
410    def test_distribute_first_order(self):
411        quantity = Conserved_quantity(self.mesh4)       
412
413        #Test centroids
414        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
415        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
416
417
418        #Extrapolate
419        quantity.extrapolate_first_order()
420
421        #Interpolate
422        quantity.interpolate_from_vertices_to_edges()       
423
424        assert allclose(quantity.vertex_values,
425                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
426        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
427                                               [3,3,3], [4, 4, 4]])
428
429
430    def test_distribute_second_order(self):
431        quantity = Conserved_quantity(self.mesh4)       
432
433        #Test centroids
434        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
435        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
436
437
438        #Extrapolate
439        quantity.extrapolate_second_order()
440
441        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
442
443
444    def test_update_explicit(self):
445        quantity = Conserved_quantity(self.mesh4)
446
447        #Test centroids
448        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
449        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
450
451        #Set explicit_update
452        quantity.explicit_update = array( [1.,1.,1.,1.] )
453
454        #Update with given timestep
455        quantity.update(0.1)
456
457        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
458        assert allclose( quantity.centroid_values, x)
459
460    def test_update_semi_implicit(self):
461        quantity = Conserved_quantity(self.mesh4)
462
463        #Test centroids
464        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
465        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
466
467        #Set semi implicit update
468        quantity.semi_implicit_update = array([1.,1.,1.,1.])
469
470        #Update with given timestep
471        timestep = 0.1
472        quantity.update(timestep)
473
474        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
475        denom = ones(4, Float)-timestep*sem
476       
477        x = array([1, 2, 3, 4])/denom
478        assert allclose( quantity.centroid_values, x)
479
480
481    def test_both_updates(self):
482        quantity = Conserved_quantity(self.mesh4)
483
484        #Test centroids
485        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
486        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
487
488        #Set explicit_update
489        quantity.explicit_update = array( [4.,3.,2.,1.] )
490       
491        #Set semi implicit update
492        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
493
494        #Update with given timestep
495        timestep = 0.1       
496        quantity.update(0.1)
497
498        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
499        denom = ones(4, Float)-timestep*sem
500
501        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
502        x /= denom
503        assert allclose( quantity.centroid_values, x)       
504
505       
506
507
508    #Test smoothing   
509    def test_smoothing(self):
510
511        from mesh_factory import rectangular
512        from shallow_water import Domain, Transmissive_boundary
513        from Numeric import zeros, Float
514        from util import mean
515
516        #Create basic mesh
517        points, vertices, boundary = rectangular(2, 2)
518
519        #Create shallow water domain
520        domain = Domain(points, vertices, boundary)
521        domain.default_order=2
522        domain.reduction = mean
523
524       
525        #Set some field values
526        domain.set_quantity('elevation', lambda x,y: x)       
527        domain.set_quantity('friction', 0.03)
528
529
530        ######################
531        # Boundary conditions
532        B = Transmissive_boundary(domain)
533        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
534       
535
536        ######################
537        #Initial condition - with jumps
538
539        bed = domain.quantities['elevation'].vertex_values
540        stage = zeros(bed.shape, Float)
541
542        h = 0.03
543        for i in range(stage.shape[0]):
544            if i % 2 == 0:           
545                stage[i,:] = bed[i,:] + h
546            else:
547                stage[i,:] = bed[i,:]
548               
549        domain.set_quantity('stage', stage)
550
551        stage = domain.quantities['stage']
552
553        #Get smoothed stage
554        A, V = stage.get_vertex_values(xy=False, smooth=True)
555        Q = stage.vertex_values
556
557       
558        assert A.shape[0] == 9 
559        assert V.shape[0] == 8 
560        assert V.shape[1] == 3                 
561       
562        #First four points
563        assert allclose(A[0], (Q[0,2] + Q[1,1])/2) 
564        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
565        assert allclose(A[2], Q[3,0])
566        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
567
568        #Center point
569        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
570                               Q[5,0] + Q[6,2] + Q[7,1])/6)
571                                 
572
573        #Check V
574        assert allclose(V[0,:], [3,4,0])
575        assert allclose(V[1,:], [1,0,4])
576        assert allclose(V[2,:], [4,5,1])                       
577        assert allclose(V[3,:], [2,1,5])
578        assert allclose(V[4,:], [6,7,3])       
579        assert allclose(V[5,:], [4,3,7])
580        assert allclose(V[6,:], [7,8,4])
581        assert allclose(V[7,:], [5,4,8])                       
582
583        #Get smoothed stage with XY
584        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
585
586        assert allclose(A, A1)
587        assert allclose(V, V1)       
588
589        #Check XY
590        assert allclose(X[4], 0.5)
591        assert allclose(Y[4], 0.5)
592
593        assert allclose(X[7], 1.0)
594        assert allclose(Y[7], 0.5)                 
595
596
597
598
599    def test_vertex_values_no_smoothing(self):
600
601        from mesh_factory import rectangular
602        from shallow_water import Domain, Transmissive_boundary
603        from Numeric import zeros, Float
604        from util import mean
605
606       
607        #Create basic mesh
608        points, vertices, boundary = rectangular(2, 2)
609
610        #Create shallow water domain
611        domain = Domain(points, vertices, boundary)
612        domain.default_order=2
613        domain.reduction = mean
614
615       
616        #Set some field values
617        domain.set_quantity('elevation', lambda x,y: x)       
618        domain.set_quantity('friction', 0.03)
619
620
621        ######################
622        #Initial condition - with jumps
623
624        bed = domain.quantities['elevation'].vertex_values
625        stage = zeros(bed.shape, Float)
626
627        h = 0.03
628        for i in range(stage.shape[0]):
629            if i % 2 == 0:           
630                stage[i,:] = bed[i,:] + h
631            else:
632                stage[i,:] = bed[i,:]
633               
634        domain.set_quantity('stage', stage)
635
636        #Get stage
637        stage = domain.quantities['stage']       
638        A, V = stage.get_vertex_values(xy=False, smooth=False)
639        Q = stage.vertex_values.flat
640
641        for k in range(8):
642            assert allclose(A[k], Q[k])
643
644           
645        for k in range(8):
646            assert V[k, 0] == 3*k
647            assert V[k, 1] == 3*k+1
648            assert V[k, 2] == 3*k+2           
649           
650
651
652        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
653
654       
655        assert allclose(A, A1)
656        assert allclose(V, V1)       
657
658        #Check XY
659        assert allclose(X[1], 0.5)
660        assert allclose(Y[1], 0.5)
661        assert allclose(X[4], 0.0)
662        assert allclose(Y[4], 0.0)
663        assert allclose(X[12], 1.0)
664        assert allclose(Y[12], 0.0)               
665       
666       
667 
668    def set_array_values_by_index(self):
669
670        from mesh_factory import rectangular
671        from shallow_water import Domain
672        from Numeric import zeros, Float
673       
674        #Create basic mesh
675        points, vertices, boundary = rectangular(1, 1)
676
677        #Create shallow water domain
678        domain = Domain(points, vertices, boundary)
679        #print "domain.number_of_elements ",domain.number_of_elements
680        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
681        value = [7]
682        indexes = [1]
683        quantity.set_array_values_by_index(value,
684                                           location = 'centroids',
685                                           indexes = indexes)
686        #print "quantity.centroid_values",quantity.centroid_values
687       
688        assert allclose(quantity.centroid_values, [1,7])
689
690        quantity.set_array_values([15,20,25], indexes = indexes)
691        assert allclose(quantity.centroid_values, [1,20])
692
693        quantity.set_array_values([15,20,25], indexes = indexes)
694        assert allclose(quantity.centroid_values, [1,20])
695       
696    def test_setting_some_vertex_values(self):
697        """
698        set values based on triangle lists.
699        """
700        from mesh_factory import rectangular
701        from shallow_water import Domain
702        from Numeric import zeros, Float
703       
704        #Create basic mesh
705        points, vertices, boundary = rectangular(1, 3)
706        #print "vertices",vertices
707        #Create shallow water domain
708        domain = Domain(points, vertices, boundary)
709        #print "domain.number_of_elements ",domain.number_of_elements
710        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
711                                    [4,4,4],[5,5,5],[6,6,6]])
712        value = [7]
713        indexes = [1]
714        quantity.set_values(value,
715                                  location = 'centroids',
716                                  indexes = indexes)
717        #print "quantity.centroid_values",quantity.centroid_values
718        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
719       
720        value = [[15,20,25]]
721        quantity.set_values(value, indexes = indexes)
722        #print "1 quantity.vertex_values",quantity.vertex_values
723        assert allclose(quantity.vertex_values[1], value[0])
724
725
726        #print "quantity",quantity.vertex_values
727        values = [10,100,50]
728        quantity.set_values(values, indexes = [0,1,5], location = 'centroids')
729        #print "2 quantity.vertex_values",quantity.vertex_values
730        assert allclose(quantity.vertex_values[0], [10,10,10])
731        assert allclose(quantity.vertex_values[5], [50,50,50])
732        #quantity.interpolate()
733        #print "quantity.centroid_values",quantity.centroid_values
734        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
735
736       
737        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
738                                    [4,4,4],[5,5,5],[6,6,6]])
739        values = [10,100,50]
740        #this will be per unique vertex, indexing the vertices
741        #print "quantity.vertex_values",quantity.vertex_values
742        quantity.set_values(values, indexes = [0,1,5]) 
743        #print "quantity.vertex_values",quantity.vertex_values
744        assert allclose(quantity.vertex_values[0], [1,50,10])
745        assert allclose(quantity.vertex_values[5], [6,6,6])
746        assert allclose(quantity.vertex_values[1], [100,10,50])
747
748        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
749                                    [4,4,4],[5,5,5],[6,6,6]])
750        values = [[31,30,29],[400,400,400],[1000,999,998]]
751        quantity.set_values(values, indexes = [3,3,5])
752        quantity.interpolate()
753        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
754
755        values = [[1,1,1],[2,2,2],[3,3,3],
756                                    [4,4,4],[5,5,5],[6,6,6]]
757        quantity.set_values(values)
758       
759        # testing the standard set values by vertex
760        # indexed by vertex_id in general_mesh.coordinates
761        values = [0,1,2,3,4,5,6,7]
762       
763        quantity.set_values(values)
764        #print "1 quantity.vertex_values",quantity.vertex_values
765        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
766                                                [ 1.,  0.,  5.],
767                                                [ 5.,  6.,  1.],
768                                                [ 2.,  1.,  6.],
769                                                [ 6.,  7.,  2.],
770                                                [ 3.,  2.,  7.]])
771   
772    def test_setting_unique_vertex_values(self):
773        """
774        set values based on unique_vertex lists.
775        """
776        from mesh_factory import rectangular
777        from shallow_water import Domain
778        from Numeric import zeros, Float
779       
780        #Create basic mesh
781        points, vertices, boundary = rectangular(1, 3)
782        #print "vertices",vertices
783        #Create shallow water domain
784        domain = Domain(points, vertices, boundary)
785        #print "domain.number_of_elements ",domain.number_of_elements
786        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
787                                    [4,4,4],[5,5,5]])
788        value = 7
789        indexes = [1,5]
790        quantity.set_values(value,
791                                  location = 'unique vertices',
792                                  indexes = indexes)
793        #print "quantity.centroid_values",quantity.centroid_values
794        assert allclose(quantity.vertex_values[0], [0,7,0])
795        assert allclose(quantity.vertex_values[1], [7,1,7])
796        assert allclose(quantity.vertex_values[2], [7,2,7])
797       
798   
799    def test_get_values(self):
800        """
801        get values based on triangle lists.
802        """
803        from mesh_factory import rectangular
804        from shallow_water import Domain
805        from Numeric import zeros, Float
806       
807        #Create basic mesh
808        points, vertices, boundary = rectangular(1, 3)
809
810        #print "points",points
811        #print "vertices",vertices
812        #print "boundary",boundary
813
814        #Create shallow water domain
815        domain = Domain(points, vertices, boundary)
816        #print "domain.number_of_elements ",domain.number_of_elements
817        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
818                                    [4,4,4],[5,5,5]])
819       
820        #print "quantity.get_values(location = 'unique vertices')", \
821        #      quantity.get_values(location = 'unique vertices')
822
823        #print "quantity.get_values(location = 'unique vertices')", \
824        #      quantity.get_values(indexes=[0,1,2,3,4,5,6,7], \
825        #                          location = 'unique vertices')
826       
827        answer = [0.5,2,4,5,0,1,3,4.5]
828        assert allclose(answer,
829                        quantity.get_values(location = 'unique vertices'))
830       
831        indexes = [0,5,3]
832        answer = [0.5,1,5]
833        assert allclose(answer,
834                        quantity.get_values(indexes=indexes, \
835                                            location = 'unique vertices'))
836        #print "quantity.centroid_values",quantity.centroid_values
837        #print "quantity.get_values(location = 'centroids') ",\
838        #      quantity.get_values(location = 'centroids')
839       
840    def test_getting_some_vertex_values(self):
841        """
842        get values based on triangle lists.
843        """
844        from mesh_factory import rectangular
845        from shallow_water import Domain
846        from Numeric import zeros, Float
847       
848        #Create basic mesh
849        points, vertices, boundary = rectangular(1, 3)
850
851        #print "points",points
852        #print "vertices",vertices
853        #print "boundary",boundary
854
855        #Create shallow water domain
856        domain = Domain(points, vertices, boundary)
857        #print "domain.number_of_elements ",domain.number_of_elements
858        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
859                                    [4,4,4],[5,5,5],[6,6,6]])
860        value = [7]
861        indexes = [1]
862        quantity.set_values(value,
863                                  location = 'centroids',
864                                  indexes = indexes)
865        #print "quantity.centroid_values",quantity.centroid_values
866        #print "quantity.get_values(location = 'centroids') ",\
867        #      quantity.get_values(location = 'centroids')
868        assert allclose(quantity.centroid_values,
869                        quantity.get_values(location = 'centroids'))
870
871       
872        value = [[15,20,25]]
873        quantity.set_values(value, indexes = indexes)
874        #print "1 quantity.vertex_values",quantity.vertex_values
875        assert allclose(quantity.vertex_values, quantity.get_values())
876
877        assert allclose(quantity.edge_values,
878                        quantity.get_values(location = 'edges'))
879
880        # get a subset of elements
881        subset = quantity.get_values(location='centroids', indexes=[0,5])
882        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
883        assert allclose(subset, answer)
884
885
886        subset = quantity.get_values(location='edges', indexes=[0,5])
887        answer = [quantity.edge_values[0],quantity.edge_values[5]]
888        #print "subset",subset
889        #print "answer",answer
890        assert allclose(subset, answer)
891
892        subset = quantity.get_values( indexes=[1,5])
893        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
894        #print "subset",subset
895        #print "answer",answer
896        assert allclose(subset, answer)
897
898       
899       
900#-------------------------------------------------------------
901if __name__ == "__main__":
902    suite = unittest.makeSuite(TestCase,'test_limiter2')
903    #print "restricted test"
904    #suite = unittest.makeSuite(TestCase,'test_set_vertex_values_subset')
905    runner = unittest.TextTestRunner()
906    runner.run(suite)
907
908
Note: See TracBrowser for help on using the repository browser.