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

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

Added another polygon test for set_values
(both polygon tests are currently disabled)

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