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

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

Working on limiters
File size: 76.6 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_limit_vertices_by_all_neighbours(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_vertices_by_all_neighbours()
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
1477    def test_limit_edges_by_all_neighbours(self):
1478        quantity = Conserved_quantity(self.mesh4)
1479
1480        #Create a deliberate overshoot (e.g. from gradient computation)
1481        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1482
1483
1484        #Limit
1485        quantity.limit_edges_by_all_neighbours()
1486
1487        #Assert that central triangle is limited by neighbours
1488        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
1489        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
1490
1491        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1492        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
1493
1494        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
1495        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1496
1497
1498
1499        #Assert that quantities are conserved
1500        from Numeric import sum
1501        for k in range(quantity.centroid_values.shape[0]):
1502            assert allclose (quantity.centroid_values[k],
1503                             sum(quantity.vertex_values[k,:])/3)
1504
1505
1506    def test_limit_edges_by_neighbour(self):
1507        quantity = Conserved_quantity(self.mesh4)
1508
1509        #Create a deliberate overshoot (e.g. from gradient computation)
1510        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
1511
1512
1513        #Limit
1514        quantity.limit_edges_by_neighbour()
1515
1516        #Assert that central triangle is limited by neighbours
1517        assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
1518        assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
1519
1520        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
1521        assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
1522
1523        assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
1524        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
1525
1526
1527
1528        #Assert that quantities are conserved
1529        from Numeric import sum
1530        for k in range(quantity.centroid_values.shape[0]):
1531            assert allclose (quantity.centroid_values[k],
1532                             sum(quantity.vertex_values[k,:])/3)
1533
1534    def test_limiter2(self):
1535        """Taken from test_shallow_water
1536        """
1537        quantity = Conserved_quantity(self.mesh4)
1538        quantity.domain.beta_w = 0.9
1539       
1540        #Test centroids
1541        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1542        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1543
1544
1545        #Extrapolate
1546        quantity.extrapolate_second_order()
1547
1548        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1549
1550        #Limit
1551        quantity.limit()
1552
1553        # limited value for beta_w = 0.9
1554       
1555        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1556        # limited values for beta_w = 0.5
1557        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1558
1559
1560        #Assert that quantities are conserved
1561        from Numeric import sum
1562        for k in range(quantity.centroid_values.shape[0]):
1563            assert allclose (quantity.centroid_values[k],
1564                             sum(quantity.vertex_values[k,:])/3)
1565
1566
1567
1568
1569
1570    def test_distribute_first_order(self):
1571        quantity = Conserved_quantity(self.mesh4)
1572
1573        #Test centroids
1574        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1575        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1576
1577
1578        #Extrapolate from centroid to vertices and edges
1579        quantity.extrapolate_first_order()
1580
1581        #Interpolate
1582        #quantity.interpolate_from_vertices_to_edges()
1583
1584        assert allclose(quantity.vertex_values,
1585                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1586        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1587                                               [3,3,3], [4, 4, 4]])
1588
1589
1590    def test_interpolate_from_vertices_to_edges(self):
1591        quantity = Quantity(self.mesh4)
1592
1593        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
1594
1595        quantity.interpolate_from_vertices_to_edges()
1596
1597        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
1598                                               [3., 2.5, 1.5],
1599                                               [3.5, 4.5, 3.],
1600                                               [2.5, 3.5, 2]])
1601
1602
1603    def test_interpolate_from_edges_to_vertices(self):
1604        quantity = Quantity(self.mesh4)
1605
1606        quantity.edge_values = array([[1., 1.5, 0.5],
1607                                [3., 2.5, 1.5],
1608                                [3.5, 4.5, 3.],
1609                                [2.5, 3.5, 2]],Float)
1610
1611        quantity.interpolate_from_edges_to_vertices()
1612
1613        assert allclose(quantity.vertex_values,
1614                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
1615
1616
1617
1618    def test_distribute_second_order(self):
1619        quantity = Conserved_quantity(self.mesh4)
1620
1621        #Test centroids
1622        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1623        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1624
1625
1626        #Extrapolate
1627        quantity.extrapolate_second_order()
1628
1629        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1630
1631
1632    def test_update_explicit(self):
1633        quantity = Conserved_quantity(self.mesh4)
1634
1635        #Test centroids
1636        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1637        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1638
1639        #Set explicit_update
1640        quantity.explicit_update = array( [1.,1.,1.,1.] )
1641
1642        #Update with given timestep
1643        quantity.update(0.1)
1644
1645        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1646        assert allclose( quantity.centroid_values, x)
1647
1648    def test_update_semi_implicit(self):
1649        quantity = Conserved_quantity(self.mesh4)
1650
1651        #Test centroids
1652        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1653        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1654
1655        #Set semi implicit update
1656        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1657
1658        #Update with given timestep
1659        timestep = 0.1
1660        quantity.update(timestep)
1661
1662        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1663        denom = ones(4, Float)-timestep*sem
1664
1665        x = array([1, 2, 3, 4])/denom
1666        assert allclose( quantity.centroid_values, x)
1667
1668
1669    def test_both_updates(self):
1670        quantity = Conserved_quantity(self.mesh4)
1671
1672        #Test centroids
1673        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1674        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1675
1676        #Set explicit_update
1677        quantity.explicit_update = array( [4.,3.,2.,1.] )
1678
1679        #Set semi implicit update
1680        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1681
1682        #Update with given timestep
1683        timestep = 0.1
1684        quantity.update(0.1)
1685
1686        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1687        denom = ones(4, Float)-timestep*sem
1688
1689        x = array([1., 2., 3., 4.])
1690        x /= denom
1691        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1692
1693        assert allclose( quantity.centroid_values, x)
1694
1695
1696
1697
1698    #Test smoothing
1699    def test_smoothing(self):
1700
1701        from mesh_factory import rectangular
1702        from shallow_water import Domain, Transmissive_boundary
1703        from Numeric import zeros, Float
1704        from anuga.utilities.numerical_tools import mean
1705
1706        #Create basic mesh
1707        points, vertices, boundary = rectangular(2, 2)
1708
1709        #Create shallow water domain
1710        domain = Domain(points, vertices, boundary)
1711        domain.default_order=2
1712        domain.reduction = mean
1713
1714
1715        #Set some field values
1716        domain.set_quantity('elevation', lambda x,y: x)
1717        domain.set_quantity('friction', 0.03)
1718
1719
1720        ######################
1721        # Boundary conditions
1722        B = Transmissive_boundary(domain)
1723        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1724
1725
1726        ######################
1727        #Initial condition - with jumps
1728
1729        bed = domain.quantities['elevation'].vertex_values
1730        stage = zeros(bed.shape, Float)
1731
1732        h = 0.03
1733        for i in range(stage.shape[0]):
1734            if i % 2 == 0:
1735                stage[i,:] = bed[i,:] + h
1736            else:
1737                stage[i,:] = bed[i,:]
1738
1739        domain.set_quantity('stage', stage)
1740
1741        stage = domain.quantities['stage']
1742
1743        #Get smoothed stage
1744        A, V = stage.get_vertex_values(xy=False, smooth=True)
1745        Q = stage.vertex_values
1746
1747
1748        assert A.shape[0] == 9
1749        assert V.shape[0] == 8
1750        assert V.shape[1] == 3
1751
1752        #First four points
1753        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1754        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1755        assert allclose(A[2], Q[3,0])
1756        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1757
1758        #Center point
1759        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1760                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1761
1762
1763        #Check V
1764        assert allclose(V[0,:], [3,4,0])
1765        assert allclose(V[1,:], [1,0,4])
1766        assert allclose(V[2,:], [4,5,1])
1767        assert allclose(V[3,:], [2,1,5])
1768        assert allclose(V[4,:], [6,7,3])
1769        assert allclose(V[5,:], [4,3,7])
1770        assert allclose(V[6,:], [7,8,4])
1771        assert allclose(V[7,:], [5,4,8])
1772
1773        #Get smoothed stage with XY
1774        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1775
1776        assert allclose(A, A1)
1777        assert allclose(V, V1)
1778
1779        #Check XY
1780        assert allclose(X[4], 0.5)
1781        assert allclose(Y[4], 0.5)
1782
1783        assert allclose(X[7], 1.0)
1784        assert allclose(Y[7], 0.5)
1785
1786
1787
1788
1789    def test_vertex_values_no_smoothing(self):
1790
1791        from mesh_factory import rectangular
1792        from shallow_water import Domain, Transmissive_boundary
1793        from Numeric import zeros, Float
1794        from anuga.utilities.numerical_tools import mean
1795
1796
1797        #Create basic mesh
1798        points, vertices, boundary = rectangular(2, 2)
1799
1800        #Create shallow water domain
1801        domain = Domain(points, vertices, boundary)
1802        domain.default_order=2
1803        domain.reduction = mean
1804
1805
1806        #Set some field values
1807        domain.set_quantity('elevation', lambda x,y: x)
1808        domain.set_quantity('friction', 0.03)
1809
1810
1811        ######################
1812        #Initial condition - with jumps
1813
1814        bed = domain.quantities['elevation'].vertex_values
1815        stage = zeros(bed.shape, Float)
1816
1817        h = 0.03
1818        for i in range(stage.shape[0]):
1819            if i % 2 == 0:
1820                stage[i,:] = bed[i,:] + h
1821            else:
1822                stage[i,:] = bed[i,:]
1823
1824        domain.set_quantity('stage', stage)
1825
1826        #Get stage
1827        stage = domain.quantities['stage']
1828        A, V = stage.get_vertex_values(xy=False, smooth=False)
1829        Q = stage.vertex_values.flat
1830
1831        for k in range(8):
1832            assert allclose(A[k], Q[k])
1833
1834
1835        for k in range(8):
1836            assert V[k, 0] == 3*k
1837            assert V[k, 1] == 3*k+1
1838            assert V[k, 2] == 3*k+2
1839
1840
1841
1842        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1843
1844
1845        assert allclose(A, A1)
1846        assert allclose(V, V1)
1847
1848        #Check XY
1849        assert allclose(X[1], 0.5)
1850        assert allclose(Y[1], 0.5)
1851        assert allclose(X[4], 0.0)
1852        assert allclose(Y[4], 0.0)
1853        assert allclose(X[12], 1.0)
1854        assert allclose(Y[12], 0.0)
1855
1856
1857
1858    def set_array_values_by_index(self):
1859
1860        from mesh_factory import rectangular
1861        from shallow_water import Domain
1862        from Numeric import zeros, Float
1863
1864        #Create basic mesh
1865        points, vertices, boundary = rectangular(1, 1)
1866
1867        #Create shallow water domain
1868        domain = Domain(points, vertices, boundary)
1869        #print "domain.number_of_elements ",domain.number_of_elements
1870        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1871        value = [7]
1872        indices = [1]
1873        quantity.set_array_values_by_index(value,
1874                                           location = 'centroids',
1875                                           indices = indices)
1876        #print "quantity.centroid_values",quantity.centroid_values
1877
1878        assert allclose(quantity.centroid_values, [1,7])
1879
1880        quantity.set_array_values([15,20,25], indices = indices)
1881        assert allclose(quantity.centroid_values, [1,20])
1882
1883        quantity.set_array_values([15,20,25], indices = indices)
1884        assert allclose(quantity.centroid_values, [1,20])
1885
1886    def test_setting_some_vertex_values(self):
1887        """
1888        set values based on triangle lists.
1889        """
1890        from mesh_factory import rectangular
1891        from shallow_water import Domain
1892        from Numeric import zeros, Float
1893
1894        #Create basic mesh
1895        points, vertices, boundary = rectangular(1, 3)
1896        #print "vertices",vertices
1897        #Create shallow water domain
1898        domain = Domain(points, vertices, boundary)
1899        #print "domain.number_of_elements ",domain.number_of_elements
1900        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1901                                    [4,4,4],[5,5,5],[6,6,6]])
1902
1903
1904        # Check that constants work
1905        value = 7
1906        indices = [1]
1907        quantity.set_values(value,
1908                            location = 'centroids',
1909                            indices = indices)
1910        #print "quantity.centroid_values",quantity.centroid_values
1911        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1912       
1913        value = [7]
1914        indices = [1]
1915        quantity.set_values(value,
1916                            location = 'centroids',
1917                            indices = indices)
1918        #print "quantity.centroid_values",quantity.centroid_values
1919        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1920
1921        value = [[15,20,25]]
1922        quantity.set_values(value, indices = indices)
1923        #print "1 quantity.vertex_values",quantity.vertex_values
1924        assert allclose(quantity.vertex_values[1], value[0])
1925
1926
1927        #print "quantity",quantity.vertex_values
1928        values = [10,100,50]
1929        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1930        #print "2 quantity.vertex_values",quantity.vertex_values
1931        assert allclose(quantity.vertex_values[0], [10,10,10])
1932        assert allclose(quantity.vertex_values[5], [50,50,50])
1933        #quantity.interpolate()
1934        #print "quantity.centroid_values",quantity.centroid_values
1935        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1936
1937
1938        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1939                                    [4,4,4],[5,5,5],[6,6,6]])
1940        values = [10,100,50]
1941        #this will be per unique vertex, indexing the vertices
1942        #print "quantity.vertex_values",quantity.vertex_values
1943        quantity.set_values(values, indices = [0,1,5])
1944        #print "quantity.vertex_values",quantity.vertex_values
1945        assert allclose(quantity.vertex_values[0], [1,50,10])
1946        assert allclose(quantity.vertex_values[5], [6,6,6])
1947        assert allclose(quantity.vertex_values[1], [100,10,50])
1948
1949        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1950                                    [4,4,4],[5,5,5],[6,6,6]])
1951        values = [[31,30,29],[400,400,400],[1000,999,998]]
1952        quantity.set_values(values, indices = [3,3,5])
1953        quantity.interpolate()
1954        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1955
1956        values = [[1,1,1],[2,2,2],[3,3,3],
1957                                    [4,4,4],[5,5,5],[6,6,6]]
1958        quantity.set_values(values)
1959
1960        # testing the standard set values by vertex
1961        # indexed by vertex_id in general_mesh.coordinates
1962        values = [0,1,2,3,4,5,6,7]
1963
1964        quantity.set_values(values)
1965        #print "1 quantity.vertex_values",quantity.vertex_values
1966        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1967                                                [ 1.,  0.,  5.],
1968                                                [ 5.,  6.,  1.],
1969                                                [ 2.,  1.,  6.],
1970                                                [ 6.,  7.,  2.],
1971                                                [ 3.,  2.,  7.]])
1972
1973    def test_setting_unique_vertex_values(self):
1974        """
1975        set values based on unique_vertex lists.
1976        """
1977        from mesh_factory import rectangular
1978        from shallow_water import Domain
1979        from Numeric import zeros, Float
1980
1981        #Create basic mesh
1982        points, vertices, boundary = rectangular(1, 3)
1983        #print "vertices",vertices
1984        #Create shallow water domain
1985        domain = Domain(points, vertices, boundary)
1986        #print "domain.number_of_elements ",domain.number_of_elements
1987        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1988                                    [4,4,4],[5,5,5]])
1989        value = 7
1990        indices = [1,5]
1991        quantity.set_values(value,
1992                            location = 'unique vertices',
1993                            indices = indices)
1994        #print "quantity.centroid_values",quantity.centroid_values
1995        assert allclose(quantity.vertex_values[0], [0,7,0])
1996        assert allclose(quantity.vertex_values[1], [7,1,7])
1997        assert allclose(quantity.vertex_values[2], [7,2,7])
1998
1999
2000    def test_get_values(self):
2001        """
2002        get values based on triangle lists.
2003        """
2004        from mesh_factory import rectangular
2005        from shallow_water import Domain
2006        from Numeric import zeros, Float
2007
2008        #Create basic mesh
2009        points, vertices, boundary = rectangular(1, 3)
2010
2011        #print "points",points
2012        #print "vertices",vertices
2013        #print "boundary",boundary
2014
2015        #Create shallow water domain
2016        domain = Domain(points, vertices, boundary)
2017        #print "domain.number_of_elements ",domain.number_of_elements
2018        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2019                                    [4,4,4],[5,5,5]])
2020
2021        #print "quantity.get_values(location = 'unique vertices')", \
2022        #      quantity.get_values(location = 'unique vertices')
2023
2024        #print "quantity.get_values(location = 'unique vertices')", \
2025        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2026        #                          location = 'unique vertices')
2027
2028        answer = [0.5,2,4,5,0,1,3,4.5]
2029        assert allclose(answer,
2030                        quantity.get_values(location = 'unique vertices'))
2031
2032        indices = [0,5,3]
2033        answer = [0.5,1,5]
2034        assert allclose(answer,
2035                        quantity.get_values(indices=indices, \
2036                                            location = 'unique vertices'))
2037        #print "quantity.centroid_values",quantity.centroid_values
2038        #print "quantity.get_values(location = 'centroids') ",\
2039        #      quantity.get_values(location = 'centroids')
2040
2041
2042
2043
2044    def test_get_values_2(self):
2045        """Different mesh (working with domain object) - also check centroids.
2046        """
2047
2048       
2049        a = [0.0, 0.0]
2050        b = [0.0, 2.0]
2051        c = [2.0,0.0]
2052        d = [0.0, 4.0]
2053        e = [2.0, 2.0]
2054        f = [4.0,0.0]
2055
2056        points = [a, b, c, d, e, f]
2057        #bac, bce, ecf, dbe
2058        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2059
2060        domain = Domain(points, vertices)
2061
2062        quantity = Quantity(domain)
2063        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2064       
2065        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
2066        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
2067
2068
2069        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
2070                                                                   [4,2,6],
2071                                                                   [6,2,4],
2072                                                                   [8,4,6]])
2073       
2074        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
2075                                                                                  [8,4,6]])
2076
2077
2078        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
2079                                                                [4,5,3],
2080                                                                [3,5,4],
2081                                                                [5,7,6]])
2082        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
2083                        [[4,5,3],
2084                         [5,7,6]])       
2085
2086        # Check averaging over vertices
2087        #a: 0
2088        #b: (4+4+4)/3
2089        #c: (2+2+2)/3
2090        #d: 8
2091        #e: (6+6+6)/3       
2092        #f: 4
2093        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
2094                                                                                 
2095       
2096
2097
2098
2099
2100    def test_get_interpolated_values(self):
2101
2102        from mesh_factory import rectangular
2103        from shallow_water import Domain
2104        from Numeric import zeros, Float
2105
2106        #Create basic mesh
2107        points, vertices, boundary = rectangular(1, 3)
2108        domain = Domain(points, vertices, boundary)
2109
2110        #Constant values
2111        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2112                                    [4,4,4],[5,5,5]])
2113
2114       
2115
2116        # Get interpolated values at centroids
2117        interpolation_points = domain.get_centroid_coordinates()
2118        answer = quantity.get_values(location='centroids')
2119
2120       
2121        #print quantity.get_values(points=interpolation_points)
2122        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
2123
2124
2125        #Arbitrary values
2126        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
2127                                    [1,4,-9],[2,5,0]])
2128
2129
2130        # Get interpolated values at centroids
2131        interpolation_points = domain.get_centroid_coordinates()
2132        answer = quantity.get_values(location='centroids')
2133        #print answer
2134        #print quantity.get_values(points=interpolation_points)
2135        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
2136                       
2137
2138        #FIXME TODO
2139        #indices = [0,5,3]
2140        #answer = [0.5,1,5]
2141        #assert allclose(answer,
2142        #                quantity.get_values(indices=indices, \
2143        #                                    location = 'unique vertices'))
2144
2145
2146
2147
2148    def test_get_interpolated_values_2(self):
2149        a = [0.0, 0.0]
2150        b = [0.0, 2.0]
2151        c = [2.0,0.0]
2152        d = [0.0, 4.0]
2153        e = [2.0, 2.0]
2154        f = [4.0,0.0]
2155
2156        points = [a, b, c, d, e, f]
2157        #bac, bce, ecf, dbe
2158        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
2159
2160        domain = Domain(points, vertices)
2161
2162        quantity = Quantity(domain)
2163        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
2164
2165        #First pick one point
2166        x, y = 2.0/3, 8.0/3
2167        v = quantity.get_values(interpolation_points = [[x,y]])
2168        assert allclose(v, 6)       
2169
2170        # Then another to test that algorithm won't blindly
2171        # reuse interpolation matrix
2172        x, y = 4.0/3, 4.0/3
2173        v = quantity.get_values(interpolation_points = [[x,y]])
2174        assert allclose(v, 4)       
2175
2176
2177
2178
2179    def test_getting_some_vertex_values(self):
2180        """
2181        get values based on triangle lists.
2182        """
2183        from mesh_factory import rectangular
2184        from shallow_water import Domain
2185        from Numeric import zeros, Float
2186
2187        #Create basic mesh
2188        points, vertices, boundary = rectangular(1, 3)
2189
2190        #print "points",points
2191        #print "vertices",vertices
2192        #print "boundary",boundary
2193
2194        #Create shallow water domain
2195        domain = Domain(points, vertices, boundary)
2196        #print "domain.number_of_elements ",domain.number_of_elements
2197        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
2198                                    [4,4,4],[5,5,5],[6,6,6]])
2199        value = [7]
2200        indices = [1]
2201        quantity.set_values(value,
2202                            location = 'centroids',
2203                            indices = indices)
2204        #print "quantity.centroid_values",quantity.centroid_values
2205        #print "quantity.get_values(location = 'centroids') ",\
2206        #      quantity.get_values(location = 'centroids')
2207        assert allclose(quantity.centroid_values,
2208                        quantity.get_values(location = 'centroids'))
2209
2210
2211        value = [[15,20,25]]
2212        quantity.set_values(value, indices = indices)
2213        #print "1 quantity.vertex_values",quantity.vertex_values
2214        assert allclose(quantity.vertex_values, quantity.get_values())
2215
2216        assert allclose(quantity.edge_values,
2217                        quantity.get_values(location = 'edges'))
2218
2219        # get a subset of elements
2220        subset = quantity.get_values(location='centroids', indices=[0,5])
2221        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
2222        assert allclose(subset, answer)
2223
2224
2225        subset = quantity.get_values(location='edges', indices=[0,5])
2226        answer = [quantity.edge_values[0],quantity.edge_values[5]]
2227        #print "subset",subset
2228        #print "answer",answer
2229        assert allclose(subset, answer)
2230
2231        subset = quantity.get_values( indices=[1,5])
2232        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
2233        #print "subset",subset
2234        #print "answer",answer
2235        assert allclose(subset, answer)
2236
2237    def test_smooth_vertex_values(self):
2238        """
2239        get values based on triangle lists.
2240        """
2241        from mesh_factory import rectangular
2242        from shallow_water import Domain
2243        from Numeric import zeros, Float
2244
2245        #Create basic mesh
2246        points, vertices, boundary = rectangular(2, 2)
2247
2248        #print "points",points
2249        #print "vertices",vertices
2250        #print "boundary",boundary
2251
2252        #Create shallow water domain
2253        domain = Domain(points, vertices, boundary)
2254        #print "domain.number_of_elements ",domain.number_of_elements
2255        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
2256                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
2257
2258        #print "quantity.get_values(location = 'unique vertices')", \
2259        #      quantity.get_values(location = 'unique vertices')
2260
2261        #print "quantity.get_values(location = 'unique vertices')", \
2262        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
2263        #                          location = 'unique vertices')
2264
2265        #print quantity.get_values(location = 'unique vertices')
2266        #print quantity.domain.number_of_triangles_per_node
2267        #print quantity.vertex_values
2268       
2269        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
2270        #assert allclose(answer,
2271        #                quantity.get_values(location = 'unique vertices'))
2272
2273        quantity.smooth_vertex_values()
2274
2275        #print quantity.vertex_values
2276
2277
2278        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
2279                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
2280       
2281        assert allclose(answer_vertex_values,
2282                        quantity.vertex_values)
2283        #print "quantity.centroid_values",quantity.centroid_values
2284        #print "quantity.get_values(location = 'centroids') ",\
2285        #      quantity.get_values(location = 'centroids')
2286
2287
2288
2289#-------------------------------------------------------------
2290if __name__ == "__main__":
2291    suite = unittest.makeSuite(Test_Quantity, 'test')
2292
2293    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_UTM')
2294    #print "restricted test"
2295    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
2296    runner = unittest.TextTestRunner()
2297    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.