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

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

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.

File size: 55.1 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(obj = True))
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.coordinates,
404                                        quantity.domain.triangles,
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(obj = True))
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(obj = True) - [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(obj = True) + [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
1015        #Test centroids
1016        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1017        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1018
1019
1020        #Extrapolate
1021        quantity.extrapolate_second_order()
1022
1023        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1024
1025        #Limit
1026        quantity.limit()
1027
1028        # limited value for beta_w = 0.9
1029        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1030        # limited values for beta_w = 0.5
1031        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1032
1033
1034        #Assert that quantities are conserved
1035        from Numeric import sum
1036        for k in range(quantity.centroid_values.shape[0]):
1037            assert allclose (quantity.centroid_values[k],
1038                             sum(quantity.vertex_values[k,:])/3)
1039
1040
1041
1042
1043
1044    def test_distribute_first_order(self):
1045        quantity = Conserved_quantity(self.mesh4)
1046
1047        #Test centroids
1048        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1049        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1050
1051
1052        #Extrapolate
1053        quantity.extrapolate_first_order()
1054
1055        #Interpolate
1056        quantity.interpolate_from_vertices_to_edges()
1057
1058        assert allclose(quantity.vertex_values,
1059                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1060        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1061                                               [3,3,3], [4, 4, 4]])
1062
1063
1064    def test_distribute_second_order(self):
1065        quantity = Conserved_quantity(self.mesh4)
1066
1067        #Test centroids
1068        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1069        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1070
1071
1072        #Extrapolate
1073        quantity.extrapolate_second_order()
1074
1075        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1076
1077
1078    def test_update_explicit(self):
1079        quantity = Conserved_quantity(self.mesh4)
1080
1081        #Test centroids
1082        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1083        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1084
1085        #Set explicit_update
1086        quantity.explicit_update = array( [1.,1.,1.,1.] )
1087
1088        #Update with given timestep
1089        quantity.update(0.1)
1090
1091        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1092        assert allclose( quantity.centroid_values, x)
1093
1094    def test_update_semi_implicit(self):
1095        quantity = Conserved_quantity(self.mesh4)
1096
1097        #Test centroids
1098        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1099        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1100
1101        #Set semi implicit update
1102        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1103
1104        #Update with given timestep
1105        timestep = 0.1
1106        quantity.update(timestep)
1107
1108        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1109        denom = ones(4, Float)-timestep*sem
1110
1111        x = array([1, 2, 3, 4])/denom
1112        assert allclose( quantity.centroid_values, x)
1113
1114
1115    def test_both_updates(self):
1116        quantity = Conserved_quantity(self.mesh4)
1117
1118        #Test centroids
1119        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1120        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1121
1122        #Set explicit_update
1123        quantity.explicit_update = array( [4.,3.,2.,1.] )
1124
1125        #Set semi implicit update
1126        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1127
1128        #Update with given timestep
1129        timestep = 0.1
1130        quantity.update(0.1)
1131
1132        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1133        denom = ones(4, Float)-timestep*sem
1134
1135        x = array([1., 2., 3., 4.])
1136        x /= denom
1137        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1138
1139        assert allclose( quantity.centroid_values, x)
1140
1141
1142
1143
1144    #Test smoothing
1145    def test_smoothing(self):
1146
1147        from mesh_factory import rectangular
1148        from shallow_water import Domain, Transmissive_boundary
1149        from Numeric import zeros, Float
1150        from anuga.utilities.numerical_tools import mean
1151
1152        #Create basic mesh
1153        points, vertices, boundary = rectangular(2, 2)
1154
1155        #Create shallow water domain
1156        domain = Domain(points, vertices, boundary)
1157        domain.default_order=2
1158        domain.reduction = mean
1159
1160
1161        #Set some field values
1162        domain.set_quantity('elevation', lambda x,y: x)
1163        domain.set_quantity('friction', 0.03)
1164
1165
1166        ######################
1167        # Boundary conditions
1168        B = Transmissive_boundary(domain)
1169        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1170
1171
1172        ######################
1173        #Initial condition - with jumps
1174
1175        bed = domain.quantities['elevation'].vertex_values
1176        stage = zeros(bed.shape, Float)
1177
1178        h = 0.03
1179        for i in range(stage.shape[0]):
1180            if i % 2 == 0:
1181                stage[i,:] = bed[i,:] + h
1182            else:
1183                stage[i,:] = bed[i,:]
1184
1185        domain.set_quantity('stage', stage)
1186
1187        stage = domain.quantities['stage']
1188
1189        #Get smoothed stage
1190        A, V = stage.get_vertex_values(xy=False, smooth=True)
1191        Q = stage.vertex_values
1192
1193
1194        assert A.shape[0] == 9
1195        assert V.shape[0] == 8
1196        assert V.shape[1] == 3
1197
1198        #First four points
1199        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1200        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1201        assert allclose(A[2], Q[3,0])
1202        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1203
1204        #Center point
1205        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1206                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1207
1208
1209        #Check V
1210        assert allclose(V[0,:], [3,4,0])
1211        assert allclose(V[1,:], [1,0,4])
1212        assert allclose(V[2,:], [4,5,1])
1213        assert allclose(V[3,:], [2,1,5])
1214        assert allclose(V[4,:], [6,7,3])
1215        assert allclose(V[5,:], [4,3,7])
1216        assert allclose(V[6,:], [7,8,4])
1217        assert allclose(V[7,:], [5,4,8])
1218
1219        #Get smoothed stage with XY
1220        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1221
1222        assert allclose(A, A1)
1223        assert allclose(V, V1)
1224
1225        #Check XY
1226        assert allclose(X[4], 0.5)
1227        assert allclose(Y[4], 0.5)
1228
1229        assert allclose(X[7], 1.0)
1230        assert allclose(Y[7], 0.5)
1231
1232
1233
1234
1235    def test_vertex_values_no_smoothing(self):
1236
1237        from mesh_factory import rectangular
1238        from shallow_water import Domain, Transmissive_boundary
1239        from Numeric import zeros, Float
1240        from anuga.utilities.numerical_tools import mean
1241
1242
1243        #Create basic mesh
1244        points, vertices, boundary = rectangular(2, 2)
1245
1246        #Create shallow water domain
1247        domain = Domain(points, vertices, boundary)
1248        domain.default_order=2
1249        domain.reduction = mean
1250
1251
1252        #Set some field values
1253        domain.set_quantity('elevation', lambda x,y: x)
1254        domain.set_quantity('friction', 0.03)
1255
1256
1257        ######################
1258        #Initial condition - with jumps
1259
1260        bed = domain.quantities['elevation'].vertex_values
1261        stage = zeros(bed.shape, Float)
1262
1263        h = 0.03
1264        for i in range(stage.shape[0]):
1265            if i % 2 == 0:
1266                stage[i,:] = bed[i,:] + h
1267            else:
1268                stage[i,:] = bed[i,:]
1269
1270        domain.set_quantity('stage', stage)
1271
1272        #Get stage
1273        stage = domain.quantities['stage']
1274        A, V = stage.get_vertex_values(xy=False, smooth=False)
1275        Q = stage.vertex_values.flat
1276
1277        for k in range(8):
1278            assert allclose(A[k], Q[k])
1279
1280
1281        for k in range(8):
1282            assert V[k, 0] == 3*k
1283            assert V[k, 1] == 3*k+1
1284            assert V[k, 2] == 3*k+2
1285
1286
1287
1288        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1289
1290
1291        assert allclose(A, A1)
1292        assert allclose(V, V1)
1293
1294        #Check XY
1295        assert allclose(X[1], 0.5)
1296        assert allclose(Y[1], 0.5)
1297        assert allclose(X[4], 0.0)
1298        assert allclose(Y[4], 0.0)
1299        assert allclose(X[12], 1.0)
1300        assert allclose(Y[12], 0.0)
1301
1302
1303
1304    def set_array_values_by_index(self):
1305
1306        from mesh_factory import rectangular
1307        from shallow_water import Domain
1308        from Numeric import zeros, Float
1309
1310        #Create basic mesh
1311        points, vertices, boundary = rectangular(1, 1)
1312
1313        #Create shallow water domain
1314        domain = Domain(points, vertices, boundary)
1315        #print "domain.number_of_elements ",domain.number_of_elements
1316        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1317        value = [7]
1318        indices = [1]
1319        quantity.set_array_values_by_index(value,
1320                                           location = 'centroids',
1321                                           indices = indices)
1322        #print "quantity.centroid_values",quantity.centroid_values
1323
1324        assert allclose(quantity.centroid_values, [1,7])
1325
1326        quantity.set_array_values([15,20,25], indices = indices)
1327        assert allclose(quantity.centroid_values, [1,20])
1328
1329        quantity.set_array_values([15,20,25], indices = indices)
1330        assert allclose(quantity.centroid_values, [1,20])
1331
1332    def test_setting_some_vertex_values(self):
1333        """
1334        set values based on triangle lists.
1335        """
1336        from mesh_factory import rectangular
1337        from shallow_water import Domain
1338        from Numeric import zeros, Float
1339
1340        #Create basic mesh
1341        points, vertices, boundary = rectangular(1, 3)
1342        #print "vertices",vertices
1343        #Create shallow water domain
1344        domain = Domain(points, vertices, boundary)
1345        #print "domain.number_of_elements ",domain.number_of_elements
1346        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1347                                    [4,4,4],[5,5,5],[6,6,6]])
1348        value = [7]
1349        indices = [1]
1350        quantity.set_values(value,
1351                            location = 'centroids',
1352                            indices = indices)
1353        #print "quantity.centroid_values",quantity.centroid_values
1354        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1355
1356        value = [[15,20,25]]
1357        quantity.set_values(value, indices = indices)
1358        #print "1 quantity.vertex_values",quantity.vertex_values
1359        assert allclose(quantity.vertex_values[1], value[0])
1360
1361
1362        #print "quantity",quantity.vertex_values
1363        values = [10,100,50]
1364        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1365        #print "2 quantity.vertex_values",quantity.vertex_values
1366        assert allclose(quantity.vertex_values[0], [10,10,10])
1367        assert allclose(quantity.vertex_values[5], [50,50,50])
1368        #quantity.interpolate()
1369        #print "quantity.centroid_values",quantity.centroid_values
1370        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1371
1372
1373        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1374                                    [4,4,4],[5,5,5],[6,6,6]])
1375        values = [10,100,50]
1376        #this will be per unique vertex, indexing the vertices
1377        #print "quantity.vertex_values",quantity.vertex_values
1378        quantity.set_values(values, indices = [0,1,5])
1379        #print "quantity.vertex_values",quantity.vertex_values
1380        assert allclose(quantity.vertex_values[0], [1,50,10])
1381        assert allclose(quantity.vertex_values[5], [6,6,6])
1382        assert allclose(quantity.vertex_values[1], [100,10,50])
1383
1384        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1385                                    [4,4,4],[5,5,5],[6,6,6]])
1386        values = [[31,30,29],[400,400,400],[1000,999,998]]
1387        quantity.set_values(values, indices = [3,3,5])
1388        quantity.interpolate()
1389        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1390
1391        values = [[1,1,1],[2,2,2],[3,3,3],
1392                                    [4,4,4],[5,5,5],[6,6,6]]
1393        quantity.set_values(values)
1394
1395        # testing the standard set values by vertex
1396        # indexed by vertex_id in general_mesh.coordinates
1397        values = [0,1,2,3,4,5,6,7]
1398
1399        quantity.set_values(values)
1400        #print "1 quantity.vertex_values",quantity.vertex_values
1401        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1402                                                [ 1.,  0.,  5.],
1403                                                [ 5.,  6.,  1.],
1404                                                [ 2.,  1.,  6.],
1405                                                [ 6.,  7.,  2.],
1406                                                [ 3.,  2.,  7.]])
1407
1408    def test_setting_unique_vertex_values(self):
1409        """
1410        set values based on unique_vertex lists.
1411        """
1412        from mesh_factory import rectangular
1413        from shallow_water import Domain
1414        from Numeric import zeros, Float
1415
1416        #Create basic mesh
1417        points, vertices, boundary = rectangular(1, 3)
1418        #print "vertices",vertices
1419        #Create shallow water domain
1420        domain = Domain(points, vertices, boundary)
1421        #print "domain.number_of_elements ",domain.number_of_elements
1422        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1423                                    [4,4,4],[5,5,5]])
1424        value = 7
1425        indices = [1,5]
1426        quantity.set_values(value,
1427                            location = 'unique vertices',
1428                            indices = indices)
1429        #print "quantity.centroid_values",quantity.centroid_values
1430        assert allclose(quantity.vertex_values[0], [0,7,0])
1431        assert allclose(quantity.vertex_values[1], [7,1,7])
1432        assert allclose(quantity.vertex_values[2], [7,2,7])
1433
1434
1435    def test_get_values(self):
1436        """
1437        get values based on triangle lists.
1438        """
1439        from mesh_factory import rectangular
1440        from shallow_water import Domain
1441        from Numeric import zeros, Float
1442
1443        #Create basic mesh
1444        points, vertices, boundary = rectangular(1, 3)
1445
1446        #print "points",points
1447        #print "vertices",vertices
1448        #print "boundary",boundary
1449
1450        #Create shallow water domain
1451        domain = Domain(points, vertices, boundary)
1452        #print "domain.number_of_elements ",domain.number_of_elements
1453        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1454                                    [4,4,4],[5,5,5]])
1455
1456        #print "quantity.get_values(location = 'unique vertices')", \
1457        #      quantity.get_values(location = 'unique vertices')
1458
1459        #print "quantity.get_values(location = 'unique vertices')", \
1460        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1461        #                          location = 'unique vertices')
1462
1463        answer = [0.5,2,4,5,0,1,3,4.5]
1464        assert allclose(answer,
1465                        quantity.get_values(location = 'unique vertices'))
1466
1467        indices = [0,5,3]
1468        answer = [0.5,1,5]
1469        assert allclose(answer,
1470                        quantity.get_values(indices=indices, \
1471                                            location = 'unique vertices'))
1472        #print "quantity.centroid_values",quantity.centroid_values
1473        #print "quantity.get_values(location = 'centroids') ",\
1474        #      quantity.get_values(location = 'centroids')
1475
1476
1477
1478
1479    def test_get_values_2(self):
1480        """Different mesh (working with domain object) - also check centroids.
1481        """
1482
1483       
1484        a = [0.0, 0.0]
1485        b = [0.0, 2.0]
1486        c = [2.0,0.0]
1487        d = [0.0, 4.0]
1488        e = [2.0, 2.0]
1489        f = [4.0,0.0]
1490
1491        points = [a, b, c, d, e, f]
1492        #bac, bce, ecf, dbe
1493        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1494
1495        domain = Domain(points, vertices)
1496
1497        quantity = Quantity(domain)
1498        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1499       
1500        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
1501        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
1502
1503
1504        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
1505                                                                   [4,2,6],
1506                                                                   [6,2,4],
1507                                                                   [8,4,6]])
1508       
1509        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
1510                                                                                  [8,4,6]])
1511
1512
1513        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
1514                                                                [4,5,3],
1515                                                                [3,5,4],
1516                                                                [5,7,6]])
1517        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
1518                        [[4,5,3],
1519                         [5,7,6]])       
1520
1521        # Check averaging over vertices
1522        #a: 0
1523        #b: (4+4+4)/3
1524        #c: (2+2+2)/3
1525        #d: 8
1526        #e: (6+6+6)/3       
1527        #f: 4
1528        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
1529                                                                                 
1530       
1531
1532
1533
1534
1535    def test_get_interpolated_values(self):
1536
1537        from mesh_factory import rectangular
1538        from shallow_water import Domain
1539        from Numeric import zeros, Float
1540
1541        #Create basic mesh
1542        points, vertices, boundary = rectangular(1, 3)
1543        domain = Domain(points, vertices, boundary)
1544
1545        #Constant values
1546        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1547                                    [4,4,4],[5,5,5]])
1548
1549       
1550
1551        # Get interpolated values at centroids
1552        interpolation_points = domain.get_centroid_coordinates()
1553        answer = quantity.get_values(location='centroids')
1554
1555       
1556        #print quantity.get_values(points=interpolation_points)
1557        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
1558
1559
1560        #Arbitrary values
1561        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
1562                                    [1,4,-9],[2,5,0]])
1563
1564
1565        # Get interpolated values at centroids
1566        interpolation_points = domain.get_centroid_coordinates()
1567        answer = quantity.get_values(location='centroids')
1568        #print answer
1569        #print quantity.get_values(points=interpolation_points)
1570        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
1571                       
1572
1573        #FIXME TODO
1574        #indices = [0,5,3]
1575        #answer = [0.5,1,5]
1576        #assert allclose(answer,
1577        #                quantity.get_values(indices=indices, \
1578        #                                    location = 'unique vertices'))
1579
1580
1581
1582
1583    def test_get_interpolated_values_2(self):
1584        a = [0.0, 0.0]
1585        b = [0.0, 2.0]
1586        c = [2.0,0.0]
1587        d = [0.0, 4.0]
1588        e = [2.0, 2.0]
1589        f = [4.0,0.0]
1590
1591        points = [a, b, c, d, e, f]
1592        #bac, bce, ecf, dbe
1593        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1594
1595        domain = Domain(points, vertices)
1596
1597        quantity = Quantity(domain)
1598        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1599
1600        #First pick one point
1601        x, y = 2.0/3, 8.0/3
1602        v = quantity.get_values(interpolation_points = [[x,y]])
1603        assert allclose(v, 6)       
1604
1605        # Then another to test that algorithm won't blindly
1606        # reuse interpolation matrix
1607        x, y = 4.0/3, 4.0/3
1608        v = quantity.get_values(interpolation_points = [[x,y]])
1609        assert allclose(v, 4)       
1610
1611
1612
1613
1614    def test_getting_some_vertex_values(self):
1615        """
1616        get values based on triangle lists.
1617        """
1618        from mesh_factory import rectangular
1619        from shallow_water import Domain
1620        from Numeric import zeros, Float
1621
1622        #Create basic mesh
1623        points, vertices, boundary = rectangular(1, 3)
1624
1625        #print "points",points
1626        #print "vertices",vertices
1627        #print "boundary",boundary
1628
1629        #Create shallow water domain
1630        domain = Domain(points, vertices, boundary)
1631        #print "domain.number_of_elements ",domain.number_of_elements
1632        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1633                                    [4,4,4],[5,5,5],[6,6,6]])
1634        value = [7]
1635        indices = [1]
1636        quantity.set_values(value,
1637                            location = 'centroids',
1638                            indices = indices)
1639        #print "quantity.centroid_values",quantity.centroid_values
1640        #print "quantity.get_values(location = 'centroids') ",\
1641        #      quantity.get_values(location = 'centroids')
1642        assert allclose(quantity.centroid_values,
1643                        quantity.get_values(location = 'centroids'))
1644
1645
1646        value = [[15,20,25]]
1647        quantity.set_values(value, indices = indices)
1648        #print "1 quantity.vertex_values",quantity.vertex_values
1649        assert allclose(quantity.vertex_values, quantity.get_values())
1650
1651        assert allclose(quantity.edge_values,
1652                        quantity.get_values(location = 'edges'))
1653
1654        # get a subset of elements
1655        subset = quantity.get_values(location='centroids', indices=[0,5])
1656        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1657        assert allclose(subset, answer)
1658
1659
1660        subset = quantity.get_values(location='edges', indices=[0,5])
1661        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1662        #print "subset",subset
1663        #print "answer",answer
1664        assert allclose(subset, answer)
1665
1666        subset = quantity.get_values( indices=[1,5])
1667        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1668        #print "subset",subset
1669        #print "answer",answer
1670        assert allclose(subset, answer)
1671
1672
1673
1674
1675#-------------------------------------------------------------
1676if __name__ == "__main__":
1677    suite = unittest.makeSuite(Test_Quantity, 'test')
1678    #print "restricted test"
1679    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1680    runner = unittest.TextTestRunner()
1681    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.