source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py @ 3956

Last change on this file since 3956 was 3956, checked in by ole, 17 years ago

More cleanup and refactoring. Also adhered to style guide for comparisons to
singletons such as None.

File size: 55.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5
6
7from quantity import *
8from anuga.config import epsilon
9from Numeric import allclose, array, ones, Float
10
11from anuga.fit_interpolate.fit import fit_to_mesh
12#from anuga.pyvolution.least_squares import fit_to_mesh         
13from domain import Domain
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.coordinate_transforms.geo_reference import Geo_reference
16
17#Aux for fit_interpolate.fit example
18def linear_function(point):
19    point = array(point)
20    return point[:,0]+point[:,1]
21
22
23class Test_Quantity(unittest.TestCase):
24    def setUp(self):
25        from domain import Domain
26
27        a = [0.0, 0.0]
28        b = [0.0, 2.0]
29        c = [2.0, 0.0]
30        d = [0.0, 4.0]
31        e = [2.0, 2.0]
32        f = [4.0, 0.0]
33
34        points = [a, b, c, d, e, f]
35
36        #bac, bce, ecf, dbe
37        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
38
39        self.mesh1 = Domain(points[:3], [elements[0]])
40        self.mesh1.check_integrity()
41
42        self.mesh4 = Domain(points, elements)
43        self.mesh4.check_integrity()
44
45    def tearDown(self):
46        pass
47        #print "  Tearing down"
48
49
50    def test_creation(self):
51
52        quantity = Quantity(self.mesh1, [[1,2,3]])
53        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
54
55        try:
56            quantity = Quantity()
57        except:
58            pass
59        else:
60            raise 'Should have raised empty quantity exception'
61
62
63        try:
64            quantity = Quantity([1,2,3])
65        except AssertionError:
66            pass
67        except:
68            raise 'Should have raised "mising mesh object" error'
69
70
71    def test_creation_zeros(self):
72
73        quantity = Quantity(self.mesh1)
74        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
75
76
77        quantity = Quantity(self.mesh4)
78        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
79                                                 [0.,0.,0.], [0.,0.,0.]])
80
81
82    def test_interpolation(self):
83        quantity = Quantity(self.mesh1, [[1,2,3]])
84        assert allclose(quantity.centroid_values, [2.0]) #Centroid
85
86        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
87
88
89    def test_interpolation2(self):
90        quantity = Conserved_quantity(self.mesh4,
91                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
92        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
93
94
95        quantity.extrapolate_second_order()
96
97        #print quantity.vertex_values
98        #assert allclose(quantity.vertex_values, [[2., 2., 2.],
99        #                                         [3.+2./3, 6.+2./3, 4.+2./3],
100        #                                         [7.5, 0.5, 1.],
101        #                                         [-5, -2.5, 7.5]])
102
103        assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
104        #FIXME: Work out the others
105
106
107        #print quantity.edge_values
108        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
109                                               [5., 5., 5.],
110                                               [4.5, 4.5, 0.],
111                                               [3.0, -1.5, -1.5]])
112
113    def test_get_maximum_1(self):
114        quantity = Conserved_quantity(self.mesh4,
115                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
116        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
117
118        v = quantity.get_maximum_value()
119        assert v == 5
120
121        i = quantity.get_maximum_index()
122        assert i == 1
123       
124        x,y = quantity.get_maximum_location()
125        xref, yref = 4.0/3, 4.0/3
126        assert x == xref
127        assert y == yref
128
129        v = quantity.get_values(interpolation_points = [[x,y]])
130        assert allclose(v, 5)
131
132    def test_get_maximum_2(self):
133
134        a = [0.0, 0.0]
135        b = [0.0, 2.0]
136        c = [2.0,0.0]
137        d = [0.0, 4.0]
138        e = [2.0, 2.0]
139        f = [4.0,0.0]
140
141        points = [a, b, c, d, e, f]
142        #bac, bce, ecf, dbe
143        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
144
145        domain = Domain(points, vertices)
146
147        quantity = Quantity(domain)
148        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
149       
150        v = quantity.get_maximum_value()
151        assert v == 6
152
153        i = quantity.get_maximum_index()
154        assert i == 3
155       
156        x,y = quantity.get_maximum_location()
157        xref, yref = 2.0/3, 8.0/3
158        assert x == xref
159        assert y == yref
160
161        v = quantity.get_values(interpolation_points = [[x,y]])
162        assert allclose(v, 6)
163
164
165        #Multiple locations for maximum -
166        #Test that the algorithm picks the first occurrence       
167        v = quantity.get_maximum_value(indices=[0,1,2])
168        assert allclose(v, 4)
169
170        i = quantity.get_maximum_index(indices=[0,1,2])
171        assert i == 1
172       
173        x,y = quantity.get_maximum_location(indices=[0,1,2])
174        xref, yref = 4.0/3, 4.0/3
175        assert x == xref
176        assert y == yref
177
178        v = quantity.get_values(interpolation_points = [[x,y]])
179        assert allclose(v, 4)       
180
181        # More test of indices......
182        v = quantity.get_maximum_value(indices=[2,3])
183        assert allclose(v, 6)
184
185        i = quantity.get_maximum_index(indices=[2,3])
186        assert i == 3
187       
188        x,y = quantity.get_maximum_location(indices=[2,3])
189        xref, yref = 2.0/3, 8.0/3
190        assert x == xref
191        assert y == yref
192
193        v = quantity.get_values(interpolation_points = [[x,y]])
194        assert allclose(v, 6)       
195
196       
197
198    def test_boundary_allocation(self):
199        quantity = Conserved_quantity(self.mesh4,
200                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
201
202        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
203
204
205    def test_set_values(self):
206        quantity = Quantity(self.mesh4)
207
208
209        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
210                            location = 'vertices')
211        assert allclose(quantity.vertex_values,
212                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
213        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
214        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
215                                               [5., 5., 5.],
216                                               [4.5, 4.5, 0.],
217                                               [3.0, -1.5, -1.5]])
218
219
220        #Test default
221        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
222        assert allclose(quantity.vertex_values,
223                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
224        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
225        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
226                                               [5., 5., 5.],
227                                               [4.5, 4.5, 0.],
228                                               [3.0, -1.5, -1.5]])
229
230        #Test centroids
231        quantity.set_values([1,2,3,4], location = 'centroids')
232        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
233
234        #Test edges
235        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
236                            location = 'edges')
237        assert allclose(quantity.edge_values,
238                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
239
240        #Test exceptions
241        try:
242            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
243                                location = 'bas kamel tuba')
244        except:
245            pass
246
247
248        try:
249            quantity.set_values([[1,2,3], [0,0,9]])
250        except AssertionError:
251            pass
252        except:
253            raise 'should have raised Assertionerror'
254
255
256
257    def test_set_values_const(self):
258        quantity = Quantity(self.mesh4)
259
260        quantity.set_values(1.0, location = 'vertices')
261        assert allclose(quantity.vertex_values,
262                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
263
264        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
265        assert allclose(quantity.edge_values, [[1, 1, 1],
266                                               [1, 1, 1],
267                                               [1, 1, 1],
268                                               [1, 1, 1]])
269
270
271        quantity.set_values(2.0, location = 'centroids')
272        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
273
274        quantity.set_values(3.0, location = 'edges')
275        assert allclose(quantity.edge_values, [[3, 3, 3],
276                                               [3, 3, 3],
277                                               [3, 3, 3],
278                                               [3, 3, 3]])
279
280
281    def test_set_values_func(self):
282        quantity = Quantity(self.mesh4)
283
284        def f(x, y):
285            return x+y
286
287        quantity.set_values(f, location = 'vertices')
288        #print "quantity.vertex_values",quantity.vertex_values
289        assert allclose(quantity.vertex_values,
290                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
291        assert allclose(quantity.centroid_values,
292                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
293        assert allclose(quantity.edge_values,
294                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
295
296
297        quantity.set_values(f, location = 'centroids')
298        assert allclose(quantity.centroid_values,
299                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
300
301
302    def test_integral(self):
303        quantity = Quantity(self.mesh4)
304
305        #Try constants first
306        const = 5
307        quantity.set_values(const, location = 'vertices')
308        #print 'Q', quantity.get_integral()
309
310        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
311
312        #Try with a linear function
313        def f(x, y):
314            return x+y
315
316        quantity.set_values(f, location = 'vertices')
317
318
319        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
320
321        assert allclose (quantity.get_integral(), ref_integral)
322
323
324
325    def test_set_vertex_values(self):
326        quantity = Quantity(self.mesh4)
327        quantity.set_vertex_values([0,1,2,3,4,5])
328
329        assert allclose(quantity.vertex_values,
330                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
331        assert allclose(quantity.centroid_values,
332                        [1., 7./3, 11./3, 8./3]) #Centroid
333        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
334                                               [3., 2.5, 1.5],
335                                               [3.5, 4.5, 3.],
336                                               [2.5, 3.5, 2]])
337
338
339    def test_set_vertex_values_subset(self):
340        quantity = Quantity(self.mesh4)
341        quantity.set_vertex_values([0,1,2,3,4,5])
342        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
343
344        assert allclose(quantity.vertex_values,
345                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
346
347
348    def test_set_vertex_values_using_general_interface(self):
349        quantity = Quantity(self.mesh4)
350
351
352        quantity.set_values([0,1,2,3,4,5])
353
354
355        assert allclose(quantity.vertex_values,
356                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
357
358        #Centroid
359        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
360
361        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
362                                               [3., 2.5, 1.5],
363                                               [3.5, 4.5, 3.],
364                                               [2.5, 3.5, 2]])
365
366
367
368
369
370    def test_set_values_using_fit(self):
371
372
373        quantity = Quantity(self.mesh4)
374
375        #Get (enough) datapoints
376        data_points = [[ 0.66666667, 0.66666667],
377                       [ 1.33333333, 1.33333333],
378                       [ 2.66666667, 0.66666667],
379                       [ 0.66666667, 2.66666667],
380                       [ 0.0, 1.0],
381                       [ 0.0, 3.0],
382                       [ 1.0, 0.0],
383                       [ 1.0, 1.0],
384                       [ 1.0, 2.0],
385                       [ 1.0, 3.0],
386                       [ 2.0, 1.0],
387                       [ 3.0, 0.0],
388                       [ 3.0, 1.0]]
389
390        z = linear_function(data_points)
391
392        #Use built-in fit_interpolate.fit
393        quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
394        #quantity.set_values(points = data_points, values = z, alpha = 0)
395
396
397        answer = linear_function(quantity.domain.get_vertex_coordinates())
398        #print quantity.vertex_values, answer
399        assert allclose(quantity.vertex_values.flat, answer)
400
401
402        #Now try by setting the same values directly
403        vertex_attributes = fit_to_mesh(quantity.domain.get_nodes(),
404                                        quantity.domain.triangles, #FIXME
405                                        data_points,
406                                        z,
407                                        alpha = 0,
408                                        verbose=False)
409
410        #print vertex_attributes
411        quantity.set_values(vertex_attributes)
412        assert allclose(quantity.vertex_values.flat, answer)
413
414
415
416
417
418    def test_test_set_values_using_fit_w_geo(self):
419
420
421        #Mesh
422        vertex_coordinates = [[0.76, 0.76],
423                              [0.76, 5.76],
424                              [5.76, 0.76]]
425        triangles = [[0,2,1]]
426
427        mesh_georef = Geo_reference(56,-0.76,-0.76)
428        mesh1 = Domain(vertex_coordinates, triangles,
429                       geo_reference = mesh_georef)
430        mesh1.check_integrity()
431
432        #Quantity
433        quantity = Quantity(mesh1)
434
435        #Data
436        data_points = [[ 201.0, 401.0],
437                       [ 201.0, 403.0],
438                       [ 203.0, 401.0]]
439
440        z = [2, 4, 4]
441
442        data_georef = Geo_reference(56,-200,-400)
443
444
445        #Reference
446        ref = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
447                          data_origin = data_georef.get_origin(),
448                          mesh_origin = mesh_georef.get_origin(),
449                          alpha = 0)
450
451        assert allclose( ref, [0,5,5] )
452
453
454        #Test set_values
455
456        quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
457
458        #quantity.set_values(points = data_points,
459        #                    values = z,
460        #                    data_georef = data_georef,
461        #                    alpha = 0)
462
463
464        #quantity.set_values(points = data_points,
465        #                    values = z,
466        #                    data_georef = data_georef,
467        #                    alpha = 0)
468        assert allclose(quantity.vertex_values.flat, ref)
469
470
471
472        #Test set_values using geospatial data object
473        quantity.vertex_values[:] = 0.0
474
475        geo = Geospatial_data(data_points, z, data_georef)
476
477
478        quantity.set_values(geospatial_data = geo, alpha = 0)
479        assert allclose(quantity.vertex_values.flat, ref)
480
481
482
483    def test_set_values_from_file1(self):
484        quantity = Quantity(self.mesh4)
485
486        #Get (enough) datapoints
487        data_points = [[ 0.66666667, 0.66666667],
488                       [ 1.33333333, 1.33333333],
489                       [ 2.66666667, 0.66666667],
490                       [ 0.66666667, 2.66666667],
491                       [ 0.0, 1.0],
492                       [ 0.0, 3.0],
493                       [ 1.0, 0.0],
494                       [ 1.0, 1.0],
495                       [ 1.0, 2.0],
496                       [ 1.0, 3.0],
497                       [ 2.0, 1.0],
498                       [ 3.0, 0.0],
499                       [ 3.0, 1.0]]
500
501        z = linear_function(data_points)
502
503
504        #Create pts file
505        from load_mesh.loadASCII import export_points_file
506        ptsfile = 'testptsfile.pts'
507        att = 'spam_and_eggs'
508        points_dict = {'pointlist': data_points,
509                       'attributelist': {att: z}}
510
511        export_points_file(ptsfile, points_dict)
512
513        #Check that values can be set from file
514        quantity.set_values(filename = ptsfile,
515                            attribute_name = att, alpha = 0)
516        answer = linear_function(quantity.domain.get_vertex_coordinates())
517
518        #print quantity.vertex_values.flat
519        #print answer
520
521
522        assert allclose(quantity.vertex_values.flat, answer)
523
524
525        #Check that values can be set from file using default attribute
526        quantity.set_values(filename = ptsfile, alpha = 0)
527        assert allclose(quantity.vertex_values.flat, answer)
528
529        #Cleanup
530        import os
531        os.remove(ptsfile)
532
533
534    def test_set_values_from_file_with_georef1(self):
535
536        #Mesh in zone 56 (absolute coords)
537
538        x0 = 314036.58727982
539        y0 = 6224951.2960092
540
541        a = [x0+0.0, y0+0.0]
542        b = [x0+0.0, y0+2.0]
543        c = [x0+2.0, y0+0.0]
544        d = [x0+0.0, y0+4.0]
545        e = [x0+2.0, y0+2.0]
546        f = [x0+4.0, y0+0.0]
547
548        points = [a, b, c, d, e, f]
549
550        #bac, bce, ecf, dbe
551        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
552
553        mesh4 = Domain(points, elements,
554                       geo_reference = Geo_reference(56, 0, 0))
555        mesh4.check_integrity()
556        quantity = Quantity(mesh4)
557
558        #Get (enough) datapoints (relative to georef)
559        data_points = [[ 0.66666667, 0.66666667],
560                       [ 1.33333333, 1.33333333],
561                       [ 2.66666667, 0.66666667],
562                       [ 0.66666667, 2.66666667],
563                       [ 0.0, 1.0],
564                       [ 0.0, 3.0],
565                       [ 1.0, 0.0],
566                       [ 1.0, 1.0],
567                       [ 1.0, 2.0],
568                       [ 1.0, 3.0],
569                       [ 2.0, 1.0],
570                       [ 3.0, 0.0],
571                       [ 3.0, 1.0]]
572
573        z = linear_function(data_points)
574
575
576        #Create pts file
577        from load_mesh.loadASCII import export_points_file
578
579        ptsfile = 'testptsfile.pts'
580        att = 'spam_and_eggs'
581
582        points_dict = {'pointlist': data_points,
583                       'attributelist': {att: z},
584                       'geo_reference': Geo_reference(zone = 56,
585                                                      xllcorner = x0,
586                                                      yllcorner = y0)}
587
588        export_points_file(ptsfile, points_dict)
589
590
591        #Check that values can be set from file
592        quantity.set_values(filename = ptsfile,
593                            attribute_name = att, alpha = 0)
594        answer = linear_function(quantity.domain.get_vertex_coordinates() - [x0, y0])
595
596        assert allclose(quantity.vertex_values.flat, answer)
597
598
599        #Check that values can be set from file using default attribute
600        quantity.set_values(filename = ptsfile, alpha = 0)
601        assert allclose(quantity.vertex_values.flat, answer)
602
603        #Cleanup
604        import os
605        os.remove(ptsfile)
606
607
608    def test_set_values_from_file_with_georef2(self):
609
610        #Mesh in zone 56 (relative coords)
611
612        x0 = 314036.58727982
613        y0 = 6224951.2960092
614
615        a = [0.0, 0.0]
616        b = [0.0, 2.0]
617        c = [2.0, 0.0]
618        d = [0.0, 4.0]
619        e = [2.0, 2.0]
620        f = [4.0, 0.0]
621
622        points = [a, b, c, d, e, f]
623
624        #bac, bce, ecf, dbe
625        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
626
627        mesh4 = Domain(points, elements,
628                       geo_reference = Geo_reference(56, x0, y0))
629        mesh4.check_integrity()
630        quantity = Quantity(mesh4)
631
632        #Get (enough) datapoints (relative to georef)
633        data_points = [[ x0+0.66666667, y0+0.66666667],
634                       [ x0+1.33333333, y0+1.33333333],
635                       [ x0+2.66666667, y0+0.66666667],
636                       [ x0+0.66666667, y0+2.66666667],
637                       [ x0+0.0, y0+1.0],
638                       [ x0+0.0, y0+3.0],
639                       [ x0+1.0, y0+0.0],
640                       [ x0+1.0, y0+1.0],
641                       [ x0+1.0, y0+2.0],
642                       [ x0+1.0, y0+3.0],
643                       [ x0+2.0, y0+1.0],
644                       [ x0+3.0, y0+0.0],
645                       [ x0+3.0, y0+1.0]]
646
647        z = linear_function(data_points)
648
649
650        #Create pts file
651        from load_mesh.loadASCII import export_points_file
652
653        ptsfile = 'testptsfile.pts'
654        att = 'spam_and_eggs'
655
656        points_dict = {'pointlist': data_points,
657                       'attributelist': {att: z},
658                       'geo_reference': Geo_reference(zone = 56,
659                                                      xllcorner = 0,
660                                                      yllcorner = 0)}
661
662        export_points_file(ptsfile, points_dict)
663
664
665        #Check that values can be set from file
666        quantity.set_values(filename = ptsfile,
667                            attribute_name = att, alpha = 0)
668        answer = linear_function(quantity.domain.get_vertex_coordinates() + [x0, y0])
669
670
671        assert allclose(quantity.vertex_values.flat, answer)
672
673
674        #Check that values can be set from file using default attribute
675        quantity.set_values(filename = ptsfile, alpha = 0)
676        assert allclose(quantity.vertex_values.flat, answer)
677
678        #Cleanup
679        import os
680        os.remove(ptsfile)
681
682
683
684
685    def test_set_values_from_quantity(self):
686
687        quantity1 = Quantity(self.mesh4)
688        quantity1.set_vertex_values([0,1,2,3,4,5])
689
690        assert allclose(quantity1.vertex_values,
691                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
692
693
694        quantity2 = Quantity(self.mesh4)
695        quantity2.set_values(quantity = quantity1)
696        assert allclose(quantity2.vertex_values,
697                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
698
699        quantity2.set_values(quantity = 2*quantity1)
700        assert allclose(quantity2.vertex_values,
701                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
702
703        quantity2.set_values(quantity = 2*quantity1 + 3)
704        assert allclose(quantity2.vertex_values,
705                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
706
707
708        #Check detection of quantity as first orgument
709        quantity2.set_values(2*quantity1 + 3)
710        assert allclose(quantity2.vertex_values,
711                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
712
713
714
715
716
717    def test_overloading(self):
718
719        quantity1 = Quantity(self.mesh4)
720        quantity1.set_vertex_values([0,1,2,3,4,5])
721
722        assert allclose(quantity1.vertex_values,
723                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
724
725
726        quantity2 = Quantity(self.mesh4)
727        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
728                             location = 'vertices')
729
730
731
732        quantity3 = Quantity(self.mesh4)
733        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
734                             location = 'vertices')
735
736
737        #Negation
738        Q = -quantity1
739        assert allclose(Q.vertex_values, -quantity1.vertex_values)
740        assert allclose(Q.centroid_values, -quantity1.centroid_values)
741        assert allclose(Q.edge_values, -quantity1.edge_values)
742
743        #Addition
744        Q = quantity1 + 7
745        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
746        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
747        assert allclose(Q.edge_values, quantity1.edge_values + 7)
748
749        Q = 7 + quantity1
750        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
751        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
752        assert allclose(Q.edge_values, quantity1.edge_values + 7)
753
754        Q = quantity1 + quantity2
755        assert allclose(Q.vertex_values,
756                        quantity1.vertex_values + quantity2.vertex_values)
757        assert allclose(Q.centroid_values,
758                        quantity1.centroid_values + quantity2.centroid_values)
759        assert allclose(Q.edge_values,
760                        quantity1.edge_values + quantity2.edge_values)
761
762
763        Q = quantity1 + quantity2 - 3
764        assert allclose(Q.vertex_values,
765                        quantity1.vertex_values + quantity2.vertex_values - 3)
766
767        Q = quantity1 - quantity2
768        assert allclose(Q.vertex_values,
769                        quantity1.vertex_values - quantity2.vertex_values)
770
771        #Scaling
772        Q = quantity1*3
773        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
774        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
775        assert allclose(Q.edge_values, quantity1.edge_values*3)
776        Q = 3*quantity1
777        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
778
779        #Multiplication
780        Q = quantity1 * quantity2
781        #print Q.vertex_values
782        #print Q.centroid_values
783        #print quantity1.centroid_values
784        #print quantity2.centroid_values
785
786        assert allclose(Q.vertex_values,
787                        quantity1.vertex_values * quantity2.vertex_values)
788
789        #Linear combinations
790        Q = 4*quantity1 + 2
791        assert allclose(Q.vertex_values,
792                        4*quantity1.vertex_values + 2)
793
794        Q = quantity1*quantity2 + 2
795        assert allclose(Q.vertex_values,
796                        quantity1.vertex_values * quantity2.vertex_values + 2)
797
798        Q = quantity1*quantity2 + quantity3
799        assert allclose(Q.vertex_values,
800                        quantity1.vertex_values * quantity2.vertex_values +
801                        quantity3.vertex_values)
802        Q = quantity1*quantity2 + 3*quantity3
803        assert allclose(Q.vertex_values,
804                        quantity1.vertex_values * quantity2.vertex_values +
805                        3*quantity3.vertex_values)
806        Q = quantity1*quantity2 + 3*quantity3 + 5.0
807        assert allclose(Q.vertex_values,
808                        quantity1.vertex_values * quantity2.vertex_values +
809                        3*quantity3.vertex_values + 5)
810
811        Q = quantity1*quantity2 - quantity3
812        assert allclose(Q.vertex_values,
813                        quantity1.vertex_values * quantity2.vertex_values -
814                        quantity3.vertex_values)
815        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
816        assert allclose(Q.vertex_values,
817                        1.5*quantity1.vertex_values * quantity2.vertex_values -
818                        3*quantity3.vertex_values + 5)
819
820        #Try combining quantities and arrays and scalars
821        Q = 1.5*quantity1*quantity2.vertex_values -\
822            3*quantity3.vertex_values + 5.0
823        assert allclose(Q.vertex_values,
824                        1.5*quantity1.vertex_values * quantity2.vertex_values -
825                        3*quantity3.vertex_values + 5)
826
827
828        #Powers
829        Q = quantity1**2
830        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
831
832        Q = quantity1**2 +quantity2**2
833        assert allclose(Q.vertex_values,
834                        quantity1.vertex_values**2 + \
835                        quantity2.vertex_values**2)
836
837        Q = (quantity1**2 +quantity2**2)**0.5
838        assert allclose(Q.vertex_values,
839                        (quantity1.vertex_values**2 + \
840                        quantity2.vertex_values**2)**0.5)
841
842
843
844
845
846
847
848    def test_gradient(self):
849        quantity = Conserved_quantity(self.mesh4)
850
851        #Set up for a gradient of (2,0) at mid triangle
852        quantity.set_values([2.0, 4.0, 6.0, 2.0],
853                            location = 'centroids')
854
855
856        #Gradients
857        a, b = quantity.compute_gradients()
858
859        #print self.mesh4.centroid_coordinates
860        #print a, b
861
862        #The central triangle (1)
863        #(using standard gradient based on neigbours controid values)
864        assert allclose(a[1], 2.0)
865        assert allclose(b[1], 0.0)
866
867
868        #Left triangle (0) using two point gradient
869        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
870        #2  = 4  + a*(-2/3)  + b*(-2/3)
871        assert allclose(a[0] + b[0], 3)
872        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
873        assert allclose(a[0] - b[0], 0)
874
875
876        #Right triangle (2) using two point gradient
877        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
878        #6  = 4  + a*(4/3)  + b*(-2/3)
879        assert allclose(2*a[2] - b[2], 3)
880        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
881        assert allclose(a[2] + 2*b[2], 0)
882
883
884        #Top triangle (3) using two point gradient
885        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
886        #2  = 4  + a*(-2/3)  + b*(4/3)
887        assert allclose(a[3] - 2*b[3], 3)
888        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
889        assert allclose(2*a[3] + b[3], 0)
890
891
892
893        #print a, b
894        quantity.extrapolate_second_order()
895
896        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
897        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
898        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
899
900
901        #a = 1.2, b=-0.6
902        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
903        assert allclose(quantity.vertex_values[2,2], 8)
904
905
906    def test_second_order_extrapolation2(self):
907        quantity = Conserved_quantity(self.mesh4)
908
909        #Set up for a gradient of (3,1), f(x) = 3x+y
910        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
911                            location = 'centroids')
912
913        #Gradients
914        a, b = quantity.compute_gradients()
915
916        #print a, b
917
918        assert allclose(a[1], 3.0)
919        assert allclose(b[1], 1.0)
920
921        #Work out the others
922
923        quantity.extrapolate_second_order()
924
925        #print quantity.vertex_values
926        assert allclose(quantity.vertex_values[1,0], 2.0)
927        assert allclose(quantity.vertex_values[1,1], 6.0)
928        assert allclose(quantity.vertex_values[1,2], 8.0)
929
930
931
932    def test_first_order_extrapolator(self):
933        quantity = Conserved_quantity(self.mesh4)
934
935        #Test centroids
936        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
937        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
938
939        #Extrapolate
940        quantity.extrapolate_first_order()
941
942        #Check vertices but not edge values
943        assert allclose(quantity.vertex_values,
944                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
945
946
947    def test_second_order_extrapolator(self):
948        quantity = Conserved_quantity(self.mesh4)
949
950        #Set up for a gradient of (3,0) at mid triangle
951        quantity.set_values([2.0, 4.0, 8.0, 2.0],
952                            location = 'centroids')
953
954
955
956        quantity.extrapolate_second_order()
957        quantity.limit()
958
959
960        #Assert that central triangle is limited by neighbours
961        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
962        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
963
964        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
965        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
966
967        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
968        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
969
970
971        #Assert that quantities are conserved
972        from Numeric import sum
973        for k in range(quantity.centroid_values.shape[0]):
974            assert allclose (quantity.centroid_values[k],
975                             sum(quantity.vertex_values[k,:])/3)
976
977
978
979
980
981    def test_limiter(self):
982        quantity = Conserved_quantity(self.mesh4)
983
984        #Create a deliberate overshoot (e.g. from gradient computation)
985        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
986
987
988        #Limit
989        quantity.limit()
990
991        #Assert that central triangle is limited by neighbours
992        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
993        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
994
995        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
996        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
997
998        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
999        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
1000
1001
1002
1003        #Assert that quantities are conserved
1004        from Numeric import sum
1005        for k in range(quantity.centroid_values.shape[0]):
1006            assert allclose (quantity.centroid_values[k],
1007                             sum(quantity.vertex_values[k,:])/3)
1008
1009
1010    def test_limiter2(self):
1011        """Taken from test_shallow_water
1012        """
1013        quantity = Conserved_quantity(self.mesh4)
1014        quantity.domain.beta_w = 0.9
1015       
1016        #Test centroids
1017        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1018        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1019
1020
1021        #Extrapolate
1022        quantity.extrapolate_second_order()
1023
1024        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1025
1026        #Limit
1027        quantity.limit()
1028
1029        # limited value for beta_w = 0.9
1030       
1031        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1032        # limited values for beta_w = 0.5
1033        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1034
1035
1036        #Assert that quantities are conserved
1037        from Numeric import sum
1038        for k in range(quantity.centroid_values.shape[0]):
1039            assert allclose (quantity.centroid_values[k],
1040                             sum(quantity.vertex_values[k,:])/3)
1041
1042
1043
1044
1045
1046    def test_distribute_first_order(self):
1047        quantity = Conserved_quantity(self.mesh4)
1048
1049        #Test centroids
1050        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1051        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1052
1053
1054        #Extrapolate
1055        quantity.extrapolate_first_order()
1056
1057        #Interpolate
1058        quantity.interpolate_from_vertices_to_edges()
1059
1060        assert allclose(quantity.vertex_values,
1061                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1062        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1063                                               [3,3,3], [4, 4, 4]])
1064
1065
1066    def test_distribute_second_order(self):
1067        quantity = Conserved_quantity(self.mesh4)
1068
1069        #Test centroids
1070        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1071        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1072
1073
1074        #Extrapolate
1075        quantity.extrapolate_second_order()
1076
1077        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1078
1079
1080    def test_update_explicit(self):
1081        quantity = Conserved_quantity(self.mesh4)
1082
1083        #Test centroids
1084        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1085        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1086
1087        #Set explicit_update
1088        quantity.explicit_update = array( [1.,1.,1.,1.] )
1089
1090        #Update with given timestep
1091        quantity.update(0.1)
1092
1093        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1094        assert allclose( quantity.centroid_values, x)
1095
1096    def test_update_semi_implicit(self):
1097        quantity = Conserved_quantity(self.mesh4)
1098
1099        #Test centroids
1100        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1101        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1102
1103        #Set semi implicit update
1104        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1105
1106        #Update with given timestep
1107        timestep = 0.1
1108        quantity.update(timestep)
1109
1110        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1111        denom = ones(4, Float)-timestep*sem
1112
1113        x = array([1, 2, 3, 4])/denom
1114        assert allclose( quantity.centroid_values, x)
1115
1116
1117    def test_both_updates(self):
1118        quantity = Conserved_quantity(self.mesh4)
1119
1120        #Test centroids
1121        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1122        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1123
1124        #Set explicit_update
1125        quantity.explicit_update = array( [4.,3.,2.,1.] )
1126
1127        #Set semi implicit update
1128        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1129
1130        #Update with given timestep
1131        timestep = 0.1
1132        quantity.update(0.1)
1133
1134        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1135        denom = ones(4, Float)-timestep*sem
1136
1137        x = array([1., 2., 3., 4.])
1138        x /= denom
1139        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1140
1141        assert allclose( quantity.centroid_values, x)
1142
1143
1144
1145
1146    #Test smoothing
1147    def test_smoothing(self):
1148
1149        from mesh_factory import rectangular
1150        from shallow_water import Domain, Transmissive_boundary
1151        from Numeric import zeros, Float
1152        from anuga.utilities.numerical_tools import mean
1153
1154        #Create basic mesh
1155        points, vertices, boundary = rectangular(2, 2)
1156
1157        #Create shallow water domain
1158        domain = Domain(points, vertices, boundary)
1159        domain.default_order=2
1160        domain.reduction = mean
1161
1162
1163        #Set some field values
1164        domain.set_quantity('elevation', lambda x,y: x)
1165        domain.set_quantity('friction', 0.03)
1166
1167
1168        ######################
1169        # Boundary conditions
1170        B = Transmissive_boundary(domain)
1171        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1172
1173
1174        ######################
1175        #Initial condition - with jumps
1176
1177        bed = domain.quantities['elevation'].vertex_values
1178        stage = zeros(bed.shape, Float)
1179
1180        h = 0.03
1181        for i in range(stage.shape[0]):
1182            if i % 2 == 0:
1183                stage[i,:] = bed[i,:] + h
1184            else:
1185                stage[i,:] = bed[i,:]
1186
1187        domain.set_quantity('stage', stage)
1188
1189        stage = domain.quantities['stage']
1190
1191        #Get smoothed stage
1192        A, V = stage.get_vertex_values(xy=False, smooth=True)
1193        Q = stage.vertex_values
1194
1195
1196        assert A.shape[0] == 9
1197        assert V.shape[0] == 8
1198        assert V.shape[1] == 3
1199
1200        #First four points
1201        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1202        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1203        assert allclose(A[2], Q[3,0])
1204        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1205
1206        #Center point
1207        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1208                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1209
1210
1211        #Check V
1212        assert allclose(V[0,:], [3,4,0])
1213        assert allclose(V[1,:], [1,0,4])
1214        assert allclose(V[2,:], [4,5,1])
1215        assert allclose(V[3,:], [2,1,5])
1216        assert allclose(V[4,:], [6,7,3])
1217        assert allclose(V[5,:], [4,3,7])
1218        assert allclose(V[6,:], [7,8,4])
1219        assert allclose(V[7,:], [5,4,8])
1220
1221        #Get smoothed stage with XY
1222        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1223
1224        assert allclose(A, A1)
1225        assert allclose(V, V1)
1226
1227        #Check XY
1228        assert allclose(X[4], 0.5)
1229        assert allclose(Y[4], 0.5)
1230
1231        assert allclose(X[7], 1.0)
1232        assert allclose(Y[7], 0.5)
1233
1234
1235
1236
1237    def test_vertex_values_no_smoothing(self):
1238
1239        from mesh_factory import rectangular
1240        from shallow_water import Domain, Transmissive_boundary
1241        from Numeric import zeros, Float
1242        from anuga.utilities.numerical_tools import mean
1243
1244
1245        #Create basic mesh
1246        points, vertices, boundary = rectangular(2, 2)
1247
1248        #Create shallow water domain
1249        domain = Domain(points, vertices, boundary)
1250        domain.default_order=2
1251        domain.reduction = mean
1252
1253
1254        #Set some field values
1255        domain.set_quantity('elevation', lambda x,y: x)
1256        domain.set_quantity('friction', 0.03)
1257
1258
1259        ######################
1260        #Initial condition - with jumps
1261
1262        bed = domain.quantities['elevation'].vertex_values
1263        stage = zeros(bed.shape, Float)
1264
1265        h = 0.03
1266        for i in range(stage.shape[0]):
1267            if i % 2 == 0:
1268                stage[i,:] = bed[i,:] + h
1269            else:
1270                stage[i,:] = bed[i,:]
1271
1272        domain.set_quantity('stage', stage)
1273
1274        #Get stage
1275        stage = domain.quantities['stage']
1276        A, V = stage.get_vertex_values(xy=False, smooth=False)
1277        Q = stage.vertex_values.flat
1278
1279        for k in range(8):
1280            assert allclose(A[k], Q[k])
1281
1282
1283        for k in range(8):
1284            assert V[k, 0] == 3*k
1285            assert V[k, 1] == 3*k+1
1286            assert V[k, 2] == 3*k+2
1287
1288
1289
1290        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1291
1292
1293        assert allclose(A, A1)
1294        assert allclose(V, V1)
1295
1296        #Check XY
1297        assert allclose(X[1], 0.5)
1298        assert allclose(Y[1], 0.5)
1299        assert allclose(X[4], 0.0)
1300        assert allclose(Y[4], 0.0)
1301        assert allclose(X[12], 1.0)
1302        assert allclose(Y[12], 0.0)
1303
1304
1305
1306    def set_array_values_by_index(self):
1307
1308        from mesh_factory import rectangular
1309        from shallow_water import Domain
1310        from Numeric import zeros, Float
1311
1312        #Create basic mesh
1313        points, vertices, boundary = rectangular(1, 1)
1314
1315        #Create shallow water domain
1316        domain = Domain(points, vertices, boundary)
1317        #print "domain.number_of_elements ",domain.number_of_elements
1318        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1319        value = [7]
1320        indices = [1]
1321        quantity.set_array_values_by_index(value,
1322                                           location = 'centroids',
1323                                           indices = indices)
1324        #print "quantity.centroid_values",quantity.centroid_values
1325
1326        assert allclose(quantity.centroid_values, [1,7])
1327
1328        quantity.set_array_values([15,20,25], indices = indices)
1329        assert allclose(quantity.centroid_values, [1,20])
1330
1331        quantity.set_array_values([15,20,25], indices = indices)
1332        assert allclose(quantity.centroid_values, [1,20])
1333
1334    def test_setting_some_vertex_values(self):
1335        """
1336        set values based on triangle lists.
1337        """
1338        from mesh_factory import rectangular
1339        from shallow_water import Domain
1340        from Numeric import zeros, Float
1341
1342        #Create basic mesh
1343        points, vertices, boundary = rectangular(1, 3)
1344        #print "vertices",vertices
1345        #Create shallow water domain
1346        domain = Domain(points, vertices, boundary)
1347        #print "domain.number_of_elements ",domain.number_of_elements
1348        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1349                                    [4,4,4],[5,5,5],[6,6,6]])
1350
1351
1352        # Check that constants work
1353        value = 7
1354        indices = [1]
1355        quantity.set_values(value,
1356                            location = 'centroids',
1357                            indices = indices)
1358        #print "quantity.centroid_values",quantity.centroid_values
1359        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1360       
1361        value = [7]
1362        indices = [1]
1363        quantity.set_values(value,
1364                            location = 'centroids',
1365                            indices = indices)
1366        #print "quantity.centroid_values",quantity.centroid_values
1367        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1368
1369        value = [[15,20,25]]
1370        quantity.set_values(value, indices = indices)
1371        #print "1 quantity.vertex_values",quantity.vertex_values
1372        assert allclose(quantity.vertex_values[1], value[0])
1373
1374
1375        #print "quantity",quantity.vertex_values
1376        values = [10,100,50]
1377        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1378        #print "2 quantity.vertex_values",quantity.vertex_values
1379        assert allclose(quantity.vertex_values[0], [10,10,10])
1380        assert allclose(quantity.vertex_values[5], [50,50,50])
1381        #quantity.interpolate()
1382        #print "quantity.centroid_values",quantity.centroid_values
1383        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1384
1385
1386        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1387                                    [4,4,4],[5,5,5],[6,6,6]])
1388        values = [10,100,50]
1389        #this will be per unique vertex, indexing the vertices
1390        #print "quantity.vertex_values",quantity.vertex_values
1391        quantity.set_values(values, indices = [0,1,5])
1392        #print "quantity.vertex_values",quantity.vertex_values
1393        assert allclose(quantity.vertex_values[0], [1,50,10])
1394        assert allclose(quantity.vertex_values[5], [6,6,6])
1395        assert allclose(quantity.vertex_values[1], [100,10,50])
1396
1397        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1398                                    [4,4,4],[5,5,5],[6,6,6]])
1399        values = [[31,30,29],[400,400,400],[1000,999,998]]
1400        quantity.set_values(values, indices = [3,3,5])
1401        quantity.interpolate()
1402        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1403
1404        values = [[1,1,1],[2,2,2],[3,3,3],
1405                                    [4,4,4],[5,5,5],[6,6,6]]
1406        quantity.set_values(values)
1407
1408        # testing the standard set values by vertex
1409        # indexed by vertex_id in general_mesh.coordinates
1410        values = [0,1,2,3,4,5,6,7]
1411
1412        quantity.set_values(values)
1413        #print "1 quantity.vertex_values",quantity.vertex_values
1414        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1415                                                [ 1.,  0.,  5.],
1416                                                [ 5.,  6.,  1.],
1417                                                [ 2.,  1.,  6.],
1418                                                [ 6.,  7.,  2.],
1419                                                [ 3.,  2.,  7.]])
1420
1421    def test_setting_unique_vertex_values(self):
1422        """
1423        set values based on unique_vertex lists.
1424        """
1425        from mesh_factory import rectangular
1426        from shallow_water import Domain
1427        from Numeric import zeros, Float
1428
1429        #Create basic mesh
1430        points, vertices, boundary = rectangular(1, 3)
1431        #print "vertices",vertices
1432        #Create shallow water domain
1433        domain = Domain(points, vertices, boundary)
1434        #print "domain.number_of_elements ",domain.number_of_elements
1435        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1436                                    [4,4,4],[5,5,5]])
1437        value = 7
1438        indices = [1,5]
1439        quantity.set_values(value,
1440                            location = 'unique vertices',
1441                            indices = indices)
1442        #print "quantity.centroid_values",quantity.centroid_values
1443        assert allclose(quantity.vertex_values[0], [0,7,0])
1444        assert allclose(quantity.vertex_values[1], [7,1,7])
1445        assert allclose(quantity.vertex_values[2], [7,2,7])
1446
1447
1448    def test_get_values(self):
1449        """
1450        get values based on triangle lists.
1451        """
1452        from mesh_factory import rectangular
1453        from shallow_water import Domain
1454        from Numeric import zeros, Float
1455
1456        #Create basic mesh
1457        points, vertices, boundary = rectangular(1, 3)
1458
1459        #print "points",points
1460        #print "vertices",vertices
1461        #print "boundary",boundary
1462
1463        #Create shallow water domain
1464        domain = Domain(points, vertices, boundary)
1465        #print "domain.number_of_elements ",domain.number_of_elements
1466        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1467                                    [4,4,4],[5,5,5]])
1468
1469        #print "quantity.get_values(location = 'unique vertices')", \
1470        #      quantity.get_values(location = 'unique vertices')
1471
1472        #print "quantity.get_values(location = 'unique vertices')", \
1473        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1474        #                          location = 'unique vertices')
1475
1476        answer = [0.5,2,4,5,0,1,3,4.5]
1477        assert allclose(answer,
1478                        quantity.get_values(location = 'unique vertices'))
1479
1480        indices = [0,5,3]
1481        answer = [0.5,1,5]
1482        assert allclose(answer,
1483                        quantity.get_values(indices=indices, \
1484                                            location = 'unique vertices'))
1485        #print "quantity.centroid_values",quantity.centroid_values
1486        #print "quantity.get_values(location = 'centroids') ",\
1487        #      quantity.get_values(location = 'centroids')
1488
1489
1490
1491
1492    def test_get_values_2(self):
1493        """Different mesh (working with domain object) - also check centroids.
1494        """
1495
1496       
1497        a = [0.0, 0.0]
1498        b = [0.0, 2.0]
1499        c = [2.0,0.0]
1500        d = [0.0, 4.0]
1501        e = [2.0, 2.0]
1502        f = [4.0,0.0]
1503
1504        points = [a, b, c, d, e, f]
1505        #bac, bce, ecf, dbe
1506        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1507
1508        domain = Domain(points, vertices)
1509
1510        quantity = Quantity(domain)
1511        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1512       
1513        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
1514        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
1515
1516
1517        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
1518                                                                   [4,2,6],
1519                                                                   [6,2,4],
1520                                                                   [8,4,6]])
1521       
1522        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
1523                                                                                  [8,4,6]])
1524
1525
1526        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
1527                                                                [4,5,3],
1528                                                                [3,5,4],
1529                                                                [5,7,6]])
1530        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
1531                        [[4,5,3],
1532                         [5,7,6]])       
1533
1534        # Check averaging over vertices
1535        #a: 0
1536        #b: (4+4+4)/3
1537        #c: (2+2+2)/3
1538        #d: 8
1539        #e: (6+6+6)/3       
1540        #f: 4
1541        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
1542                                                                                 
1543       
1544
1545
1546
1547
1548    def test_get_interpolated_values(self):
1549
1550        from mesh_factory import rectangular
1551        from shallow_water import Domain
1552        from Numeric import zeros, Float
1553
1554        #Create basic mesh
1555        points, vertices, boundary = rectangular(1, 3)
1556        domain = Domain(points, vertices, boundary)
1557
1558        #Constant values
1559        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1560                                    [4,4,4],[5,5,5]])
1561
1562       
1563
1564        # Get interpolated values at centroids
1565        interpolation_points = domain.get_centroid_coordinates()
1566        answer = quantity.get_values(location='centroids')
1567
1568       
1569        #print quantity.get_values(points=interpolation_points)
1570        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
1571
1572
1573        #Arbitrary values
1574        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
1575                                    [1,4,-9],[2,5,0]])
1576
1577
1578        # Get interpolated values at centroids
1579        interpolation_points = domain.get_centroid_coordinates()
1580        answer = quantity.get_values(location='centroids')
1581        #print answer
1582        #print quantity.get_values(points=interpolation_points)
1583        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
1584                       
1585
1586        #FIXME TODO
1587        #indices = [0,5,3]
1588        #answer = [0.5,1,5]
1589        #assert allclose(answer,
1590        #                quantity.get_values(indices=indices, \
1591        #                                    location = 'unique vertices'))
1592
1593
1594
1595
1596    def test_get_interpolated_values_2(self):
1597        a = [0.0, 0.0]
1598        b = [0.0, 2.0]
1599        c = [2.0,0.0]
1600        d = [0.0, 4.0]
1601        e = [2.0, 2.0]
1602        f = [4.0,0.0]
1603
1604        points = [a, b, c, d, e, f]
1605        #bac, bce, ecf, dbe
1606        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1607
1608        domain = Domain(points, vertices)
1609
1610        quantity = Quantity(domain)
1611        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1612
1613        #First pick one point
1614        x, y = 2.0/3, 8.0/3
1615        v = quantity.get_values(interpolation_points = [[x,y]])
1616        assert allclose(v, 6)       
1617
1618        # Then another to test that algorithm won't blindly
1619        # reuse interpolation matrix
1620        x, y = 4.0/3, 4.0/3
1621        v = quantity.get_values(interpolation_points = [[x,y]])
1622        assert allclose(v, 4)       
1623
1624
1625
1626
1627    def test_getting_some_vertex_values(self):
1628        """
1629        get values based on triangle lists.
1630        """
1631        from mesh_factory import rectangular
1632        from shallow_water import Domain
1633        from Numeric import zeros, Float
1634
1635        #Create basic mesh
1636        points, vertices, boundary = rectangular(1, 3)
1637
1638        #print "points",points
1639        #print "vertices",vertices
1640        #print "boundary",boundary
1641
1642        #Create shallow water domain
1643        domain = Domain(points, vertices, boundary)
1644        #print "domain.number_of_elements ",domain.number_of_elements
1645        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1646                                    [4,4,4],[5,5,5],[6,6,6]])
1647        value = [7]
1648        indices = [1]
1649        quantity.set_values(value,
1650                            location = 'centroids',
1651                            indices = indices)
1652        #print "quantity.centroid_values",quantity.centroid_values
1653        #print "quantity.get_values(location = 'centroids') ",\
1654        #      quantity.get_values(location = 'centroids')
1655        assert allclose(quantity.centroid_values,
1656                        quantity.get_values(location = 'centroids'))
1657
1658
1659        value = [[15,20,25]]
1660        quantity.set_values(value, indices = indices)
1661        #print "1 quantity.vertex_values",quantity.vertex_values
1662        assert allclose(quantity.vertex_values, quantity.get_values())
1663
1664        assert allclose(quantity.edge_values,
1665                        quantity.get_values(location = 'edges'))
1666
1667        # get a subset of elements
1668        subset = quantity.get_values(location='centroids', indices=[0,5])
1669        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1670        assert allclose(subset, answer)
1671
1672
1673        subset = quantity.get_values(location='edges', indices=[0,5])
1674        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1675        #print "subset",subset
1676        #print "answer",answer
1677        assert allclose(subset, answer)
1678
1679        subset = quantity.get_values( indices=[1,5])
1680        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1681        #print "subset",subset
1682        #print "answer",answer
1683        assert allclose(subset, answer)
1684
1685
1686
1687
1688#-------------------------------------------------------------
1689if __name__ == "__main__":
1690    suite = unittest.makeSuite(Test_Quantity, 'test')
1691    #print "restricted test"
1692    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1693    runner = unittest.TextTestRunner()
1694    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.