source: anuga_work/development/anuga_1d/test_quantity.py @ 5858

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