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

Last change on this file since 6145 was 6145, checked in by rwilson, 16 years ago

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