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

Last change on this file since 5897 was 5897, checked in by ole, 16 years ago

Reverted numpy changes to the trunk that should have been made to the branch.
The command was svn merge -r 5895:5890 .
File size: 84.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(interpolation_points=interpolation_points)
2313        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points,
2314                                                    verbose=False))       
2315                       
2316
2317        #FIXME TODO
2318        #indices = [0,5,3]
2319        #answer = [0.5,1,5]
2320        #assert allclose(answer,
2321        #                quantity.get_values(indices=indices, \
2322        #                                    location = 'unique vertices'))
2323
2324
2325
2326
2327    def test_get_interpolated_values_2(self):
2328        a = [0.0, 0.0]
2329        b = [0.0, 2.0]
2330        c = [2.0,0.0]
2331        d = [0.0, 4.0]
2332        e = [2.0, 2.0]
2333        f = [4.0,0.0]
2334
2335        points = [a, b, c, d, e, f]
2336        #bac, bce, ecf, dbe
2337        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2338
2339        domain = Domain(points, vertices)
2340
2341        quantity = Quantity(domain)
2342        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2343
2344        #First pick one point
2345        x, y = 2.0/3, 8.0/3
2346        v = quantity.get_values(interpolation_points = [[x,y]])
2347        assert allclose(v, 6)       
2348
2349        # Then another to test that algorithm won't blindly
2350        # reuse interpolation matrix
2351        x, y = 4.0/3, 4.0/3
2352        v = quantity.get_values(interpolation_points = [[x,y]])
2353        assert allclose(v, 4)       
2354
2355
2356
2357    def test_get_interpolated_values_with_georef(self):
2358   
2359        zone = 56
2360        xllcorner = 308500
2361        yllcorner = 6189000
2362        a = [0.0, 0.0]
2363        b = [0.0, 2.0]
2364        c = [2.0,0.0]
2365        d = [0.0, 4.0]
2366        e = [2.0, 2.0]
2367        f = [4.0,0.0]
2368
2369        points = [a, b, c, d, e, f]
2370        #bac, bce, ecf, dbe
2371        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2372
2373        domain = Domain(points, vertices,
2374                        geo_reference=Geo_reference(zone,xllcorner,yllcorner))
2375
2376        quantity = Quantity(domain)
2377        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2378
2379        #First pick one point (and turn it into absolute coordinates)
2380        x, y = 2.0/3, 8.0/3
2381        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
2382        assert allclose(v, 6)
2383       
2384
2385        # Then another to test that algorithm won't blindly
2386        # reuse interpolation matrix
2387        x, y = 4.0/3, 4.0/3
2388        v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]])
2389        assert allclose(v, 4)       
2390       
2391        # Try two points
2392        pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 
2393               [4.0/3 + xllcorner, 4.0/3 + yllcorner]]         
2394        v = quantity.get_values(interpolation_points=pts)
2395        assert allclose(v, [6, 4])               
2396       
2397        # Test it using the geospatial data format with absolute input points and default georef
2398        pts = Geospatial_data(data_points=pts)
2399        v = quantity.get_values(interpolation_points=pts)
2400        assert allclose(v, [6, 4])                               
2401       
2402       
2403        # Test it using the geospatial data format with relative input points
2404        pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 
2405                              geo_reference=Geo_reference(zone,xllcorner,yllcorner))
2406        v = quantity.get_values(interpolation_points=pts)
2407        assert allclose(v, [6, 4])                       
2408       
2409       
2410       
2411
2412    def test_getting_some_vertex_values(self):
2413        """
2414        get values based on triangle lists.
2415        """
2416        from mesh_factory import rectangular
2417        from shallow_water import Domain
2418        from Numeric import zeros, Float
2419
2420        #Create basic mesh
2421        points, vertices, boundary = rectangular(1, 3)
2422
2423        #print "points",points
2424        #print "vertices",vertices
2425        #print "boundary",boundary
2426
2427        #Create shallow water domain
2428        domain = Domain(points, vertices, boundary)
2429        #print "domain.number_of_elements ",domain.number_of_elements
2430        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2431                                    [4,4,4],[5,5,5],[6,6,6]])
2432        value = [7]
2433        indices = [1]
2434        quantity.set_values(value,
2435                            location = 'centroids',
2436                            indices = indices)
2437        #print "quantity.centroid_values",quantity.centroid_values
2438        #print "quantity.get_values(location = 'centroids') ",\
2439        #      quantity.get_values(location = 'centroids')
2440        assert allclose(quantity.centroid_values,
2441                        quantity.get_values(location = 'centroids'))
2442
2443
2444        value = [[15,20,25]]
2445        quantity.set_values(value, indices = indices)
2446        #print "1 quantity.vertex_values",quantity.vertex_values
2447        assert allclose(quantity.vertex_values, quantity.get_values())
2448
2449        assert allclose(quantity.edge_values,
2450                        quantity.get_values(location = 'edges'))
2451
2452        # get a subset of elements
2453        subset = quantity.get_values(location='centroids', indices=[0,5])
2454        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
2455        assert allclose(subset, answer)
2456
2457
2458        subset = quantity.get_values(location='edges', indices=[0,5])
2459        answer = [quantity.edge_values[0],quantity.edge_values[5]]
2460        #print "subset",subset
2461        #print "answer",answer
2462        assert allclose(subset, answer)
2463
2464        subset = quantity.get_values( indices=[1,5])
2465        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
2466        #print "subset",subset
2467        #print "answer",answer
2468        assert allclose(subset, answer)
2469
2470    def test_smooth_vertex_values(self):
2471        """
2472        get values based on triangle lists.
2473        """
2474        from mesh_factory import rectangular
2475        from shallow_water import Domain
2476        from Numeric import zeros, Float
2477
2478        #Create basic mesh
2479        points, vertices, boundary = rectangular(2, 2)
2480
2481        #print "points",points
2482        #print "vertices",vertices
2483        #print "boundary",boundary
2484
2485        #Create shallow water domain
2486        domain = Domain(points, vertices, boundary)
2487        #print "domain.number_of_elements ",domain.number_of_elements
2488        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2489                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
2490
2491        #print "quantity.get_values(location = 'unique vertices')", \
2492        #      quantity.get_values(location = 'unique vertices')
2493
2494        #print "quantity.get_values(location = 'unique vertices')", \
2495        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2496        #                          location = 'unique vertices')
2497
2498        #print quantity.get_values(location = 'unique vertices')
2499        #print quantity.domain.number_of_triangles_per_node
2500        #print quantity.vertex_values
2501       
2502        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
2503        #assert allclose(answer,
2504        #                quantity.get_values(location = 'unique vertices'))
2505
2506        quantity.smooth_vertex_values()
2507
2508        #print quantity.vertex_values
2509
2510
2511        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
2512                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
2513       
2514        assert allclose(answer_vertex_values,
2515                        quantity.vertex_values)
2516        #print "quantity.centroid_values",quantity.centroid_values
2517        #print "quantity.get_values(location = 'centroids') ",\
2518        #      quantity.get_values(location = 'centroids')
2519
2520
2521
2522#-------------------------------------------------------------
2523if __name__ == "__main__":
2524    suite = unittest.makeSuite(Test_Quantity, 'test')   
2525    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
2526
2527    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
2528    #print "restricted test"
2529    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
2530    runner = unittest.TextTestRunner()
2531    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.