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

Last change on this file since 5521 was 5521, checked in by ole, 14 years ago

Deprecated 'edges' as an option for location in set_values.
Validation tests pass as well as all unit tests except for two references to
this functionality.

File size: 82.8 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt, pi
5import tempfile
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
16from anuga.utilities.polygon import *
17
18#Aux for fit_interpolate.fit example
19def linear_function(point):
20    point = array(point)
21    return point[:,0]+point[:,1]
22
23
24class Test_Quantity(unittest.TestCase):
25    def setUp(self):
26        from domain import Domain
27
28        a = [0.0, 0.0]
29        b = [0.0, 2.0]
30        c = [2.0, 0.0]
31        d = [0.0, 4.0]
32        e = [2.0, 2.0]
33        f = [4.0, 0.0]
34
35        points = [a, b, c, d, e, f]
36
37        #bac, bce, ecf, dbe
38        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
39
40        self.mesh1 = Domain(points[:3], [elements[0]])
41        self.mesh1.check_integrity()
42
43        self.mesh4 = Domain(points, elements)
44        self.mesh4.check_integrity()
45
46        # UTM round Onslow
47        a = [240000, 7620000]
48        b = [240000, 7680000]
49        c = [300000, 7620000]
50
51        points = [a, b, c]
52        elements = [[0,2,1]]
53       
54        self.mesh_onslow = Domain(points, elements)
55        self.mesh_onslow.check_integrity()
56       
57    def tearDown(self):
58        pass
59        #print "  Tearing down"
60
61
62    def test_creation(self):
63
64        quantity = Quantity(self.mesh1, [[1,2,3]])
65        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
66
67        try:
68            quantity = Quantity()
69        except:
70            pass
71        else:
72            raise 'Should have raised empty quantity exception'
73
74
75        try:
76            quantity = Quantity([1,2,3])
77        except AssertionError:
78            pass
79        except:
80            raise 'Should have raised "mising mesh object" error'
81
82
83    def test_creation_zeros(self):
84
85        quantity = Quantity(self.mesh1)
86        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
87
88
89        quantity = Quantity(self.mesh4)
90        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
91                                                 [0.,0.,0.], [0.,0.,0.]])
92
93
94    def test_interpolation(self):
95        quantity = Quantity(self.mesh1, [[1,2,3]])
96        assert allclose(quantity.centroid_values, [2.0]) #Centroid
97
98        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
99
100
101    def test_interpolation2(self):
102        quantity = Quantity(self.mesh4,
103                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
104        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
105
106
107        quantity.extrapolate_second_order()
108
109        #print quantity.vertex_values
110        assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
111                                                 [3.+2./3, 6.+2./3, 4.+2./3],
112                                                 [4.6, 3.4, 1.],
113                                                 [-5.0, 1.0, 4.0]])
114
115        #print quantity.edge_values
116        assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
117                                               [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
118                                               [2.2, 2.8, 4.0],
119                                               [2.5, -0.5, -2.0]])
120
121
122        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
123        #                                       [5., 5., 5.],
124        #                                       [4.5, 4.5, 0.],
125        #                                       [3.0, -1.5, -1.5]])
126
127    def test_get_extrema_1(self):
128        quantity = Quantity(self.mesh4,
129                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
130        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
131
132        v = quantity.get_maximum_value()
133        assert v == 5
134
135        v = quantity.get_minimum_value()
136        assert v == 0       
137
138        i = quantity.get_maximum_index()
139        assert i == 1
140
141        i = quantity.get_minimum_index()
142        assert i == 3       
143       
144        x,y = quantity.get_maximum_location()
145        xref, yref = 4.0/3, 4.0/3
146        assert x == xref
147        assert y == yref
148
149        v = quantity.get_values(interpolation_points = [[x,y]])
150        assert allclose(v, 5)
151
152
153        x,y = quantity.get_minimum_location()
154        v = quantity.get_values(interpolation_points = [[x,y]])
155        assert allclose(v, 0)
156
157
158    def test_get_maximum_2(self):
159
160        a = [0.0, 0.0]
161        b = [0.0, 2.0]
162        c = [2.0,0.0]
163        d = [0.0, 4.0]
164        e = [2.0, 2.0]
165        f = [4.0,0.0]
166
167        points = [a, b, c, d, e, f]
168        #bac, bce, ecf, dbe
169        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
170
171        domain = Domain(points, vertices)
172
173        quantity = Quantity(domain)
174        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
175       
176        v = quantity.get_maximum_value()
177        assert v == 6
178
179        v = quantity.get_minimum_value()
180        assert v == 2       
181
182        i = quantity.get_maximum_index()
183        assert i == 3
184
185        i = quantity.get_minimum_index()
186        assert i == 0       
187       
188        x,y = quantity.get_maximum_location()
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        x,y = quantity.get_minimum_location()       
197        v = quantity.get_values(interpolation_points = [[x,y]])
198        assert allclose(v, 2)
199
200        #Multiple locations for maximum -
201        #Test that the algorithm picks the first occurrence       
202        v = quantity.get_maximum_value(indices=[0,1,2])
203        assert allclose(v, 4)
204
205        i = quantity.get_maximum_index(indices=[0,1,2])
206        assert i == 1
207       
208        x,y = quantity.get_maximum_location(indices=[0,1,2])
209        xref, yref = 4.0/3, 4.0/3
210        assert x == xref
211        assert y == yref
212
213        v = quantity.get_values(interpolation_points = [[x,y]])
214        assert allclose(v, 4)       
215
216        # More test of indices......
217        v = quantity.get_maximum_value(indices=[2,3])
218        assert allclose(v, 6)
219
220        i = quantity.get_maximum_index(indices=[2,3])
221        assert i == 3
222       
223        x,y = quantity.get_maximum_location(indices=[2,3])
224        xref, yref = 2.0/3, 8.0/3
225        assert x == xref
226        assert y == yref
227
228        v = quantity.get_values(interpolation_points = [[x,y]])
229        assert allclose(v, 6)       
230
231       
232
233    def test_boundary_allocation(self):
234        quantity = Quantity(self.mesh4,
235                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
236
237        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
238
239
240    def test_set_values(self):
241        quantity = Quantity(self.mesh4)
242
243
244        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
245                            location = 'vertices')
246        assert allclose(quantity.vertex_values,
247                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
248        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
249        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
250                                               [5., 5., 5.],
251                                               [4.5, 4.5, 0.],
252                                               [3.0, -1.5, -1.5]])
253
254
255        # Test default
256        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
257        assert allclose(quantity.vertex_values,
258                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
259        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
260        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
261                                               [5., 5., 5.],
262                                               [4.5, 4.5, 0.],
263                                               [3.0, -1.5, -1.5]])
264
265        # Test centroids
266        quantity.set_values([1,2,3,4], location = 'centroids')
267        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
268
269        # Test exceptions
270        try:
271            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
272                                location = 'bas kamel tuba')
273        except:
274            pass
275
276
277        try:
278            quantity.set_values([[1,2,3], [0,0,9]])
279        except AssertionError:
280            pass
281        except:
282            raise 'should have raised Assertionerror'
283
284
285
286    def test_set_values_const(self):
287        quantity = Quantity(self.mesh4)
288
289        quantity.set_values(1.0, location = 'vertices')
290        assert allclose(quantity.vertex_values,
291                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
292
293        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
294        assert allclose(quantity.edge_values, [[1, 1, 1],
295                                               [1, 1, 1],
296                                               [1, 1, 1],
297                                               [1, 1, 1]])
298
299
300        quantity.set_values(2.0, location = 'centroids')
301        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
302
303
304    def test_set_values_func(self):
305        quantity = Quantity(self.mesh4)
306
307        def f(x, y):
308            return x+y
309
310        quantity.set_values(f, location = 'vertices')
311        #print "quantity.vertex_values",quantity.vertex_values
312        assert allclose(quantity.vertex_values,
313                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
314        assert allclose(quantity.centroid_values,
315                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
316        assert allclose(quantity.edge_values,
317                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
318
319
320        quantity.set_values(f, location = 'centroids')
321        assert allclose(quantity.centroid_values,
322                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
323
324
325    def test_integral(self):
326        quantity = Quantity(self.mesh4)
327
328        #Try constants first
329        const = 5
330        quantity.set_values(const, location = 'vertices')
331        #print 'Q', quantity.get_integral()
332
333        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
334
335        #Try with a linear function
336        def f(x, y):
337            return x+y
338
339        quantity.set_values(f, location = 'vertices')
340
341
342        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
343
344        assert allclose (quantity.get_integral(), ref_integral)
345
346
347
348    def test_set_vertex_values(self):
349        quantity = Quantity(self.mesh4)
350        quantity.set_vertex_values([0,1,2,3,4,5])
351
352        assert allclose(quantity.vertex_values,
353                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
354        assert allclose(quantity.centroid_values,
355                        [1., 7./3, 11./3, 8./3]) #Centroid
356        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
357                                               [3., 2.5, 1.5],
358                                               [3.5, 4.5, 3.],
359                                               [2.5, 3.5, 2]])
360
361
362    def test_set_vertex_values_subset(self):
363        quantity = Quantity(self.mesh4)
364        quantity.set_vertex_values([0,1,2,3,4,5])
365        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
366
367        assert allclose(quantity.vertex_values,
368                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
369
370
371    def test_set_vertex_values_using_general_interface(self):
372        quantity = Quantity(self.mesh4)
373
374
375        quantity.set_values([0,1,2,3,4,5])
376
377
378        assert allclose(quantity.vertex_values,
379                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
380
381        #Centroid
382        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
383
384        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
385                                               [3., 2.5, 1.5],
386                                               [3.5, 4.5, 3.],
387                                               [2.5, 3.5, 2]])
388
389
390
391    def test_set_vertex_values_using_general_interface_with_subset(self):
392        """test_set_vertex_values_using_general_interface_with_subset(self):
393       
394        Test that indices and polygon works (for constants values)
395        """
396       
397        quantity = Quantity(self.mesh4)
398
399
400        quantity.set_values([0,2,3,5], indices=[0,2,3,5])
401        assert allclose(quantity.vertex_values,
402                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
403
404
405        # Constant
406        quantity.set_values(0.0)
407        quantity.set_values(3.14, indices=[0,2], location='vertices')
408
409        # Indices refer to triangle numbers
410        assert allclose(quantity.vertex_values,
411                        [[3.14,3.14,3.14], [0,0,0],
412                         [3.14,3.14,3.14], [0,0,0]])       
413       
414
415
416        # Now try with polygon (pick points where y>2)
417        polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
418        quantity.set_values(0.0)
419        quantity.set_values(3.14, polygon=polygon)
420       
421        assert allclose(quantity.vertex_values,
422                        [[0,0,0], [0,0,0], [0,0,0],
423                         [3.14,3.14,3.14]])               
424
425
426        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
427        # using centroids
428        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
429        quantity.set_values(0.0)
430        quantity.set_values(3.14, location='centroids', polygon=polygon)
431        assert allclose(quantity.vertex_values,
432                        [[0,0,0],
433                         [3.14,3.14,3.14],
434                         [3.14,3.14,3.14],                         
435                         [0,0,0]])               
436
437
438        # Same polygon now use vertices (default)
439        polygon = [[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]]
440        quantity.set_values(0.0)
441        #print 'Here 2'
442        quantity.set_values(3.14, polygon=polygon)
443        assert allclose(quantity.vertex_values,
444                        [[0,0,0],
445                         [3.14,3.14,3.14],
446                         [3.14,3.14,3.14],                         
447                         [0,0,0]])               
448       
449
450        # Test input checking
451        try:
452            quantity.set_values(3.14, polygon=polygon, indices = [0,2])
453        except:
454            pass
455        else:
456            msg = 'Should have caught this'
457            raise msg
458
459
460
461
462
463    def test_set_vertex_values_using_general_interface_subset_and_geo(self):
464        """test_set_vertex_values_using_general_interface_with_subset(self):
465        Test that indices and polygon works using georeferencing
466        """
467       
468        quantity = Quantity(self.mesh4)
469        G = Geo_reference(56, 10, 100)
470        quantity.domain.geo_reference = G
471
472        #print quantity.domain.get_nodes(absolute=True)
473
474
475        # Constant
476        quantity.set_values(0.0)
477        quantity.set_values(3.14, indices=[0,2], location='vertices')
478
479        # Indices refer to triangle numbers here - not vertices (why?)
480        assert allclose(quantity.vertex_values,
481                        [[3.14,3.14,3.14], [0,0,0],
482                         [3.14,3.14,3.14], [0,0,0]])       
483       
484
485
486        # Now try with polygon (pick points where y>2)
487        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
488        polygon += [G.xllcorner, G.yllcorner]
489       
490        quantity.set_values(0.0)
491        quantity.set_values(3.14, polygon=polygon, location='centroids')
492       
493        assert allclose(quantity.vertex_values,
494                        [[0,0,0], [0,0,0], [0,0,0],
495                         [3.14,3.14,3.14]])               
496
497
498        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
499        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
500        polygon += [G.xllcorner, G.yllcorner]
501       
502        quantity.set_values(0.0)
503        quantity.set_values(3.14, polygon=polygon)
504
505        assert allclose(quantity.vertex_values,
506                        [[0,0,0],
507                         [3.14,3.14,3.14],
508                         [3.14,3.14,3.14],                         
509                         [0,0,0]])               
510
511
512
513    def test_set_values_using_fit(self):
514
515
516        quantity = Quantity(self.mesh4)
517
518        #Get (enough) datapoints
519        data_points = [[ 0.66666667, 0.66666667],
520                       [ 1.33333333, 1.33333333],
521                       [ 2.66666667, 0.66666667],
522                       [ 0.66666667, 2.66666667],
523                       [ 0.0, 1.0],
524                       [ 0.0, 3.0],
525                       [ 1.0, 0.0],
526                       [ 1.0, 1.0],
527                       [ 1.0, 2.0],
528                       [ 1.0, 3.0],
529                       [ 2.0, 1.0],
530                       [ 3.0, 0.0],
531                       [ 3.0, 1.0]]
532
533        z = linear_function(data_points)
534
535        #Use built-in fit_interpolate.fit
536        quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
537        #quantity.set_values(points = data_points, values = z, alpha = 0)
538
539
540        answer = linear_function(quantity.domain.get_vertex_coordinates())
541        #print quantity.vertex_values, answer
542        assert allclose(quantity.vertex_values.flat, answer)
543
544
545        #Now try by setting the same values directly
546        vertex_attributes = fit_to_mesh(data_points,
547                                        quantity.domain.get_nodes(),
548                                        quantity.domain.triangles, #FIXME
549                                        point_attributes=z,
550                                        alpha = 0,
551                                        verbose=False)
552
553        #print vertex_attributes
554        quantity.set_values(vertex_attributes)
555        assert allclose(quantity.vertex_values.flat, answer)
556
557
558
559
560
561    def test_test_set_values_using_fit_w_geo(self):
562
563
564        #Mesh
565        vertex_coordinates = [[0.76, 0.76],
566                              [0.76, 5.76],
567                              [5.76, 0.76]]
568        triangles = [[0,2,1]]
569
570        mesh_georef = Geo_reference(56,-0.76,-0.76)
571        mesh1 = Domain(vertex_coordinates, triangles,
572                       geo_reference = mesh_georef)
573        mesh1.check_integrity()
574
575        #Quantity
576        quantity = Quantity(mesh1)
577
578        #Data
579        data_points = [[ 201.0, 401.0],
580                       [ 201.0, 403.0],
581                       [ 203.0, 401.0]]
582
583        z = [2, 4, 4]
584
585        data_georef = Geo_reference(56,-200,-400)
586
587
588        #Reference
589        ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
590                          point_attributes=z,
591                          data_origin = data_georef.get_origin(),
592                          mesh_origin = mesh_georef.get_origin(),
593                          alpha = 0)
594
595        assert allclose( ref, [0,5,5] )
596
597
598        #Test set_values
599
600        quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
601
602        #quantity.set_values(points = data_points,
603        #                    values = z,
604        #                    data_georef = data_georef,
605        #                    alpha = 0)
606
607
608        #quantity.set_values(points = data_points,
609        #                    values = z,
610        #                    data_georef = data_georef,
611        #                    alpha = 0)
612        assert allclose(quantity.vertex_values.flat, ref)
613
614
615
616        #Test set_values using geospatial data object
617        quantity.vertex_values[:] = 0.0
618
619        geo = Geospatial_data(data_points, z, data_georef)
620
621
622        quantity.set_values(geospatial_data = geo, alpha = 0)
623        assert allclose(quantity.vertex_values.flat, ref)
624
625
626
627    def test_set_values_from_file1(self):
628        quantity = Quantity(self.mesh4)
629
630        #Get (enough) datapoints
631        data_points = [[ 0.66666667, 0.66666667],
632                       [ 1.33333333, 1.33333333],
633                       [ 2.66666667, 0.66666667],
634                       [ 0.66666667, 2.66666667],
635                       [ 0.0, 1.0],
636                       [ 0.0, 3.0],
637                       [ 1.0, 0.0],
638                       [ 1.0, 1.0],
639                       [ 1.0, 2.0],
640                       [ 1.0, 3.0],
641                       [ 2.0, 1.0],
642                       [ 3.0, 0.0],
643                       [ 3.0, 1.0]]
644
645        data_geo_spatial = Geospatial_data(data_points,
646                         geo_reference = Geo_reference(56, 0, 0))
647        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
648        attributes = linear_function(data_points_absolute)
649        att = 'spam_and_eggs'
650       
651        #Create .txt file
652        ptsfile = tempfile.mktemp(".txt")
653        file = open(ptsfile,"w")
654        file.write(" x,y," + att + " \n")
655        for data_point, attribute in map(None, data_points_absolute
656                                         ,attributes):
657            row = str(data_point[0]) + ',' + str(data_point[1]) \
658                  + ',' + str(attribute)
659            file.write(row + "\n")
660        file.close()
661
662
663        #Check that values can be set from file
664        quantity.set_values(filename = ptsfile,
665                            attribute_name = att, alpha = 0)
666        answer = linear_function(quantity.domain.get_vertex_coordinates())
667
668        #print quantity.vertex_values.flat
669        #print answer
670
671
672        assert allclose(quantity.vertex_values.flat, answer)
673
674
675        #Check that values can be set from file using default attribute
676        quantity.set_values(filename = ptsfile, alpha = 0)
677        assert allclose(quantity.vertex_values.flat, answer)
678
679        #Cleanup
680        import os
681        os.remove(ptsfile)
682
683
684
685    def Xtest_set_values_from_file_using_polygon(self):
686        """test_set_values_from_file_using_polygon(self):
687       
688        Test that polygon restriction works for general points data
689        """
690       
691        quantity = Quantity(self.mesh4)
692
693        #Get (enough) datapoints
694        data_points = [[ 0.66666667, 0.66666667],
695                       [ 1.33333333, 1.33333333],
696                       [ 2.66666667, 0.66666667],
697                       [ 0.66666667, 2.66666667],
698                       [ 0.0, 1.0],
699                       [ 0.0, 3.0],
700                       [ 1.0, 0.0],
701                       [ 1.0, 1.0],
702                       [ 1.0, 2.0],
703                       [ 1.0, 3.0],
704                       [ 2.0, 1.0],
705                       [ 3.0, 0.0],
706                       [ 3.0, 1.0]]
707
708        data_geo_spatial = Geospatial_data(data_points,
709                         geo_reference = Geo_reference(56, 0, 0))
710        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
711        attributes = linear_function(data_points_absolute)
712        att = 'spam_and_eggs'
713       
714        #Create .txt file
715        ptsfile = tempfile.mktemp(".txt")
716        file = open(ptsfile,"w")
717        file.write(" x,y," + att + " \n")
718        for data_point, attribute in map(None, data_points_absolute
719                                         ,attributes):
720            row = str(data_point[0]) + ',' + str(data_point[1]) \
721                  + ',' + str(attribute)
722            file.write(row + "\n")
723        file.close()
724
725        # Create restricting polygon (containing node #4 (2,2) and
726        # centroid of triangle #1 (bce)
727        polygon = [[1.0, 1.0], [4.0, 1.0],
728                   [4.0, 4.0], [1.0, 4.0]]
729
730        #print self.mesh4.nodes
731        #print inside_polygon(self.mesh4.nodes, polygon)
732        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)
733
734        #print quantity.domain.get_vertex_coordinates()
735        #print quantity.domain.get_nodes()       
736       
737        # Check that values can be set from file
738        quantity.set_values(filename=ptsfile,
739                            polygon=polygon,
740                            location='unique vertices',
741                            alpha=0)
742
743        # Get indices for vertex coordinates in polygon
744        indices = inside_polygon(quantity.domain.get_vertex_coordinates(), 
745                                 polygon)
746        points = take(quantity.domain.get_vertex_coordinates(), indices)
747       
748        answer = linear_function(points)
749
750        #print quantity.vertex_values.flat
751        #print answer
752
753        # Check vertices in polygon have been set
754        assert allclose(take(quantity.vertex_values.flat, indices),
755                             answer)
756
757        # Check vertices outside polygon are zero
758        indices = outside_polygon(quantity.domain.get_vertex_coordinates(), 
759                                  polygon)       
760        assert allclose(take(quantity.vertex_values.flat, indices),
761                             0.0)       
762
763        #Cleanup
764        import os
765        os.remove(ptsfile)
766
767
768       
769
770    def test_cache_test_set_values_from_file(self):
771        # FIXME (Ole): What is this about?
772        # I don't think it checks anything new
773        quantity = Quantity(self.mesh4)
774
775        #Get (enough) datapoints
776        data_points = [[ 0.66666667, 0.66666667],
777                       [ 1.33333333, 1.33333333],
778                       [ 2.66666667, 0.66666667],
779                       [ 0.66666667, 2.66666667],
780                       [ 0.0, 1.0],
781                       [ 0.0, 3.0],
782                       [ 1.0, 0.0],
783                       [ 1.0, 1.0],
784                       [ 1.0, 2.0],
785                       [ 1.0, 3.0],
786                       [ 2.0, 1.0],
787                       [ 3.0, 0.0],
788                       [ 3.0, 1.0]]
789
790        georef = Geo_reference(56, 0, 0)
791        data_geo_spatial = Geospatial_data(data_points,
792                                           geo_reference=georef)
793                                           
794        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
795        attributes = linear_function(data_points_absolute)
796        att = 'spam_and_eggs'
797       
798        # Create .txt file
799        ptsfile = tempfile.mktemp(".txt")
800        file = open(ptsfile,"w")
801        file.write(" x,y," + att + " \n")
802        for data_point, attribute in map(None, data_points_absolute
803                                         ,attributes):
804            row = str(data_point[0]) + ',' + str(data_point[1]) \
805                  + ',' + str(attribute)
806            file.write(row + "\n")
807        file.close()
808
809
810        # Check that values can be set from file
811        quantity.set_values(filename=ptsfile,
812                            attribute_name=att, 
813                            alpha=0, 
814                            use_cache=True,
815                            verbose=False)
816        answer = linear_function(quantity.domain.get_vertex_coordinates())
817        assert allclose(quantity.vertex_values.flat, answer)
818
819
820        # Check that values can be set from file using default attribute
821        quantity.set_values(filename=ptsfile, 
822                            alpha=0)
823        assert allclose(quantity.vertex_values.flat, answer)
824
825        # Check cache
826        quantity.set_values(filename=ptsfile,
827                            attribute_name=att, 
828                            alpha=0, 
829                            use_cache=True,
830                            verbose=False)
831       
832       
833        #Cleanup
834        import os
835        os.remove(ptsfile)
836
837    def test_set_values_from_lat_long(self):
838        quantity = Quantity(self.mesh_onslow)
839
840        #Get (enough) datapoints
841        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
842                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
843
844        data_geo_spatial = Geospatial_data(data_points,
845                                           points_are_lats_longs=True)
846        points_UTM = data_geo_spatial.get_data_points(absolute=True)
847        attributes = linear_function(points_UTM)
848        att = 'elevation'
849       
850        #Create .txt file
851        txt_file = tempfile.mktemp(".txt")
852        file = open(txt_file,"w")
853        file.write(" lat,long," + att + " \n")
854        for data_point, attribute in map(None, data_points, attributes):
855            row = str(data_point[0]) + ',' + str(data_point[1]) \
856                  + ',' + str(attribute)
857            #print "row", row
858            file.write(row + "\n")
859        file.close()
860
861
862        #Check that values can be set from file
863        quantity.set_values(filename=txt_file,
864                            attribute_name=att, 
865                            alpha=0)
866        answer = linear_function(quantity.domain.get_vertex_coordinates())
867
868        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
869        #print "answer",answer
870
871        assert allclose(quantity.vertex_values.flat, answer)
872
873
874        #Check that values can be set from file using default attribute
875        quantity.set_values(filename=txt_file, alpha=0)
876        assert allclose(quantity.vertex_values.flat, answer)
877
878        #Cleanup
879        import os
880        os.remove(txt_file)
881         
882    def test_set_values_from_lat_long(self):
883        quantity = Quantity(self.mesh_onslow)
884
885        #Get (enough) datapoints
886        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
887                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
888
889        data_geo_spatial = Geospatial_data(data_points,
890                                           points_are_lats_longs=True)
891        points_UTM = data_geo_spatial.get_data_points(absolute=True)
892        attributes = linear_function(points_UTM)
893        att = 'elevation'
894       
895        #Create .txt file
896        txt_file = tempfile.mktemp(".txt")
897        file = open(txt_file,"w")
898        file.write(" lat,long," + att + " \n")
899        for data_point, attribute in map(None, data_points, attributes):
900            row = str(data_point[0]) + ',' + str(data_point[1]) \
901                  + ',' + str(attribute)
902            #print "row", row
903            file.write(row + "\n")
904        file.close()
905
906
907        #Check that values can be set from file
908        quantity.set_values(filename=txt_file,
909                            attribute_name=att, alpha=0)
910        answer = linear_function(quantity.domain.get_vertex_coordinates())
911
912        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
913        #print "answer",answer
914
915        assert allclose(quantity.vertex_values.flat, answer)
916
917
918        #Check that values can be set from file using default attribute
919        quantity.set_values(filename=txt_file, alpha=0)
920        assert allclose(quantity.vertex_values.flat, answer)
921
922        #Cleanup
923        import os
924        os.remove(txt_file)
925       
926    def test_set_values_from_UTM_pts(self):
927        quantity = Quantity(self.mesh_onslow)
928
929        #Get (enough) datapoints
930        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
931                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6]]
932
933        data_geo_spatial = Geospatial_data(data_points,
934                                           points_are_lats_longs=True)
935        points_UTM = data_geo_spatial.get_data_points(absolute=True)
936        attributes = linear_function(points_UTM)
937        att = 'elevation'
938       
939        #Create .txt file
940        txt_file = tempfile.mktemp(".txt")
941        file = open(txt_file,"w")
942        file.write(" x,y," + att + " \n")
943        for data_point, attribute in map(None, points_UTM, attributes):
944            row = str(data_point[0]) + ',' + str(data_point[1]) \
945                  + ',' + str(attribute)
946            #print "row", row
947            file.write(row + "\n")
948        file.close()
949
950
951        pts_file = tempfile.mktemp(".pts")       
952        convert = Geospatial_data(txt_file)
953        convert.export_points_file(pts_file)
954
955        #Check that values can be set from file
956        quantity.set_values_from_file(pts_file, att, 0,
957                                      'vertices', None)
958        answer = linear_function(quantity.domain.get_vertex_coordinates())
959        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
960        #print "answer",answer
961        assert allclose(quantity.vertex_values.flat, answer)
962
963        #Check that values can be set from file
964        quantity.set_values(filename=pts_file,
965                            attribute_name=att, alpha=0)
966        answer = linear_function(quantity.domain.get_vertex_coordinates())
967        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
968        #print "answer",answer
969        assert allclose(quantity.vertex_values.flat, answer)
970
971
972        #Check that values can be set from file using default attribute
973        quantity.set_values(filename=txt_file, alpha=0)
974        assert allclose(quantity.vertex_values.flat, answer)
975
976        #Cleanup
977        import os
978        os.remove(txt_file)
979        os.remove(pts_file)
980       
981    def verbose_test_set_values_from_UTM_pts(self):
982        quantity = Quantity(self.mesh_onslow)
983
984        #Get (enough) datapoints
985        data_points = [[-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
986                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
987                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
988                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
989                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
990                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
991                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
992                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
993                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
994                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
995                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
996                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
997                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
998                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
999                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1000                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1001                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1002                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1003                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1004                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1005                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1006                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1007                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1008                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1009                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1010                       [-21.5, 114.5],[-21.4, 114.6],[-21.45,114.65],
1011                       [-21.35, 114.65],[-21.45, 114.55],[-21.45,114.6],
1012                       ]
1013
1014        data_geo_spatial = Geospatial_data(data_points,
1015                                           points_are_lats_longs=True)
1016        points_UTM = data_geo_spatial.get_data_points(absolute=True)
1017        attributes = linear_function(points_UTM)
1018        att = 'elevation'
1019       
1020        #Create .txt file
1021        txt_file = tempfile.mktemp(".txt")
1022        file = open(txt_file,"w")
1023        file.write(" x,y," + att + " \n")
1024        for data_point, attribute in map(None, points_UTM, attributes):
1025            row = str(data_point[0]) + ',' + str(data_point[1]) \
1026                  + ',' + str(attribute)
1027            #print "row", row
1028            file.write(row + "\n")
1029        file.close()
1030
1031
1032        pts_file = tempfile.mktemp(".pts")       
1033        convert = Geospatial_data(txt_file)
1034        convert.export_points_file(pts_file)
1035
1036        #Check that values can be set from file
1037        quantity.set_values_from_file(pts_file, att, 0,
1038                                      'vertices', None, verbose = True,
1039                                      max_read_lines=2)
1040        answer = linear_function(quantity.domain.get_vertex_coordinates())
1041        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
1042        #print "answer",answer
1043        assert allclose(quantity.vertex_values.flat, answer)
1044
1045        #Check that values can be set from file
1046        quantity.set_values(filename=pts_file,
1047                            attribute_name=att, alpha=0)
1048        answer = linear_function(quantity.domain.get_vertex_coordinates())
1049        #print "quantity.vertex_values.flat", quantity.vertex_values.flat
1050        #print "answer",answer
1051        assert allclose(quantity.vertex_values.flat, answer)
1052
1053
1054        #Check that values can be set from file using default attribute
1055        quantity.set_values(filename=txt_file, alpha=0)
1056        assert allclose(quantity.vertex_values.flat, answer)
1057
1058        #Cleanup
1059        import os
1060        os.remove(txt_file)
1061        os.remove(pts_file)
1062       
1063    def test_set_values_from_file_with_georef1(self):
1064
1065        #Mesh in zone 56 (absolute coords)
1066
1067        x0 = 314036.58727982
1068        y0 = 6224951.2960092
1069
1070        a = [x0+0.0, y0+0.0]
1071        b = [x0+0.0, y0+2.0]
1072        c = [x0+2.0, y0+0.0]
1073        d = [x0+0.0, y0+4.0]
1074        e = [x0+2.0, y0+2.0]
1075        f = [x0+4.0, y0+0.0]
1076
1077        points = [a, b, c, d, e, f]
1078
1079        #bac, bce, ecf, dbe
1080        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1081
1082        #absolute going in ..
1083        mesh4 = Domain(points, elements,
1084                       geo_reference = Geo_reference(56, 0, 0))
1085        mesh4.check_integrity()
1086        quantity = Quantity(mesh4)
1087
1088        #Get (enough) datapoints (relative to georef)
1089        data_points_rel = [[ 0.66666667, 0.66666667],
1090                       [ 1.33333333, 1.33333333],
1091                       [ 2.66666667, 0.66666667],
1092                       [ 0.66666667, 2.66666667],
1093                       [ 0.0, 1.0],
1094                       [ 0.0, 3.0],
1095                       [ 1.0, 0.0],
1096                       [ 1.0, 1.0],
1097                       [ 1.0, 2.0],
1098                       [ 1.0, 3.0],
1099                       [ 2.0, 1.0],
1100                       [ 3.0, 0.0],
1101                       [ 3.0, 1.0]]
1102
1103        data_geo_spatial = Geospatial_data(data_points_rel,
1104                         geo_reference = Geo_reference(56, x0, y0))
1105        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1106        attributes = linear_function(data_points_absolute)
1107        att = 'spam_and_eggs'
1108       
1109        #Create .txt file
1110        ptsfile = tempfile.mktemp(".txt")
1111        file = open(ptsfile,"w")
1112        file.write(" x,y," + att + " \n")
1113        for data_point, attribute in map(None, data_points_absolute
1114                                         ,attributes):
1115            row = str(data_point[0]) + ',' + str(data_point[1]) \
1116                  + ',' + str(attribute)
1117            file.write(row + "\n")
1118        file.close()
1119
1120        #file = open(ptsfile, 'r')
1121        #lines = file.readlines()
1122        #file.close()
1123     
1124
1125        #Check that values can be set from file
1126        quantity.set_values(filename=ptsfile,
1127                            attribute_name=att, alpha=0)
1128        answer = linear_function(quantity.domain.get_vertex_coordinates())
1129
1130        assert allclose(quantity.vertex_values.flat, answer)
1131
1132
1133        #Check that values can be set from file using default attribute
1134        quantity.set_values(filename=ptsfile, alpha=0)
1135        assert allclose(quantity.vertex_values.flat, answer)
1136
1137        #Cleanup
1138        import os
1139        os.remove(ptsfile)
1140
1141
1142    def test_set_values_from_file_with_georef2(self):
1143
1144        #Mesh in zone 56 (relative coords)
1145
1146        x0 = 314036.58727982
1147        y0 = 6224951.2960092
1148        #x0 = 0.0
1149        #y0 = 0.0
1150
1151        a = [0.0, 0.0]
1152        b = [0.0, 2.0]
1153        c = [2.0, 0.0]
1154        d = [0.0, 4.0]
1155        e = [2.0, 2.0]
1156        f = [4.0, 0.0]
1157
1158        points = [a, b, c, d, e, f]
1159
1160        #bac, bce, ecf, dbe
1161        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
1162
1163        mesh4 = Domain(points, elements,
1164                       geo_reference = Geo_reference(56, x0, y0))
1165        mesh4.check_integrity()
1166        quantity = Quantity(mesh4)
1167
1168        #Get (enough) datapoints
1169        data_points = [[ x0+0.66666667, y0+0.66666667],
1170                       [ x0+1.33333333, y0+1.33333333],
1171                       [ x0+2.66666667, y0+0.66666667],
1172                       [ x0+0.66666667, y0+2.66666667],
1173                       [ x0+0.0, y0+1.0],
1174                       [ x0+0.0, y0+3.0],
1175                       [ x0+1.0, y0+0.0],
1176                       [ x0+1.0, y0+1.0],
1177                       [ x0+1.0, y0+2.0],
1178                       [ x0+1.0, y0+3.0],
1179                       [ x0+2.0, y0+1.0],
1180                       [ x0+3.0, y0+0.0],
1181                       [ x0+3.0, y0+1.0]]
1182
1183
1184        data_geo_spatial = Geospatial_data(data_points,
1185                         geo_reference = Geo_reference(56, 0, 0))
1186        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
1187        attributes = linear_function(data_points_absolute)
1188        att = 'spam_and_eggs'
1189       
1190        #Create .txt file
1191        ptsfile = tempfile.mktemp(".txt")
1192        file = open(ptsfile,"w")
1193        file.write(" x,y," + att + " \n")
1194        for data_point, attribute in map(None, data_points_absolute
1195                                         ,attributes):
1196            row = str(data_point[0]) + ',' + str(data_point[1]) \
1197                  + ',' + str(attribute)
1198            file.write(row + "\n")
1199        file.close()
1200
1201
1202        #Check that values can be set from file
1203        quantity.set_values(filename=ptsfile,
1204                            attribute_name=att, alpha=0)
1205        answer = linear_function(quantity.domain. \
1206                                 get_vertex_coordinates(absolute=True))
1207
1208
1209        assert allclose(quantity.vertex_values.flat, answer)
1210
1211
1212        #Check that values can be set from file using default attribute
1213        quantity.set_values(filename=ptsfile, alpha=0)
1214        assert allclose(quantity.vertex_values.flat, answer)
1215
1216        #Cleanup
1217        import os
1218        os.remove(ptsfile)
1219
1220
1221
1222
1223    def test_set_values_from_quantity(self):
1224
1225        quantity1 = Quantity(self.mesh4)
1226        quantity1.set_vertex_values([0,1,2,3,4,5])
1227
1228        assert allclose(quantity1.vertex_values,
1229                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1230
1231
1232        quantity2 = Quantity(self.mesh4)
1233        quantity2.set_values(quantity=quantity1)
1234        assert allclose(quantity2.vertex_values,
1235                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1236
1237        quantity2.set_values(quantity = 2*quantity1)
1238        assert allclose(quantity2.vertex_values,
1239                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
1240
1241        quantity2.set_values(quantity = 2*quantity1 + 3)
1242        assert allclose(quantity2.vertex_values,
1243                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
1244
1245
1246        #Check detection of quantity as first orgument
1247        quantity2.set_values(2*quantity1 + 3)
1248        assert allclose(quantity2.vertex_values,
1249                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
1250
1251
1252
1253    def Xtest_set_values_from_quantity_using_polygon(self):
1254        """test_set_values_from_quantity_using_polygon(self):
1255       
1256        Check that polygon can be used to restrict set_values when
1257        using another quantity as argument.
1258        """
1259       
1260        # Create restricting polygon (containing node #4 (2,2) and
1261        # centroid of triangle #1 (bce)
1262        polygon = [[1.0, 1.0], [4.0, 1.0],
1263                   [4.0, 4.0], [1.0, 4.0]]
1264        assert allclose(inside_polygon(self.mesh4.nodes, polygon), 4)                   
1265       
1266        quantity1 = Quantity(self.mesh4)
1267        quantity1.set_vertex_values([0,1,2,3,4,5])
1268
1269        assert allclose(quantity1.vertex_values,
1270                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1271
1272
1273        quantity2 = Quantity(self.mesh4)
1274        quantity2.set_values(quantity=quantity1,
1275                             polygon=polygon)
1276                             
1277        msg = 'Only node #4(e) at (2,2) should have values applied '
1278        assert allclose(quantity2.vertex_values,
1279                        [[0,0,0], [0,0,4], [4,0,0], [0,0,4]]), msg       
1280                        #bac,     bce,     ecf,     dbe
1281                       
1282
1283
1284    def test_overloading(self):
1285
1286        quantity1 = Quantity(self.mesh4)
1287        quantity1.set_vertex_values([0,1,2,3,4,5])
1288
1289        assert allclose(quantity1.vertex_values,
1290                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1291
1292
1293        quantity2 = Quantity(self.mesh4)
1294        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
1295                             location = 'vertices')
1296
1297
1298
1299        quantity3 = Quantity(self.mesh4)
1300        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
1301                             location = 'vertices')
1302
1303
1304        # Negation
1305        Q = -quantity1
1306        assert allclose(Q.vertex_values, -quantity1.vertex_values)
1307        assert allclose(Q.centroid_values, -quantity1.centroid_values)
1308        assert allclose(Q.edge_values, -quantity1.edge_values)
1309
1310        # Addition
1311        Q = quantity1 + 7
1312        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
1313        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
1314        assert allclose(Q.edge_values, quantity1.edge_values + 7)
1315
1316        Q = 7 + quantity1
1317        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
1318        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
1319        assert allclose(Q.edge_values, quantity1.edge_values + 7)
1320
1321        Q = quantity1 + quantity2
1322        assert allclose(Q.vertex_values,
1323                        quantity1.vertex_values + quantity2.vertex_values)
1324        assert allclose(Q.centroid_values,
1325                        quantity1.centroid_values + quantity2.centroid_values)
1326        assert allclose(Q.edge_values,
1327                        quantity1.edge_values + quantity2.edge_values)
1328
1329
1330        Q = quantity1 + quantity2 - 3
1331        assert allclose(Q.vertex_values,
1332                        quantity1.vertex_values + quantity2.vertex_values - 3)
1333
1334        Q = quantity1 - quantity2
1335        assert allclose(Q.vertex_values,
1336                        quantity1.vertex_values - quantity2.vertex_values)
1337
1338        #Scaling
1339        Q = quantity1*3
1340        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
1341        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
1342        assert allclose(Q.edge_values, quantity1.edge_values*3)
1343        Q = 3*quantity1
1344        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
1345
1346        #Multiplication
1347        Q = quantity1 * quantity2
1348        #print Q.vertex_values
1349        #print Q.centroid_values
1350        #print quantity1.centroid_values
1351        #print quantity2.centroid_values
1352
1353        assert allclose(Q.vertex_values,
1354                        quantity1.vertex_values * quantity2.vertex_values)
1355
1356        #Linear combinations
1357        Q = 4*quantity1 + 2
1358        assert allclose(Q.vertex_values,
1359                        4*quantity1.vertex_values + 2)
1360
1361        Q = quantity1*quantity2 + 2
1362        assert allclose(Q.vertex_values,
1363                        quantity1.vertex_values * quantity2.vertex_values + 2)
1364
1365        Q = quantity1*quantity2 + quantity3
1366        assert allclose(Q.vertex_values,
1367                        quantity1.vertex_values * quantity2.vertex_values +
1368                        quantity3.vertex_values)
1369        Q = quantity1*quantity2 + 3*quantity3
1370        assert allclose(Q.vertex_values,
1371                        quantity1.vertex_values * quantity2.vertex_values +
1372                        3*quantity3.vertex_values)
1373        Q = quantity1*quantity2 + 3*quantity3 + 5.0
1374        assert allclose(Q.vertex_values,
1375                        quantity1.vertex_values * quantity2.vertex_values +
1376                        3*quantity3.vertex_values + 5)
1377
1378        Q = quantity1*quantity2 - quantity3
1379        assert allclose(Q.vertex_values,
1380                        quantity1.vertex_values * quantity2.vertex_values -
1381                        quantity3.vertex_values)
1382        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
1383        assert allclose(Q.vertex_values,
1384                        1.5*quantity1.vertex_values * quantity2.vertex_values -
1385                        3*quantity3.vertex_values + 5)
1386
1387        #Try combining quantities and arrays and scalars
1388        Q = 1.5*quantity1*quantity2.vertex_values -\
1389            3*quantity3.vertex_values + 5.0
1390        assert allclose(Q.vertex_values,
1391                        1.5*quantity1.vertex_values * quantity2.vertex_values -
1392                        3*quantity3.vertex_values + 5)
1393
1394
1395        #Powers
1396        Q = quantity1**2
1397        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
1398
1399        Q = quantity1**2 +quantity2**2
1400        assert allclose(Q.vertex_values,
1401                        quantity1.vertex_values**2 + \
1402                        quantity2.vertex_values**2)
1403
1404        Q = (quantity1**2 +quantity2**2)**0.5
1405        assert allclose(Q.vertex_values,
1406                        (quantity1.vertex_values**2 + \
1407                        quantity2.vertex_values**2)**0.5)
1408
1409
1410
1411
1412
1413
1414
1415    def test_compute_gradient(self):
1416        quantity = Quantity(self.mesh4)
1417
1418        #Set up for a gradient of (2,0) at mid triangle
1419        quantity.set_values([2.0, 4.0, 6.0, 2.0],
1420                            location = 'centroids')
1421
1422
1423        #Gradients
1424        quantity.compute_gradients()
1425
1426        a = quantity.x_gradient
1427        b = quantity.y_gradient
1428        #print self.mesh4.centroid_coordinates
1429        #print a, b
1430
1431        #The central triangle (1)
1432        #(using standard gradient based on neigbours controid values)
1433        assert allclose(a[1], 2.0)
1434        assert allclose(b[1], 0.0)
1435
1436
1437        #Left triangle (0) using two point gradient
1438        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
1439        #2  = 4  + a*(-2/3)  + b*(-2/3)
1440        assert allclose(a[0] + b[0], 3)
1441        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
1442        assert allclose(a[0] - b[0], 0)
1443
1444
1445        #Right triangle (2) using two point gradient
1446        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
1447        #6  = 4  + a*(4/3)  + b*(-2/3)
1448        assert allclose(2*a[2] - b[2], 3)
1449        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
1450        assert allclose(a[2] + 2*b[2], 0)
1451
1452
1453        #Top triangle (3) using two point gradient
1454        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
1455        #2  = 4  + a*(-2/3)  + b*(4/3)
1456        assert allclose(a[3] - 2*b[3], 3)
1457        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
1458        assert allclose(2*a[3] + b[3], 0)
1459
1460
1461
1462        #print a, b
1463        quantity.extrapolate_second_order()
1464
1465        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
1466        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
1467        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
1468
1469
1470        #a = 1.2, b=-0.6
1471        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
1472        assert allclose(quantity.vertex_values[2,2], 8)
1473
1474    def test_get_gradients(self):
1475        quantity = Quantity(self.mesh4)
1476
1477        #Set up for a gradient of (2,0) at mid triangle
1478        quantity.set_values([2.0, 4.0, 6.0, 2.0],
1479                            location = 'centroids')
1480
1481
1482        #Gradients
1483        quantity.compute_gradients()
1484
1485        a, b = quantity.get_gradients()
1486        #print self.mesh4.centroid_coordinates
1487        #print a, b
1488
1489        #The central triangle (1)
1490        #(using standard gradient based on neigbours controid values)
1491        assert allclose(a[1], 2.0)
1492        assert allclose(b[1], 0.0)
1493
1494
1495        #Left triangle (0) using two point gradient
1496        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
1497        #2  = 4  + a*(-2/3)  + b*(-2/3)
1498        assert allclose(a[0] + b[0], 3)
1499        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
1500        assert allclose(a[0] - b[0], 0)
1501
1502
1503        #Right triangle (2) using two point gradient
1504        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
1505        #6  = 4  + a*(4/3)  + b*(-2/3)
1506        assert allclose(2*a[2] - b[2], 3)
1507        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
1508        assert allclose(a[2] + 2*b[2], 0)
1509
1510
1511        #Top triangle (3) using two point gradient
1512        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
1513        #2  = 4  + a*(-2/3)  + b*(4/3)
1514        assert allclose(a[3] - 2*b[3], 3)
1515        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
1516        assert allclose(2*a[3] + b[3], 0)
1517
1518
1519    def test_second_order_extrapolation2(self):
1520        quantity = Quantity(self.mesh4)
1521
1522        #Set up for a gradient of (3,1), f(x) = 3x+y
1523        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
1524                            location = 'centroids')
1525
1526        #Gradients
1527        quantity.compute_gradients()
1528
1529        a = quantity.x_gradient
1530        b = quantity.y_gradient
1531       
1532        #print a, b
1533
1534        assert allclose(a[1], 3.0)
1535        assert allclose(b[1], 1.0)
1536
1537        #Work out the others
1538
1539        quantity.extrapolate_second_order()
1540
1541        #print quantity.vertex_values
1542        assert allclose(quantity.vertex_values[1,0], 2.0)
1543        assert allclose(quantity.vertex_values[1,1], 6.0)
1544        assert allclose(quantity.vertex_values[1,2], 8.0)
1545
1546
1547
1548    def test_backup_saxpy_centroid_values(self):
1549        quantity = Quantity(self.mesh4)
1550
1551        #Set up for a gradient of (3,1), f(x) = 3x+y
1552        c_values = array([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3])
1553        d_values = array([1.0, 2.0, 3.0, 4.0])
1554        quantity.set_values(c_values, location = 'centroids')
1555
1556        #Backup
1557        quantity.backup_centroid_values()
1558
1559        #print quantity.vertex_values
1560        assert allclose(quantity.centroid_values, quantity.centroid_backup_values)
1561
1562
1563        quantity.set_values(d_values, location = 'centroids')
1564
1565        quantity.saxpy_centroid_values(2.0, 3.0)
1566
1567        assert(quantity.centroid_values, 2.0*d_values + 3.0*c_values)
1568
1569
1570
1571    def test_first_order_extrapolator(self):
1572        quantity = Quantity(self.mesh4)
1573
1574        #Test centroids
1575        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1576        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1577
1578        #Extrapolate
1579        quantity.extrapolate_first_order()
1580
1581        #Check that gradient is zero
1582        a,b = quantity.get_gradients()
1583        assert allclose(a, [0,0,0,0])
1584        assert allclose(b, [0,0,0,0])
1585
1586        #Check vertices but not edge values
1587        assert allclose(quantity.vertex_values,
1588                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1589
1590
1591    def test_second_order_extrapolator(self):
1592        quantity = Quantity(self.mesh4)
1593
1594        #Set up for a gradient of (3,0) at mid triangle
1595        quantity.set_values([2.0, 4.0, 8.0, 2.0],
1596                            location = 'centroids')
1597
1598
1599
1600        quantity.extrapolate_second_order()
1601        quantity.limit()
1602
1603
1604        #Assert that central triangle is limited by neighbours
1605        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
1606        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
1607
1608        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
1609        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
1610
1611        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
1612        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
1613
1614
1615        #Assert that quantities are conserved
1616        from Numeric import sum
1617        for k in range(quantity.centroid_values.shape[0]):
1618            assert allclose (quantity.centroid_values[k],
1619                             sum(quantity.vertex_values[k,:])/3)
1620
1621
1622
1623
1624
1625    def test_limit_vertices_by_all_neighbours(self):
1626        quantity = Quantity(self.mesh4)
1627
1628        #Create a deliberate overshoot (e.g. from gradient computation)
1629        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1630
1631
1632        #Limit
1633        quantity.limit_vertices_by_all_neighbours()
1634
1635        #Assert that central triangle is limited by neighbours
1636        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
1637        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
1638
1639        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
1640        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
1641
1642        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
1643        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
1644
1645
1646
1647        #Assert that quantities are conserved
1648        from Numeric import sum
1649        for k in range(quantity.centroid_values.shape[0]):
1650            assert allclose (quantity.centroid_values[k],
1651                             sum(quantity.vertex_values[k,:])/3)
1652
1653
1654
1655    def test_limit_edges_by_all_neighbours(self):
1656        quantity = Quantity(self.mesh4)
1657
1658        #Create a deliberate overshoot (e.g. from gradient computation)
1659        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1660
1661
1662        #Limit
1663        quantity.limit_edges_by_all_neighbours()
1664
1665        #Assert that central triangle is limited by neighbours
1666        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
1667        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
1668
1669        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1670        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
1671
1672        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
1673        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1674
1675
1676
1677        #Assert that quantities are conserved
1678        from Numeric import sum
1679        for k in range(quantity.centroid_values.shape[0]):
1680            assert allclose (quantity.centroid_values[k],
1681                             sum(quantity.vertex_values[k,:])/3)
1682
1683
1684    def test_limit_edges_by_neighbour(self):
1685        quantity = Quantity(self.mesh4)
1686
1687        #Create a deliberate overshoot (e.g. from gradient computation)
1688        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1689
1690
1691        #Limit
1692        quantity.limit_edges_by_neighbour()
1693
1694        #Assert that central triangle is limited by neighbours
1695        assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
1696        assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
1697
1698        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1699        assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
1700
1701        assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
1702        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1703
1704
1705
1706        #Assert that quantities are conserved
1707        from Numeric import sum
1708        for k in range(quantity.centroid_values.shape[0]):
1709            assert allclose (quantity.centroid_values[k],
1710                             sum(quantity.vertex_values[k,:])/3)
1711
1712    def test_limiter2(self):
1713        """Taken from test_shallow_water
1714        """
1715        quantity = Quantity(self.mesh4)
1716        quantity.domain.beta_w = 0.9
1717       
1718        #Test centroids
1719        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1720        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1721
1722
1723        #Extrapolate
1724        quantity.extrapolate_second_order()
1725
1726        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1727
1728        #Limit
1729        quantity.limit()
1730
1731        # limited value for beta_w = 0.9
1732       
1733        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1734        # limited values for beta_w = 0.5
1735        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1736
1737
1738        #Assert that quantities are conserved
1739        from Numeric import sum
1740        for k in range(quantity.centroid_values.shape[0]):
1741            assert allclose (quantity.centroid_values[k],
1742                             sum(quantity.vertex_values[k,:])/3)
1743
1744
1745
1746
1747
1748    def test_distribute_first_order(self):
1749        quantity = Quantity(self.mesh4)
1750
1751        #Test centroids
1752        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1753        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1754
1755
1756        #Extrapolate from centroid to vertices and edges
1757        quantity.extrapolate_first_order()
1758
1759        #Interpolate
1760        #quantity.interpolate_from_vertices_to_edges()
1761
1762        assert allclose(quantity.vertex_values,
1763                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1764        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1765                                               [3,3,3], [4, 4, 4]])
1766
1767
1768    def test_interpolate_from_vertices_to_edges(self):
1769        quantity = Quantity(self.mesh4)
1770
1771        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
1772
1773        quantity.interpolate_from_vertices_to_edges()
1774
1775        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
1776                                               [3., 2.5, 1.5],
1777                                               [3.5, 4.5, 3.],
1778                                               [2.5, 3.5, 2]])
1779
1780
1781    def test_interpolate_from_edges_to_vertices(self):
1782        quantity = Quantity(self.mesh4)
1783
1784        quantity.edge_values = array([[1., 1.5, 0.5],
1785                                [3., 2.5, 1.5],
1786                                [3.5, 4.5, 3.],
1787                                [2.5, 3.5, 2]],Float)
1788
1789        quantity.interpolate_from_edges_to_vertices()
1790
1791        assert allclose(quantity.vertex_values,
1792                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1793
1794
1795
1796    def test_distribute_second_order(self):
1797        quantity = Quantity(self.mesh4)
1798
1799        #Test centroids
1800        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1801        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1802
1803
1804        #Extrapolate
1805        quantity.extrapolate_second_order()
1806
1807        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1808
1809
1810    def test_update_explicit(self):
1811        quantity = Quantity(self.mesh4)
1812
1813        #Test centroids
1814        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1815        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1816
1817        #Set explicit_update
1818        quantity.explicit_update = array( [1.,1.,1.,1.] )
1819
1820        #Update with given timestep
1821        quantity.update(0.1)
1822
1823        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1824        assert allclose( quantity.centroid_values, x)
1825
1826    def test_update_semi_implicit(self):
1827        quantity = Quantity(self.mesh4)
1828
1829        #Test centroids
1830        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1831        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1832
1833        #Set semi implicit update
1834        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1835
1836        #Update with given timestep
1837        timestep = 0.1
1838        quantity.update(timestep)
1839
1840        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1841        denom = ones(4, Float)-timestep*sem
1842
1843        x = array([1, 2, 3, 4])/denom
1844        assert allclose( quantity.centroid_values, x)
1845
1846
1847    def test_both_updates(self):
1848        quantity = Quantity(self.mesh4)
1849
1850        #Test centroids
1851        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1852        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1853
1854        #Set explicit_update
1855        quantity.explicit_update = array( [4.,3.,2.,1.] )
1856
1857        #Set semi implicit update
1858        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1859
1860        #Update with given timestep
1861        timestep = 0.1
1862        quantity.update(0.1)
1863
1864        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1865        denom = ones(4, Float)-timestep*sem
1866
1867        x = array([1., 2., 3., 4.])
1868        x /= denom
1869        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1870
1871        assert allclose( quantity.centroid_values, x)
1872
1873
1874
1875
1876    #Test smoothing
1877    def test_smoothing(self):
1878
1879        from mesh_factory import rectangular
1880        from shallow_water import Domain, Transmissive_boundary
1881        from Numeric import zeros, Float
1882        from anuga.utilities.numerical_tools import mean
1883
1884        #Create basic mesh
1885        points, vertices, boundary = rectangular(2, 2)
1886
1887        #Create shallow water domain
1888        domain = Domain(points, vertices, boundary)
1889        domain.default_order=2
1890        domain.reduction = mean
1891
1892
1893        #Set some field values
1894        domain.set_quantity('elevation', lambda x,y: x)
1895        domain.set_quantity('friction', 0.03)
1896
1897
1898        ######################
1899        # Boundary conditions
1900        B = Transmissive_boundary(domain)
1901        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1902
1903
1904        ######################
1905        #Initial condition - with jumps
1906
1907        bed = domain.quantities['elevation'].vertex_values
1908        stage = zeros(bed.shape, Float)
1909
1910        h = 0.03
1911        for i in range(stage.shape[0]):
1912            if i % 2 == 0:
1913                stage[i,:] = bed[i,:] + h
1914            else:
1915                stage[i,:] = bed[i,:]
1916
1917        domain.set_quantity('stage', stage)
1918
1919        stage = domain.quantities['stage']
1920
1921        #Get smoothed stage
1922        A, V = stage.get_vertex_values(xy=False, smooth=True)
1923        Q = stage.vertex_values
1924
1925
1926        assert A.shape[0] == 9
1927        assert V.shape[0] == 8
1928        assert V.shape[1] == 3
1929
1930        #First four points
1931        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1932        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1933        assert allclose(A[2], Q[3,0])
1934        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1935
1936        #Center point
1937        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1938                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1939
1940
1941        #Check V
1942        assert allclose(V[0,:], [3,4,0])
1943        assert allclose(V[1,:], [1,0,4])
1944        assert allclose(V[2,:], [4,5,1])
1945        assert allclose(V[3,:], [2,1,5])
1946        assert allclose(V[4,:], [6,7,3])
1947        assert allclose(V[5,:], [4,3,7])
1948        assert allclose(V[6,:], [7,8,4])
1949        assert allclose(V[7,:], [5,4,8])
1950
1951        #Get smoothed stage with XY
1952        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1953
1954        assert allclose(A, A1)
1955        assert allclose(V, V1)
1956
1957        #Check XY
1958        assert allclose(X[4], 0.5)
1959        assert allclose(Y[4], 0.5)
1960
1961        assert allclose(X[7], 1.0)
1962        assert allclose(Y[7], 0.5)
1963
1964
1965
1966
1967    def test_vertex_values_no_smoothing(self):
1968
1969        from mesh_factory import rectangular
1970        from shallow_water import Domain, Transmissive_boundary
1971        from Numeric import zeros, Float
1972        from anuga.utilities.numerical_tools import mean
1973
1974
1975        #Create basic mesh
1976        points, vertices, boundary = rectangular(2, 2)
1977
1978        #Create shallow water domain
1979        domain = Domain(points, vertices, boundary)
1980        domain.default_order=2
1981        domain.reduction = mean
1982
1983
1984        #Set some field values
1985        domain.set_quantity('elevation', lambda x,y: x)
1986        domain.set_quantity('friction', 0.03)
1987
1988
1989        ######################
1990        #Initial condition - with jumps
1991
1992        bed = domain.quantities['elevation'].vertex_values
1993        stage = zeros(bed.shape, Float)
1994
1995        h = 0.03
1996        for i in range(stage.shape[0]):
1997            if i % 2 == 0:
1998                stage[i,:] = bed[i,:] + h
1999            else:
2000                stage[i,:] = bed[i,:]
2001
2002        domain.set_quantity('stage', stage)
2003
2004        #Get stage
2005        stage = domain.quantities['stage']
2006        A, V = stage.get_vertex_values(xy=False, smooth=False)
2007        Q = stage.vertex_values.flat
2008
2009        for k in range(8):
2010            assert allclose(A[k], Q[k])
2011
2012
2013        for k in range(8):
2014            assert V[k, 0] == 3*k
2015            assert V[k, 1] == 3*k+1
2016            assert V[k, 2] == 3*k+2
2017
2018
2019
2020        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
2021
2022
2023        assert allclose(A, A1)
2024        assert allclose(V, V1)
2025
2026        #Check XY
2027        assert allclose(X[1], 0.5)
2028        assert allclose(Y[1], 0.5)
2029        assert allclose(X[4], 0.0)
2030        assert allclose(Y[4], 0.0)
2031        assert allclose(X[12], 1.0)
2032        assert allclose(Y[12], 0.0)
2033
2034
2035
2036    def set_array_values_by_index(self):
2037
2038        from mesh_factory import rectangular
2039        from shallow_water import Domain
2040        from Numeric import zeros, Float
2041
2042        #Create basic mesh
2043        points, vertices, boundary = rectangular(1, 1)
2044
2045        #Create shallow water domain
2046        domain = Domain(points, vertices, boundary)
2047        #print "domain.number_of_elements ",domain.number_of_elements
2048        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
2049        value = [7]
2050        indices = [1]
2051        quantity.set_array_values_by_index(value,
2052                                           location = 'centroids',
2053                                           indices = indices)
2054        #print "quantity.centroid_values",quantity.centroid_values
2055
2056        assert allclose(quantity.centroid_values, [1,7])
2057
2058        quantity.set_array_values([15,20,25], indices = indices)
2059        assert allclose(quantity.centroid_values, [1,20])
2060
2061        quantity.set_array_values([15,20,25], indices = indices)
2062        assert allclose(quantity.centroid_values, [1,20])
2063
2064    def test_setting_some_vertex_values(self):
2065        """
2066        set values based on triangle lists.
2067        """
2068        from mesh_factory import rectangular
2069        from shallow_water import Domain
2070        from Numeric import zeros, Float
2071
2072        #Create basic mesh
2073        points, vertices, boundary = rectangular(1, 3)
2074        #print "vertices",vertices
2075        #Create shallow water domain
2076        domain = Domain(points, vertices, boundary)
2077        #print "domain.number_of_elements ",domain.number_of_elements
2078        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2079                                    [4,4,4],[5,5,5],[6,6,6]])
2080
2081
2082        # Check that constants work
2083        value = 7
2084        indices = [1]
2085        quantity.set_values(value,
2086                            location = 'centroids',
2087                            indices = indices)
2088        #print "quantity.centroid_values",quantity.centroid_values
2089        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
2090       
2091        value = [7]
2092        indices = [1]
2093        quantity.set_values(value,
2094                            location = 'centroids',
2095                            indices = indices)
2096        #print "quantity.centroid_values",quantity.centroid_values
2097        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
2098
2099        value = [[15,20,25]]
2100        quantity.set_values(value, indices = indices)
2101        #print "1 quantity.vertex_values",quantity.vertex_values
2102        assert allclose(quantity.vertex_values[1], value[0])
2103
2104
2105        #print "quantity",quantity.vertex_values
2106        values = [10,100,50]
2107        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
2108        #print "2 quantity.vertex_values",quantity.vertex_values
2109        assert allclose(quantity.vertex_values[0], [10,10,10])
2110        assert allclose(quantity.vertex_values[5], [50,50,50])
2111        #quantity.interpolate()
2112        #print "quantity.centroid_values",quantity.centroid_values
2113        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
2114
2115
2116        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2117                                    [4,4,4],[5,5,5],[6,6,6]])
2118        values = [10,100,50]
2119        #this will be per unique vertex, indexing the vertices
2120        #print "quantity.vertex_values",quantity.vertex_values
2121        quantity.set_values(values, indices = [0,1,5])
2122        #print "quantity.vertex_values",quantity.vertex_values
2123        assert allclose(quantity.vertex_values[0], [1,50,10])
2124        assert allclose(quantity.vertex_values[5], [6,6,6])
2125        assert allclose(quantity.vertex_values[1], [100,10,50])
2126
2127        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2128                                    [4,4,4],[5,5,5],[6,6,6]])
2129        values = [[31,30,29],[400,400,400],[1000,999,998]]
2130        quantity.set_values(values, indices = [3,3,5])
2131        quantity.interpolate()
2132        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
2133
2134        values = [[1,1,1],[2,2,2],[3,3,3],
2135                                    [4,4,4],[5,5,5],[6,6,6]]
2136        quantity.set_values(values)
2137
2138        # testing the standard set values by vertex
2139        # indexed by vertex_id in general_mesh.coordinates
2140        values = [0,1,2,3,4,5,6,7]
2141
2142        quantity.set_values(values)
2143        #print "1 quantity.vertex_values",quantity.vertex_values
2144        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
2145                                                [ 1.,  0.,  5.],
2146                                                [ 5.,  6.,  1.],
2147                                                [ 2.,  1.,  6.],
2148                                                [ 6.,  7.,  2.],
2149                                                [ 3.,  2.,  7.]])
2150
2151    def test_setting_unique_vertex_values(self):
2152        """
2153        set values based on unique_vertex lists.
2154        """
2155        from mesh_factory import rectangular
2156        from shallow_water import Domain
2157        from Numeric import zeros, Float
2158
2159        #Create basic mesh
2160        points, vertices, boundary = rectangular(1, 3)
2161        #print "vertices",vertices
2162        #Create shallow water domain
2163        domain = Domain(points, vertices, boundary)
2164        #print "domain.number_of_elements ",domain.number_of_elements
2165        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2166                                    [4,4,4],[5,5,5]])
2167        value = 7
2168        indices = [1,5]
2169        quantity.set_values(value,
2170                            location = 'unique vertices',
2171                            indices = indices)
2172        #print "quantity.centroid_values",quantity.centroid_values
2173        assert allclose(quantity.vertex_values[0], [0,7,0])
2174        assert allclose(quantity.vertex_values[1], [7,1,7])
2175        assert allclose(quantity.vertex_values[2], [7,2,7])
2176
2177
2178    def test_get_values(self):
2179        """
2180        get values based on triangle lists.
2181        """
2182        from mesh_factory import rectangular
2183        from shallow_water import Domain
2184        from Numeric import zeros, Float
2185
2186        #Create basic mesh
2187        points, vertices, boundary = rectangular(1, 3)
2188
2189        #print "points",points
2190        #print "vertices",vertices
2191        #print "boundary",boundary
2192
2193        #Create shallow water domain
2194        domain = Domain(points, vertices, boundary)
2195        #print "domain.number_of_elements ",domain.number_of_elements
2196        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2197                                    [4,4,4],[5,5,5]])
2198
2199        #print "quantity.get_values(location = 'unique vertices')", \
2200        #      quantity.get_values(location = 'unique vertices')
2201
2202        #print "quantity.get_values(location = 'unique vertices')", \
2203        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2204        #                          location = 'unique vertices')
2205
2206        answer = [0.5,2,4,5,0,1,3,4.5]
2207        assert allclose(answer,
2208                        quantity.get_values(location = 'unique vertices'))
2209
2210        indices = [0,5,3]
2211        answer = [0.5,1,5]
2212        assert allclose(answer,
2213                        quantity.get_values(indices=indices, \
2214                                            location = 'unique vertices'))
2215        #print "quantity.centroid_values",quantity.centroid_values
2216        #print "quantity.get_values(location = 'centroids') ",\
2217        #      quantity.get_values(location = 'centroids')
2218
2219
2220
2221
2222    def test_get_values_2(self):
2223        """Different mesh (working with domain object) - also check centroids.
2224        """
2225
2226       
2227        a = [0.0, 0.0]
2228        b = [0.0, 2.0]
2229        c = [2.0,0.0]
2230        d = [0.0, 4.0]
2231        e = [2.0, 2.0]
2232        f = [4.0,0.0]
2233
2234        points = [a, b, c, d, e, f]
2235        #bac, bce, ecf, dbe
2236        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2237
2238        domain = Domain(points, vertices)
2239
2240        quantity = Quantity(domain)
2241        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2242       
2243        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
2244        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
2245
2246
2247        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
2248                                                                   [4,2,6],
2249                                                                   [6,2,4],
2250                                                                   [8,4,6]])
2251       
2252        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
2253                                                                                  [8,4,6]])
2254
2255
2256        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
2257                                                                [4,5,3],
2258                                                                [3,5,4],
2259                                                                [5,7,6]])
2260        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
2261                        [[4,5,3],
2262                         [5,7,6]])       
2263
2264        # Check averaging over vertices
2265        #a: 0
2266        #b: (4+4+4)/3
2267        #c: (2+2+2)/3
2268        #d: 8
2269        #e: (6+6+6)/3       
2270        #f: 4
2271        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
2272                                                                                 
2273       
2274
2275
2276
2277
2278    def test_get_interpolated_values(self):
2279
2280        from mesh_factory import rectangular
2281        from shallow_water import Domain
2282        from Numeric import zeros, Float
2283
2284        #Create basic mesh
2285        points, vertices, boundary = rectangular(1, 3)
2286        domain = Domain(points, vertices, boundary)
2287
2288        #Constant values
2289        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2290                                    [4,4,4],[5,5,5]])
2291
2292       
2293
2294        # Get interpolated values at centroids
2295        interpolation_points = domain.get_centroid_coordinates()
2296        answer = quantity.get_values(location='centroids')
2297
2298       
2299        #print quantity.get_values(points=interpolation_points)
2300        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
2301
2302
2303        #Arbitrary values
2304        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
2305                                    [1,4,-9],[2,5,0]])
2306
2307
2308        # Get interpolated values at centroids
2309        interpolation_points = domain.get_centroid_coordinates()
2310        answer = quantity.get_values(location='centroids')
2311        #print answer
2312        #print quantity.get_values(points=interpolation_points)
2313        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
2314                       
2315
2316        #FIXME TODO
2317        #indices = [0,5,3]
2318        #answer = [0.5,1,5]
2319        #assert allclose(answer,
2320        #                quantity.get_values(indices=indices, \
2321        #                                    location = 'unique vertices'))
2322
2323
2324
2325
2326    def test_get_interpolated_values_2(self):
2327        a = [0.0, 0.0]
2328        b = [0.0, 2.0]
2329        c = [2.0,0.0]
2330        d = [0.0, 4.0]
2331        e = [2.0, 2.0]
2332        f = [4.0,0.0]
2333
2334        points = [a, b, c, d, e, f]
2335        #bac, bce, ecf, dbe
2336        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2337
2338        domain = Domain(points, vertices)
2339
2340        quantity = Quantity(domain)
2341        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2342
2343        #First pick one point
2344        x, y = 2.0/3, 8.0/3
2345        v = quantity.get_values(interpolation_points = [[x,y]])
2346        assert allclose(v, 6)       
2347
2348        # Then another to test that algorithm won't blindly
2349        # reuse interpolation matrix
2350        x, y = 4.0/3, 4.0/3
2351        v = quantity.get_values(interpolation_points = [[x,y]])
2352        assert allclose(v, 4)       
2353
2354
2355
2356
2357    def test_getting_some_vertex_values(self):
2358        """
2359        get values based on triangle lists.
2360        """
2361        from mesh_factory import rectangular
2362        from shallow_water import Domain
2363        from Numeric import zeros, Float
2364
2365        #Create basic mesh
2366        points, vertices, boundary = rectangular(1, 3)
2367
2368        #print "points",points
2369        #print "vertices",vertices
2370        #print "boundary",boundary
2371
2372        #Create shallow water domain
2373        domain = Domain(points, vertices, boundary)
2374        #print "domain.number_of_elements ",domain.number_of_elements
2375        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2376                                    [4,4,4],[5,5,5],[6,6,6]])
2377        value = [7]
2378        indices = [1]
2379        quantity.set_values(value,
2380                            location = 'centroids',
2381                            indices = indices)
2382        #print "quantity.centroid_values",quantity.centroid_values
2383        #print "quantity.get_values(location = 'centroids') ",\
2384        #      quantity.get_values(location = 'centroids')
2385        assert allclose(quantity.centroid_values,
2386                        quantity.get_values(location = 'centroids'))
2387
2388
2389        value = [[15,20,25]]
2390        quantity.set_values(value, indices = indices)
2391        #print "1 quantity.vertex_values",quantity.vertex_values
2392        assert allclose(quantity.vertex_values, quantity.get_values())
2393
2394        assert allclose(quantity.edge_values,
2395                        quantity.get_values(location = 'edges'))
2396
2397        # get a subset of elements
2398        subset = quantity.get_values(location='centroids', indices=[0,5])
2399        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
2400        assert allclose(subset, answer)
2401
2402
2403        subset = quantity.get_values(location='edges', indices=[0,5])
2404        answer = [quantity.edge_values[0],quantity.edge_values[5]]
2405        #print "subset",subset
2406        #print "answer",answer
2407        assert allclose(subset, answer)
2408
2409        subset = quantity.get_values( indices=[1,5])
2410        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
2411        #print "subset",subset
2412        #print "answer",answer
2413        assert allclose(subset, answer)
2414
2415    def test_smooth_vertex_values(self):
2416        """
2417        get values based on triangle lists.
2418        """
2419        from mesh_factory import rectangular
2420        from shallow_water import Domain
2421        from Numeric import zeros, Float
2422
2423        #Create basic mesh
2424        points, vertices, boundary = rectangular(2, 2)
2425
2426        #print "points",points
2427        #print "vertices",vertices
2428        #print "boundary",boundary
2429
2430        #Create shallow water domain
2431        domain = Domain(points, vertices, boundary)
2432        #print "domain.number_of_elements ",domain.number_of_elements
2433        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2434                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
2435
2436        #print "quantity.get_values(location = 'unique vertices')", \
2437        #      quantity.get_values(location = 'unique vertices')
2438
2439        #print "quantity.get_values(location = 'unique vertices')", \
2440        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2441        #                          location = 'unique vertices')
2442
2443        #print quantity.get_values(location = 'unique vertices')
2444        #print quantity.domain.number_of_triangles_per_node
2445        #print quantity.vertex_values
2446       
2447        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
2448        #assert allclose(answer,
2449        #                quantity.get_values(location = 'unique vertices'))
2450
2451        quantity.smooth_vertex_values()
2452
2453        #print quantity.vertex_values
2454
2455
2456        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
2457                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
2458       
2459        assert allclose(answer_vertex_values,
2460                        quantity.vertex_values)
2461        #print "quantity.centroid_values",quantity.centroid_values
2462        #print "quantity.get_values(location = 'centroids') ",\
2463        #      quantity.get_values(location = 'centroids')
2464
2465
2466
2467#-------------------------------------------------------------
2468if __name__ == "__main__":
2469    suite = unittest.makeSuite(Test_Quantity, 'test')   
2470    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
2471
2472    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
2473    #print "restricted test"
2474    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
2475    runner = unittest.TextTestRunner()
2476    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.