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

Last change on this file since 6304 was 6304, checked in by rwilson, 15 years ago

Initial commit of numpy changes. Still a long way to go.

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