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

Last change on this file since 4883 was 4883, checked in by steve, 16 years ago

Starting to develop edge limiting and vertex limiting in
quantity (and later in shallow_water_domain

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