source: branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py @ 6883

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

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