source: anuga_work/development/anuga_1d/test_quantity.py @ 7793

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

Setting up test function for Quantity

File size: 29.7 KB
Line 
1import quantity
2#!/usr/bin/env python
3
4import unittest
5from math import sqrt, pi
6
7
8from quantity import *
9from domain import *
10
11
12from Numeric import allclose, array, ones, Float, maximum, zeros
13
14
15class Test_Quantity(unittest.TestCase):
16    def setUp(self):
17        self.points3        = [0.0, 1.0, 2.0, 3.0]
18        self.vertex_values3 = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
19        self.domain3        = Domain(self.points3)
20
21
22       
23        self.points4          = [0.0, 1.0, 2.5, 3.0, 5.0]
24        self.vertex_values4   = [[1.0,2.0],[4.0,5.0],[-1.0,2.0],[3.0,6.0]]
25        self.centroid_values4 = [1.5, 4.5, 0.5, 4.5]
26        self.boundary4        = {(0, 0): 'left', (3, 1): 'right'}
27        self.domain4          = Domain(self.points4,self.boundary4)
28
29        self.points10 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
30        self.domain10 = Domain(self.points10)
31
32        self.points6 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
33        self.domain6 = Domain(self.points6)
34
35
36    def tearDown(self):
37        pass
38        #print "  Tearing down"
39
40
41    def test_creat_with_boundary(self):
42
43        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
44       
45    def test_creation(self):
46
47        quantity = Quantity(self.domain3)
48        assert allclose(quantity.vertex_values, [[0.0,0.0],[0.0,0.0],[0.0,0.0]])
49
50       
51        try:
52            quantity = Quantity()
53        except:
54            pass
55        else:
56            raise 'Should have raised empty quantity exception'
57
58
59        try:
60            quantity = Quantity([1,2,3])
61        except AssertionError:
62            pass
63        except:
64            raise 'Should have raised "mising domain object" error'
65
66
67    def test_creation_zeros(self):
68
69        quantity = Quantity(self.domain3)
70        assert allclose(quantity.centroid_values, [[0.,0.,0.]])
71
72
73        quantity = Quantity(self.domain4)
74        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
75                                                 [0.,0.], [0.,0.]])
76
77
78    def test_interpolation(self):
79        quantity = Quantity(self.domain4, self.vertex_values4)
80        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
81
82
83
84    def test_interpolation2(self):
85        quantity = Quantity(self.domain4, self.vertex_values4)
86        assert allclose(quantity.centroid_values, self.centroid_values4) #Centroid
87
88        quantity.extrapolate_second_order()
89
90        #print quantity.vertex_values
91        assert allclose(quantity.vertex_values,[[ 0.3,  2.7],
92                                                [ 4.5,  4.5],
93                                                [ 0.5,  0.5],
94                                                [ 1.3,  7.7]])
95
96 
97
98 
99    def test_boundary_allocation(self):
100        quantity = Quantity(self.domain4,
101                            [[1,2], [5,5], [0,9], [-6, 3]])
102
103
104        assert quantity.boundary_values.shape[0] == len(self.domain4.boundary)
105
106
107    def test_set_values(self):
108        quantity = Quantity(self.domain4)
109
110
111        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
112                            location = 'vertices')
113        assert allclose(quantity.vertex_values,
114                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
115        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
116        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
117                                               [5., 5., 5.],
118                                               [4.5, 4.5, 0.],
119                                               [3.0, -1.5, -1.5]])
120
121
122        # Test default
123        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
124        assert allclose(quantity.vertex_values,
125                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
126        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
127        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
128                                               [5., 5., 5.],
129                                               [4.5, 4.5, 0.],
130                                               [3.0, -1.5, -1.5]])
131
132        # Test centroids
133        quantity.set_values([1,2,3,4], location = 'centroids')
134        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
135
136        # Test exceptions
137        try:
138            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
139                                location = 'bas kamel tuba')
140        except:
141            pass
142
143
144        try:
145            quantity.set_values([[1,2,3], [0,0,9]])
146        except AssertionError:
147            pass
148        except:
149            raise 'should have raised Assertionerror'
150
151
152
153    def test_set_values_const(self):
154        quantity = Quantity(self.domain4)
155
156        quantity.set_values(1.0, location = 'vertices')
157        assert allclose(quantity.vertex_values,
158                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
159
160        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
161        assert allclose(quantity.edge_values, [[1, 1, 1],
162                                               [1, 1, 1],
163                                               [1, 1, 1],
164                                               [1, 1, 1]])
165
166
167        quantity.set_values(2.0, location = 'centroids')
168        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
169
170
171    def test_set_values_func(self):
172        quantity = Quantity(self.domain4)
173
174        def f(x, y):
175            return x+y
176
177        quantity.set_values(f, location = 'vertices')
178        #print "quantity.vertex_values",quantity.vertex_values
179        assert allclose(quantity.vertex_values,
180                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
181        assert allclose(quantity.centroid_values,
182                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
183        assert allclose(quantity.edge_values,
184                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
185
186
187        quantity.set_values(f, location = 'centroids')
188        assert allclose(quantity.centroid_values,
189                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
190
191
192    def test_integral(self):
193        quantity = Quantity(self.domain4)
194
195        #Try constants first
196        const = 5
197        quantity.set_values(const, location = 'vertices')
198        #print 'Q', quantity.get_integral()
199
200        assert allclose(quantity.get_integral(), self.domain4.get_area() * const)
201
202        #Try with a linear function
203        def f(x):
204            return x
205
206        quantity.set_values(f, location = 'vertices')
207
208 
209        assert allclose (quantity.centroid_values,
210        [ 0.5,   1.75,  2.75,  4.  ])
211
212        assert allclose (quantity.vertex_values, [[ 0.,   1. ],
213                                                  [ 1.,   2.5],
214                                                  [ 2.5,  3. ],
215                                                  [ 3.,   5. ]])
216
217
218        ref_integral = 0.5 + 1.5*1.75 + 0.5*2.75 + 2.0*4.0
219
220        assert allclose (quantity.get_integral(), ref_integral)
221
222
223
224    def test_set_vertex_values(self):
225        quantity = Quantity(self.domain4)
226        quantity.set_vertex_values([0,1,2,3,4,5])
227
228        assert allclose(quantity.vertex_values,
229                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
230        assert allclose(quantity.centroid_values,
231                        [1., 7./3, 11./3, 8./3]) #Centroid
232        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
233                                               [3., 2.5, 1.5],
234                                               [3.5, 4.5, 3.],
235                                               [2.5, 3.5, 2]])
236
237
238    def test_set_vertex_values_subset(self):
239        quantity = Quantity(self.domain4)
240        quantity.set_vertex_values([0,1,2,3,4,5])
241        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
242
243        assert allclose(quantity.vertex_values,
244                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
245
246
247    def test_set_vertex_values_using_general_interface(self):
248        quantity = Quantity(self.domain4)
249
250
251        quantity.set_values([0,1,2,3,4,5])
252
253
254        assert allclose(quantity.vertex_values,
255                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
256
257        #Centroid
258        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
259
260        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
261                                               [3., 2.5, 1.5],
262                                               [3.5, 4.5, 3.],
263                                               [2.5, 3.5, 2]])
264
265
266
267    def test_set_vertex_values_using_general_interface_with_subset(self):
268        """test_set_vertex_values_using_general_interface_with_subset(self):
269       
270        Test that indices and polygon works (for constants values)
271        """
272       
273        quantity = Quantity(self.domain4)
274
275
276        quantity.set_values([0,2,3,5], indices=[0,2,3,5])
277        assert allclose(quantity.vertex_values,
278                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
279
280
281        # Constant
282        quantity.set_values(0.0)
283        quantity.set_values(3.14, indices=[0,2], location='vertices')
284
285        # Indices refer to triangle numbers
286        assert allclose(quantity.vertex_values,
287                        [[3.14,3.14,3.14], [0,0,0],
288                         [3.14,3.14,3.14], [0,0,0]])       
289       
290
291
292        # Now try with polygon (pick points where y>2)
293        polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
294        quantity.set_values(0.0)
295        quantity.set_values(3.14, polygon=polygon)
296       
297        assert allclose(quantity.vertex_values,
298                        [[0,0,0], [0,0,0], [0,0,0],
299                         [3.14,3.14,3.14]])               
300
301
302        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
303        # using centroids
304        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
305        quantity.set_values(0.0)
306        quantity.set_values(3.14, location='centroids', polygon=polygon)
307        assert allclose(quantity.vertex_values,
308                        [[0,0,0],
309                         [3.14,3.14,3.14],
310                         [3.14,3.14,3.14],                         
311                         [0,0,0]])               
312
313
314        # Same polygon now use vertices (default)
315        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
316        quantity.set_values(0.0)
317        #print 'Here 2'
318        quantity.set_values(3.14, polygon=polygon)
319        assert allclose(quantity.vertex_values,
320                        [[0,0,0],
321                         [3.14,3.14,3.14],
322                         [3.14,3.14,3.14],                         
323                         [0,0,0]])               
324       
325
326        # Test input checking
327        try:
328            quantity.set_values(3.14, polygon=polygon, indices = [0,2])
329        except:
330            pass
331        else:
332            msg = 'Should have caught this'
333            raise msg
334
335
336
337
338
339    def test_set_vertex_values_using_general_interface_subset_and_geo(self):
340        """test_set_vertex_values_using_general_interface_with_subset(self):
341        Test that indices and polygon works using georeferencing
342        """
343       
344        quantity = Quantity(self.domain4)
345        G = Geo_reference(56, 10, 100)
346        quantity.domain.geo_reference = G
347
348        #print quantity.domain.get_nodes(absolute=True)
349
350
351        # Constant
352        quantity.set_values(0.0)
353        quantity.set_values(3.14, indices=[0,2], location='vertices')
354
355        # Indices refer to triangle numbers here - not vertices (why?)
356        assert allclose(quantity.vertex_values,
357                        [[3.14,3.14,3.14], [0,0,0],
358                         [3.14,3.14,3.14], [0,0,0]])       
359       
360
361
362        # Now try with polygon (pick points where y>2)
363        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
364        polygon += [G.xllcorner, G.yllcorner]
365       
366        quantity.set_values(0.0)
367        quantity.set_values(3.14, polygon=polygon, location='centroids')
368       
369        assert allclose(quantity.vertex_values,
370                        [[0,0,0], [0,0,0], [0,0,0],
371                         [3.14,3.14,3.14]])               
372
373
374        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
375        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
376        polygon += [G.xllcorner, G.yllcorner]
377       
378        quantity.set_values(0.0)
379        quantity.set_values(3.14, polygon=polygon)
380
381        assert allclose(quantity.vertex_values,
382                        [[0,0,0],
383                         [3.14,3.14,3.14],
384                         [3.14,3.14,3.14],                         
385                         [0,0,0]])               
386
387
388
389
390    def test_set_values_from_quantity(self):
391
392        quantity1 = Quantity(self.domain4)
393        quantity1.set_vertex_values([0,1,2,3,4,5])
394
395        assert allclose(quantity1.vertex_values,
396                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
397
398
399        quantity2 = Quantity(self.domain4)
400        quantity2.set_values(quantity=quantity1)
401        assert allclose(quantity2.vertex_values,
402                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
403
404        quantity2.set_values(quantity = 2*quantity1)
405        assert allclose(quantity2.vertex_values,
406                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
407
408        quantity2.set_values(quantity = 2*quantity1 + 3)
409        assert allclose(quantity2.vertex_values,
410                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
411
412
413        #Check detection of quantity as first orgument
414        quantity2.set_values(2*quantity1 + 3)
415        assert allclose(quantity2.vertex_values,
416                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
417
418
419
420
421
422    def test_overloading(self):
423
424        quantity1 = Quantity(self.domain4)
425        quantity1.set_values( [[0,1],[1,2],[2,3],[3,4]],
426                            location = 'vertices')
427
428        assert allclose(quantity1.vertex_values,
429                        [[0,1], [1,2], [2,3], [3,4]])
430
431
432        quantity2 = Quantity(self.domain4)
433        quantity2.set_values([[1,2], [5,5], [0,9], [-6, 3]],
434                             location = 'vertices')
435
436
437
438        quantity3 = Quantity(self.domain4)
439        quantity3.set_values([[2,2], [7,8], [7,6], [3, -8]],
440                             location = 'vertices')
441
442
443        # Negation
444        Q = -quantity1
445        assert allclose(Q.vertex_values, -quantity1.vertex_values)
446        assert allclose(Q.centroid_values, -quantity1.centroid_values)
447
448
449        # Addition
450        Q = quantity1 + 7
451        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
452        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
453
454        Q = 7 + quantity1
455        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
456        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
457
458        Q = quantity1 + quantity2
459        assert allclose(Q.vertex_values,
460                        quantity1.vertex_values + quantity2.vertex_values)
461        assert allclose(Q.centroid_values,
462                        quantity1.centroid_values + quantity2.centroid_values)
463
464        Q = quantity1 + quantity2 - 3
465        assert allclose(Q.vertex_values,
466                        quantity1.vertex_values + quantity2.vertex_values - 3)
467
468        Q = quantity1 - quantity2
469        assert allclose(Q.vertex_values,
470                        quantity1.vertex_values - quantity2.vertex_values)
471
472        #Scaling
473        Q = quantity1*3
474        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
475        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
476
477        Q = 3*quantity1
478        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
479
480        #Multiplication
481        Q = quantity1 * quantity2
482        assert allclose(Q.vertex_values,
483                        quantity1.vertex_values * quantity2.vertex_values)
484
485        #Linear combinations
486        Q = 4*quantity1 + 2
487        assert allclose(Q.vertex_values,
488                        4*quantity1.vertex_values + 2)
489
490        Q = quantity1*quantity2 + 2
491        assert allclose(Q.vertex_values,
492                        quantity1.vertex_values * quantity2.vertex_values + 2)
493
494        Q = quantity1*quantity2 + quantity3
495        assert allclose(Q.vertex_values,
496                        quantity1.vertex_values * quantity2.vertex_values +
497                        quantity3.vertex_values)
498        Q = quantity1*quantity2 + 3*quantity3
499        assert allclose(Q.vertex_values,
500                        quantity1.vertex_values * quantity2.vertex_values +
501                        3*quantity3.vertex_values)
502        Q = quantity1*quantity2 + 3*quantity3 + 5.0
503        assert allclose(Q.vertex_values,
504                        quantity1.vertex_values * quantity2.vertex_values +
505                        3*quantity3.vertex_values + 5)
506
507        Q = quantity1*quantity2 - quantity3
508        assert allclose(Q.vertex_values,
509                        quantity1.vertex_values * quantity2.vertex_values -
510                        quantity3.vertex_values)
511        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
512        assert allclose(Q.vertex_values,
513                        1.5*quantity1.vertex_values * quantity2.vertex_values -
514                        3*quantity3.vertex_values + 5)
515
516        #Try combining quantities and arrays and scalars
517        Q = 1.5*quantity1*quantity2.vertex_values -\
518            3*quantity3.vertex_values + 5.0
519        assert allclose(Q.vertex_values,
520                        1.5*quantity1.vertex_values * quantity2.vertex_values -
521                        3*quantity3.vertex_values + 5)
522
523
524        #Powers
525        Q = quantity1**2
526        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
527
528        Q = quantity1**2 +quantity2**2
529        assert allclose(Q.vertex_values,
530                        quantity1.vertex_values**2 + \
531                        quantity2.vertex_values**2)
532
533        Q = (quantity1**2 +quantity2**2)**0.5
534        assert allclose(Q.vertex_values,
535                        (quantity1.vertex_values**2 + \
536                        quantity2.vertex_values**2)**0.5)
537
538    def test_compute_gradient(self):
539        quantity = Quantity(self.domain6)
540
541        #Set up for a gradient of (2,0) at mid triangle
542        quantity.set_values([2.0, 4.0, 4.0, 5.0, 10.0, 12.0],
543                            location = 'centroids')
544
545
546        #Gradients
547        quantity.compute_gradients()
548
549        a = quantity.gradients
550
551        assert allclose(a, [ 3., 1., 0.5, 3., 3.5, 0.5])
552       
553        quantity.extrapolate_second_order()
554
555
556        assert allclose(quantity.vertex_values, [[ 1., 3. ],
557                                                 [ 4., 4. ],
558                                                 [ 4., 4. ],
559                                                 [ 4., 6.],
560                                                 [ 8.25, 11.75],
561                                                 [ 11.,  13. ]])
562
563                                                 
564
565    def test_second_order_extrapolation2(self):
566        quantity = Quantity(self.domain4)
567
568        #Set up for a gradient of (3,1), f(x) = 3x+y
569        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
570                            location = 'centroids')
571
572        #Gradients
573        quantity.compute_gradients()
574
575        a = quantity.x_gradient
576        b = quantity.y_gradient
577       
578        #print a, b
579
580        assert allclose(a[1], 3.0)
581        assert allclose(b[1], 1.0)
582
583        #Work out the others
584
585        quantity.extrapolate_second_order()
586
587        #print quantity.vertex_values
588        assert allclose(quantity.vertex_values[1,0], 2.0)
589        assert allclose(quantity.vertex_values[1,1], 6.0)
590        assert allclose(quantity.vertex_values[1,2], 8.0)
591
592
593
594    def test_backup_saxpy_centroid_values(self):
595        quantity = Quantity(self.domain4)
596
597        #Set up for a gradient of (3,1), f(x) = 3x+y
598        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
599        d_values = array([1.0, 2.0, 3.0, 4.0])
600        quantity.set_values(c_values, location = 'centroids')
601
602        #Backup
603        quantity.backup_centroid_values()
604
605        #print quantity.vertex_values
606        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
607
608
609        quantity.set_values(d_values, location = 'centroids')
610
611        quantity.saxpy_centroid_values(2.0, 3.0)
612
613        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
614
615
616
617    def test_first_order_extrapolator(self):
618        quantity = Quantity(self.domain4)
619
620        centroid_values = array([1.,2.,3.,4.])
621        quantity.set_values(centroid_values, location = 'centroids')
622        assert allclose(quantity.centroid_values, centroid_values) #Centroid
623
624        #Extrapolate
625        quantity.extrapolate_first_order()
626
627        #Check that gradient is zero
628        a = quantity.gradients
629        assert allclose(a, [0,0,0,0])
630       
631
632        #Check vertices but not edge values
633        assert allclose(quantity.vertex_values,
634                        [[1,1], [2,2], [3,3], [4,4]])
635
636
637    def test_second_order_extrapolator(self):
638        quantity = Quantity(self.domain4)
639
640        #Set up for a gradient of (3,0) at mid triangle
641        quantity.set_values([2.0, 4.0, 8.0, 2.0],
642                            location = 'centroids')
643
644
645
646        quantity.extrapolate_second_order()
647        quantity.limit()
648
649
650        #Assert that central triangle is limited by neighbours
651        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
652        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
653
654        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
655        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
656
657        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
658        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
659
660
661        #Assert that quantities are conserved
662        from Numeric import sum
663        for k in range(quantity.centroid_values.shape[0]):
664            assert allclose (quantity.centroid_values[k],
665                             sum(quantity.vertex_values[k,:])/3)
666
667
668    def test_limit(self):
669        quantity = Quantity(self.domain4)
670
671        #Create a deliberate overshoot (e.g. from gradient computation)
672        quantity.set_values([[0,0], [2,20], [-20,3], [8,3]])
673
674        #Limit
675        quantity.limit_minmod()
676
677       
678        #cells 1 and 2 are local max and min
679        assert quantity.vertex_values[1][0] == quantity.centroid_values[1]
680        assert quantity.vertex_values[1][1] == quantity.centroid_values[1]
681
682        assert quantity.vertex_values[2][0] == quantity.centroid_values[2]
683        assert quantity.vertex_values[2][1] == quantity.centroid_values[2]
684
685
686
687    def test_distribute_first_order(self):
688        quantity = Quantity(self.domain4)
689
690        #Test centroids
691        centroid_values = array([1.,2.,3.,4.])
692        quantity.set_values(centroid_values, location = 'centroids')
693        assert allclose(quantity.centroid_values, centroid_values) #Centroid
694
695
696        #Extrapolate from centroid to vertices and edges
697        quantity.extrapolate_first_order()
698       
699        assert allclose(quantity.vertex_values,[[ 1.,  1.],
700                                                [ 2.,  2.],
701                                                [ 3.,  3.],
702                                                [ 4.,  4.]])
703 
704
705
706    def test_distribute_second_order(self):
707        quantity = Quantity(self.domain4)
708
709        #Test centroids
710        centroid_values = array([2.,4.,8.,2.])
711        quantity.set_values(centroid_values, location = 'centroids')
712        assert allclose(quantity.centroid_values, centroid_values) #Centroid
713
714
715        #Extrapolate
716        quantity.extrapolate_second_order()
717
718        assert allclose(quantity.vertex_values, [[ 1.2,  2.8],
719                                                 [ 2.,  6. ],
720                                                 [ 8.,   8. ],
721                                                 [ 6.8, -2.8]])
722
723
724    def test_update_explicit(self):
725        quantity = Quantity(self.domain4)
726
727        #Test centroids
728        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
729        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
730
731        #Set explicit_update
732        quantity.explicit_update = array( [1.,1.,1.,1.] )
733
734        #Update with given timestep
735        quantity.update(0.1)
736
737        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
738        assert allclose( quantity.centroid_values, x)
739
740    def test_update_semi_implicit(self):
741        quantity = Quantity(self.domain4)
742
743        #Test centroids
744        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
745        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
746
747        #Set semi implicit update
748        quantity.semi_implicit_update = array([1.,1.,1.,1.])
749
750        #Update with given timestep
751        timestep = 0.1
752        quantity.update(timestep)
753
754        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
755        denom = ones(4, Float)-timestep*sem
756
757        x = array([1, 2, 3, 4])/denom
758        assert allclose( quantity.centroid_values, x)
759
760
761    def test_both_updates(self):
762        quantity = Quantity(self.domain4)
763
764        #Test centroids
765        centroid_values = array( [1, 2, 3, 4] )
766        quantity.set_values(centroid_values, location = 'centroids')
767        assert allclose(quantity.centroid_values, centroid_values) #Centroid
768
769        #Set explicit_update
770        explicit_update = array( [4.,3.,2.,1.] )
771        quantity.explicit_update[:] = explicit_update
772
773        #Set semi implicit update
774        semi_implicit_update = array( [1.,1.,1.,1.] )
775        quantity.semi_implicit_update[:] = semi_implicit_update
776
777        #Update with given timestep
778        timestep = 0.1
779        quantity.update(0.1)
780
781        denom = 1.0-timestep*semi_implicit_update
782        x = (centroid_values + timestep*explicit_update)/denom
783 
784        assert allclose( quantity.centroid_values, x)
785
786    #Test smoothing
787    def test_smoothing(self):
788
789       
790        from shallow_water import Domain, Transmissive_boundary
791        from Numeric import zeros, Float
792        from anuga.utilities.numerical_tools import mean
793
794
795        #Create shallow water domain
796        domain = Domain(points10)
797        domain.default_order=2
798        domain.reduction = mean
799
800
801        #Set some field values
802        domain.set_quantity('elevation', lambda x: x)
803        domain.set_quantity('friction', 0.03)
804
805
806        ######################
807        # Boundary conditions
808        B = Transmissive_boundary(domain)
809        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
810
811
812        ######################
813        #Initial condition - with jumps
814
815        bed = domain.quantities['elevation'].vertex_values
816        stage = zeros(bed.shape, Float)
817
818        h = 0.03
819        for i in range(stage.shape[0]):
820            if i % 2 == 0:
821                stage[i,:] = bed[i,:] + h
822            else:
823                stage[i,:] = bed[i,:]
824
825        domain.set_quantity('stage', stage)
826
827        stage = domain.quantities['stage']
828
829        #Get smoothed stage
830        A, V = stage.get_vertex_values(xy=False, smooth=True)
831        Q = stage.vertex_values
832
833
834        assert A.shape[0] == 9
835        assert V.shape[0] == 8
836        assert V.shape[1] == 3
837
838        #First four points
839        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
840        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
841        assert allclose(A[2], Q[3,0])
842        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
843
844        #Center point
845        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
846                               Q[5,0] + Q[6,2] + Q[7,1])/6)
847
848
849        #Check V
850        assert allclose(V[0,:], [3,4,0])
851        assert allclose(V[1,:], [1,0,4])
852        assert allclose(V[2,:], [4,5,1])
853        assert allclose(V[3,:], [2,1,5])
854        assert allclose(V[4,:], [6,7,3])
855        assert allclose(V[5,:], [4,3,7])
856        assert allclose(V[6,:], [7,8,4])
857        assert allclose(V[7,:], [5,4,8])
858
859        #Get smoothed stage with XY
860        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
861
862        assert allclose(A, A1)
863        assert allclose(V, V1)
864
865        #Check XY
866        assert allclose(X[4], 0.5)
867        assert allclose(Y[4], 0.5)
868
869        assert allclose(X[7], 1.0)
870        assert allclose(Y[7], 0.5)
871
872
873
874
875    def test_vertex_values_no_smoothing(self):
876
877        from mesh_factory import rectangular
878        from shallow_water import Domain, Transmissive_boundary
879        from Numeric import zeros, Float
880        from anuga.utilities.numerical_tools import mean
881
882
883        #Create basic mesh
884        points, vertices, boundary = rectangular(2, 2)
885
886        #Create shallow water domain
887        domain = Domain(points, vertices, boundary)
888        domain.default_order=2
889        domain.reduction = mean
890
891
892        #Set some field values
893        domain.set_quantity('elevation', lambda x,y: x)
894        domain.set_quantity('friction', 0.03)
895
896
897        ######################
898        #Initial condition - with jumps
899
900        bed = domain.quantities['elevation'].vertex_values
901        stage = zeros(bed.shape, Float)
902
903        h = 0.03
904        for i in range(stage.shape[0]):
905            if i % 2 == 0:
906                stage[i,:] = bed[i,:] + h
907            else:
908                stage[i,:] = bed[i,:]
909
910        domain.set_quantity('stage', stage)
911
912        #Get stage
913        stage = domain.quantities['stage']
914        A, V = stage.get_vertex_values(xy=False, smooth=False)
915        Q = stage.vertex_values.flat
916
917        for k in range(8):
918            assert allclose(A[k], Q[k])
919
920
921        for k in range(8):
922            assert V[k, 0] == 3*k
923            assert V[k, 1] == 3*k+1
924            assert V[k, 2] == 3*k+2
925
926
927
928        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
929
930
931        assert allclose(A, A1)
932        assert allclose(V, V1)
933
934        #Check XY
935        assert allclose(X[1], 0.5)
936        assert allclose(Y[1], 0.5)
937        assert allclose(X[4], 0.0)
938        assert allclose(Y[4], 0.0)
939        assert allclose(X[12], 1.0)
940        assert allclose(Y[12], 0.0)
941
942
943
944 
945#-------------------------------------------------------------
946if __name__ == "__main__":
947    suite = unittest.makeSuite(Test_Quantity, 'test')
948    runner = unittest.TextTestRunner()
949    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.