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

Last change on this file since 4130 was 4130, checked in by duncan, 17 years ago

Added first iteration of blocking to quantity.
File size: 56.4 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    def tearDown(self):
46        pass
47        #print "  Tearing down"
48
49
50    def test_creation(self):
51
52        quantity = Quantity(self.mesh1, [[1,2,3]])
53        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
54
55        try:
56            quantity = Quantity()
57        except:
58            pass
59        else:
60            raise 'Should have raised empty quantity exception'
61
62
63        try:
64            quantity = Quantity([1,2,3])
65        except AssertionError:
66            pass
67        except:
68            raise 'Should have raised "mising mesh object" error'
69
70
71    def test_creation_zeros(self):
72
73        quantity = Quantity(self.mesh1)
74        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
75
76
77        quantity = Quantity(self.mesh4)
78        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
79                                                 [0.,0.,0.], [0.,0.,0.]])
80
81
82    def test_interpolation(self):
83        quantity = Quantity(self.mesh1, [[1,2,3]])
84        assert allclose(quantity.centroid_values, [2.0]) #Centroid
85
86        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
87
88
89    def test_interpolation2(self):
90        quantity = Conserved_quantity(self.mesh4,
91                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
92        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
93
94
95        quantity.extrapolate_second_order()
96
97        #print quantity.vertex_values
98        #assert allclose(quantity.vertex_values, [[2., 2., 2.],
99        #                                         [3.+2./3, 6.+2./3, 4.+2./3],
100        #                                         [7.5, 0.5, 1.],
101        #                                         [-5, -2.5, 7.5]])
102
103        assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
104        #FIXME: Work out the others
105
106
107        #print quantity.edge_values
108        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
109                                               [5., 5., 5.],
110                                               [4.5, 4.5, 0.],
111                                               [3.0, -1.5, -1.5]])
112
113    def test_get_maximum_1(self):
114        quantity = Conserved_quantity(self.mesh4,
115                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
116        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
117
118        v = quantity.get_maximum_value()
119        assert v == 5
120
121        i = quantity.get_maximum_index()
122        assert i == 1
123       
124        x,y = quantity.get_maximum_location()
125        xref, yref = 4.0/3, 4.0/3
126        assert x == xref
127        assert y == yref
128
129        v = quantity.get_values(interpolation_points = [[x,y]])
130        assert allclose(v, 5)
131
132    def test_get_maximum_2(self):
133
134        a = [0.0, 0.0]
135        b = [0.0, 2.0]
136        c = [2.0,0.0]
137        d = [0.0, 4.0]
138        e = [2.0, 2.0]
139        f = [4.0,0.0]
140
141        points = [a, b, c, d, e, f]
142        #bac, bce, ecf, dbe
143        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
144
145        domain = Domain(points, vertices)
146
147        quantity = Quantity(domain)
148        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
149       
150        v = quantity.get_maximum_value()
151        assert v == 6
152
153        i = quantity.get_maximum_index()
154        assert i == 3
155       
156        x,y = quantity.get_maximum_location()
157        xref, yref = 2.0/3, 8.0/3
158        assert x == xref
159        assert y == yref
160
161        v = quantity.get_values(interpolation_points = [[x,y]])
162        assert allclose(v, 6)
163
164
165        #Multiple locations for maximum -
166        #Test that the algorithm picks the first occurrence       
167        v = quantity.get_maximum_value(indices=[0,1,2])
168        assert allclose(v, 4)
169
170        i = quantity.get_maximum_index(indices=[0,1,2])
171        assert i == 1
172       
173        x,y = quantity.get_maximum_location(indices=[0,1,2])
174        xref, yref = 4.0/3, 4.0/3
175        assert x == xref
176        assert y == yref
177
178        v = quantity.get_values(interpolation_points = [[x,y]])
179        assert allclose(v, 4)       
180
181        # More test of indices......
182        v = quantity.get_maximum_value(indices=[2,3])
183        assert allclose(v, 6)
184
185        i = quantity.get_maximum_index(indices=[2,3])
186        assert i == 3
187       
188        x,y = quantity.get_maximum_location(indices=[2,3])
189        xref, yref = 2.0/3, 8.0/3
190        assert x == xref
191        assert y == yref
192
193        v = quantity.get_values(interpolation_points = [[x,y]])
194        assert allclose(v, 6)       
195
196       
197
198    def test_boundary_allocation(self):
199        quantity = Conserved_quantity(self.mesh4,
200                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
201
202        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
203
204
205    def test_set_values(self):
206        quantity = Quantity(self.mesh4)
207
208
209        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
210                            location = 'vertices')
211        assert allclose(quantity.vertex_values,
212                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
213        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
214        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
215                                               [5., 5., 5.],
216                                               [4.5, 4.5, 0.],
217                                               [3.0, -1.5, -1.5]])
218
219
220        #Test default
221        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
222        assert allclose(quantity.vertex_values,
223                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
224        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
225        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
226                                               [5., 5., 5.],
227                                               [4.5, 4.5, 0.],
228                                               [3.0, -1.5, -1.5]])
229
230        #Test centroids
231        quantity.set_values([1,2,3,4], location = 'centroids')
232        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
233
234        #Test edges
235        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
236                            location = 'edges')
237        assert allclose(quantity.edge_values,
238                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
239
240        #Test exceptions
241        try:
242            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
243                                location = 'bas kamel tuba')
244        except:
245            pass
246
247
248        try:
249            quantity.set_values([[1,2,3], [0,0,9]])
250        except AssertionError:
251            pass
252        except:
253            raise 'should have raised Assertionerror'
254
255
256
257    def test_set_values_const(self):
258        quantity = Quantity(self.mesh4)
259
260        quantity.set_values(1.0, location = 'vertices')
261        assert allclose(quantity.vertex_values,
262                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
263
264        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
265        assert allclose(quantity.edge_values, [[1, 1, 1],
266                                               [1, 1, 1],
267                                               [1, 1, 1],
268                                               [1, 1, 1]])
269
270
271        quantity.set_values(2.0, location = 'centroids')
272        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
273
274        quantity.set_values(3.0, location = 'edges')
275        assert allclose(quantity.edge_values, [[3, 3, 3],
276                                               [3, 3, 3],
277                                               [3, 3, 3],
278                                               [3, 3, 3]])
279
280
281    def test_set_values_func(self):
282        quantity = Quantity(self.mesh4)
283
284        def f(x, y):
285            return x+y
286
287        quantity.set_values(f, location = 'vertices')
288        #print "quantity.vertex_values",quantity.vertex_values
289        assert allclose(quantity.vertex_values,
290                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
291        assert allclose(quantity.centroid_values,
292                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
293        assert allclose(quantity.edge_values,
294                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
295
296
297        quantity.set_values(f, location = 'centroids')
298        assert allclose(quantity.centroid_values,
299                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
300
301
302    def test_integral(self):
303        quantity = Quantity(self.mesh4)
304
305        #Try constants first
306        const = 5
307        quantity.set_values(const, location = 'vertices')
308        #print 'Q', quantity.get_integral()
309
310        assert allclose(quantity.get_integral(), self.mesh4.get_area() * const)
311
312        #Try with a linear function
313        def f(x, y):
314            return x+y
315
316        quantity.set_values(f, location = 'vertices')
317
318
319        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
320
321        assert allclose (quantity.get_integral(), ref_integral)
322
323
324
325    def test_set_vertex_values(self):
326        quantity = Quantity(self.mesh4)
327        quantity.set_vertex_values([0,1,2,3,4,5])
328
329        assert allclose(quantity.vertex_values,
330                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
331        assert allclose(quantity.centroid_values,
332                        [1., 7./3, 11./3, 8./3]) #Centroid
333        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
334                                               [3., 2.5, 1.5],
335                                               [3.5, 4.5, 3.],
336                                               [2.5, 3.5, 2]])
337
338
339    def test_set_vertex_values_subset(self):
340        quantity = Quantity(self.mesh4)
341        quantity.set_vertex_values([0,1,2,3,4,5])
342        quantity.set_vertex_values([0,20,30,50], indices = [0,2,3,5])
343
344        assert allclose(quantity.vertex_values,
345                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
346
347
348    def test_set_vertex_values_using_general_interface(self):
349        quantity = Quantity(self.mesh4)
350
351
352        quantity.set_values([0,1,2,3,4,5])
353
354
355        assert allclose(quantity.vertex_values,
356                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
357
358        #Centroid
359        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
360
361        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
362                                               [3., 2.5, 1.5],
363                                               [3.5, 4.5, 3.],
364                                               [2.5, 3.5, 2]])
365
366
367
368
369
370    def test_set_values_using_fit(self):
371
372
373        quantity = Quantity(self.mesh4)
374
375        #Get (enough) datapoints
376        data_points = [[ 0.66666667, 0.66666667],
377                       [ 1.33333333, 1.33333333],
378                       [ 2.66666667, 0.66666667],
379                       [ 0.66666667, 2.66666667],
380                       [ 0.0, 1.0],
381                       [ 0.0, 3.0],
382                       [ 1.0, 0.0],
383                       [ 1.0, 1.0],
384                       [ 1.0, 2.0],
385                       [ 1.0, 3.0],
386                       [ 2.0, 1.0],
387                       [ 3.0, 0.0],
388                       [ 3.0, 1.0]]
389
390        z = linear_function(data_points)
391
392        #Use built-in fit_interpolate.fit
393        quantity.set_values( Geospatial_data(data_points, z), alpha = 0 )
394        #quantity.set_values(points = data_points, values = z, alpha = 0)
395
396
397        answer = linear_function(quantity.domain.get_vertex_coordinates())
398        #print quantity.vertex_values, answer
399        assert allclose(quantity.vertex_values.flat, answer)
400
401
402        #Now try by setting the same values directly
403        vertex_attributes = fit_to_mesh(quantity.domain.get_nodes(),
404                                        quantity.domain.triangles, #FIXME
405                                        data_points,
406                                        z,
407                                        alpha = 0,
408                                        verbose=False)
409
410        #print vertex_attributes
411        quantity.set_values(vertex_attributes)
412        assert allclose(quantity.vertex_values.flat, answer)
413
414
415
416
417
418    def test_test_set_values_using_fit_w_geo(self):
419
420
421        #Mesh
422        vertex_coordinates = [[0.76, 0.76],
423                              [0.76, 5.76],
424                              [5.76, 0.76]]
425        triangles = [[0,2,1]]
426
427        mesh_georef = Geo_reference(56,-0.76,-0.76)
428        mesh1 = Domain(vertex_coordinates, triangles,
429                       geo_reference = mesh_georef)
430        mesh1.check_integrity()
431
432        #Quantity
433        quantity = Quantity(mesh1)
434
435        #Data
436        data_points = [[ 201.0, 401.0],
437                       [ 201.0, 403.0],
438                       [ 203.0, 401.0]]
439
440        z = [2, 4, 4]
441
442        data_georef = Geo_reference(56,-200,-400)
443
444
445        #Reference
446        ref = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
447                          data_origin = data_georef.get_origin(),
448                          mesh_origin = mesh_georef.get_origin(),
449                          alpha = 0)
450
451        assert allclose( ref, [0,5,5] )
452
453
454        #Test set_values
455
456        quantity.set_values( Geospatial_data(data_points, z, data_georef), alpha = 0 )
457
458        #quantity.set_values(points = data_points,
459        #                    values = z,
460        #                    data_georef = data_georef,
461        #                    alpha = 0)
462
463
464        #quantity.set_values(points = data_points,
465        #                    values = z,
466        #                    data_georef = data_georef,
467        #                    alpha = 0)
468        assert allclose(quantity.vertex_values.flat, ref)
469
470
471
472        #Test set_values using geospatial data object
473        quantity.vertex_values[:] = 0.0
474
475        geo = Geospatial_data(data_points, z, data_georef)
476
477
478        quantity.set_values(geospatial_data = geo, alpha = 0)
479        assert allclose(quantity.vertex_values.flat, ref)
480
481
482
483    def test_set_values_from_file1(self):
484        quantity = Quantity(self.mesh4)
485
486        #Get (enough) datapoints
487        data_points = [[ 0.66666667, 0.66666667],
488                       [ 1.33333333, 1.33333333],
489                       [ 2.66666667, 0.66666667],
490                       [ 0.66666667, 2.66666667],
491                       [ 0.0, 1.0],
492                       [ 0.0, 3.0],
493                       [ 1.0, 0.0],
494                       [ 1.0, 1.0],
495                       [ 1.0, 2.0],
496                       [ 1.0, 3.0],
497                       [ 2.0, 1.0],
498                       [ 3.0, 0.0],
499                       [ 3.0, 1.0]]
500
501        data_geo_spatial = Geospatial_data(data_points,
502                         geo_reference = Geo_reference(56, 0, 0))
503        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
504        attributes = linear_function(data_points_absolute)
505        att = 'spam_and_eggs'
506       
507        #Create .txt file
508        ptsfile = tempfile.mktemp(".txt")
509        file = open(ptsfile,"w")
510        file.write(" x,y," + att + " \n")
511        for data_point, attribute in map(None, data_points_absolute
512                                         ,attributes):
513            row = str(data_point[0]) + ',' + str(data_point[1]) \
514                  + ',' + str(attribute)
515            file.write(row + "\n")
516        file.close()
517
518
519        #Check that values can be set from file
520        quantity.set_values(filename = ptsfile,
521                            attribute_name = att, alpha = 0)
522        answer = linear_function(quantity.domain.get_vertex_coordinates())
523
524        #print quantity.vertex_values.flat
525        #print answer
526
527
528        assert allclose(quantity.vertex_values.flat, answer)
529
530
531        #Check that values can be set from file using default attribute
532        quantity.set_values(filename = ptsfile, alpha = 0)
533        assert allclose(quantity.vertex_values.flat, answer)
534
535        #Cleanup
536        import os
537        os.remove(ptsfile)
538
539
540    def test_set_values_from_file_with_georef1(self):
541
542        #Mesh in zone 56 (absolute coords)
543
544        x0 = 314036.58727982
545        y0 = 6224951.2960092
546
547        a = [x0+0.0, y0+0.0]
548        b = [x0+0.0, y0+2.0]
549        c = [x0+2.0, y0+0.0]
550        d = [x0+0.0, y0+4.0]
551        e = [x0+2.0, y0+2.0]
552        f = [x0+4.0, y0+0.0]
553
554        points = [a, b, c, d, e, f]
555
556        #bac, bce, ecf, dbe
557        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
558
559        #absolute going in ..
560        mesh4 = Domain(points, elements,
561                       geo_reference = Geo_reference(56, 0, 0))
562        mesh4.check_integrity()
563        quantity = Quantity(mesh4)
564
565        #Get (enough) datapoints (relative to georef)
566        data_points_rel = [[ 0.66666667, 0.66666667],
567                       [ 1.33333333, 1.33333333],
568                       [ 2.66666667, 0.66666667],
569                       [ 0.66666667, 2.66666667],
570                       [ 0.0, 1.0],
571                       [ 0.0, 3.0],
572                       [ 1.0, 0.0],
573                       [ 1.0, 1.0],
574                       [ 1.0, 2.0],
575                       [ 1.0, 3.0],
576                       [ 2.0, 1.0],
577                       [ 3.0, 0.0],
578                       [ 3.0, 1.0]]
579
580        data_geo_spatial = Geospatial_data(data_points_rel,
581                         geo_reference = Geo_reference(56, x0, y0))
582        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
583        attributes = linear_function(data_points_absolute)
584        att = 'spam_and_eggs'
585       
586        #Create .txt file
587        ptsfile = tempfile.mktemp(".txt")
588        file = open(ptsfile,"w")
589        file.write(" x,y," + att + " \n")
590        for data_point, attribute in map(None, data_points_absolute
591                                         ,attributes):
592            row = str(data_point[0]) + ',' + str(data_point[1]) \
593                  + ',' + str(attribute)
594            file.write(row + "\n")
595        file.close()
596
597        #file = open(ptsfile, 'r')
598        #lines = file.readlines()
599        #file.close()
600     
601
602        #Check that values can be set from file
603        quantity.set_values(filename = ptsfile,
604                            attribute_name = att, alpha = 0)
605        answer = linear_function(quantity.domain.get_vertex_coordinates())
606
607        assert allclose(quantity.vertex_values.flat, answer)
608
609
610        #Check that values can be set from file using default attribute
611        quantity.set_values(filename = ptsfile, alpha = 0)
612        assert allclose(quantity.vertex_values.flat, answer)
613
614        #Cleanup
615        import os
616        os.remove(ptsfile)
617
618
619    def test_set_values_from_file_with_georef2(self):
620
621        #Mesh in zone 56 (relative coords)
622
623        x0 = 314036.58727982
624        y0 = 6224951.2960092
625
626        a = [0.0, 0.0]
627        b = [0.0, 2.0]
628        c = [2.0, 0.0]
629        d = [0.0, 4.0]
630        e = [2.0, 2.0]
631        f = [4.0, 0.0]
632
633        points = [a, b, c, d, e, f]
634
635        #bac, bce, ecf, dbe
636        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
637
638        mesh4 = Domain(points, elements,
639                       geo_reference = Geo_reference(56, x0, y0))
640        mesh4.check_integrity()
641        quantity = Quantity(mesh4)
642
643        #Get (enough) datapoints
644        data_points = [[ x0+0.66666667, y0+0.66666667],
645                       [ x0+1.33333333, y0+1.33333333],
646                       [ x0+2.66666667, y0+0.66666667],
647                       [ x0+0.66666667, y0+2.66666667],
648                       [ x0+0.0, y0+1.0],
649                       [ x0+0.0, y0+3.0],
650                       [ x0+1.0, y0+0.0],
651                       [ x0+1.0, y0+1.0],
652                       [ x0+1.0, y0+2.0],
653                       [ x0+1.0, y0+3.0],
654                       [ x0+2.0, y0+1.0],
655                       [ x0+3.0, y0+0.0],
656                       [ x0+3.0, y0+1.0]]
657
658
659        data_geo_spatial = Geospatial_data(data_points,
660                         geo_reference = Geo_reference(56, 0, 0))
661        data_points_absolute = data_geo_spatial.get_data_points(absolute=True)
662        attributes = linear_function(data_points_absolute)
663        att = 'spam_and_eggs'
664       
665        #Create .txt file
666        ptsfile = tempfile.mktemp(".txt")
667        file = open(ptsfile,"w")
668        file.write(" x,y," + att + " \n")
669        for data_point, attribute in map(None, data_points_absolute
670                                         ,attributes):
671            row = str(data_point[0]) + ',' + str(data_point[1]) \
672                  + ',' + str(attribute)
673            file.write(row + "\n")
674        file.close()
675
676
677        #Check that values can be set from file
678        quantity.set_values(filename = ptsfile,
679                            attribute_name = att, alpha = 0)
680        answer = linear_function(quantity.domain. \
681                                 get_vertex_coordinates(absolute=True))
682
683
684        assert allclose(quantity.vertex_values.flat, answer)
685
686
687        #Check that values can be set from file using default attribute
688        quantity.set_values(filename = ptsfile, alpha = 0)
689        assert allclose(quantity.vertex_values.flat, answer)
690
691        #Cleanup
692        import os
693        os.remove(ptsfile)
694
695
696
697
698    def test_set_values_from_quantity(self):
699
700        quantity1 = Quantity(self.mesh4)
701        quantity1.set_vertex_values([0,1,2,3,4,5])
702
703        assert allclose(quantity1.vertex_values,
704                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
705
706
707        quantity2 = Quantity(self.mesh4)
708        quantity2.set_values(quantity = quantity1)
709        assert allclose(quantity2.vertex_values,
710                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
711
712        quantity2.set_values(quantity = 2*quantity1)
713        assert allclose(quantity2.vertex_values,
714                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
715
716        quantity2.set_values(quantity = 2*quantity1 + 3)
717        assert allclose(quantity2.vertex_values,
718                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
719
720
721        #Check detection of quantity as first orgument
722        quantity2.set_values(2*quantity1 + 3)
723        assert allclose(quantity2.vertex_values,
724                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
725
726
727
728
729
730    def test_overloading(self):
731
732        quantity1 = Quantity(self.mesh4)
733        quantity1.set_vertex_values([0,1,2,3,4,5])
734
735        assert allclose(quantity1.vertex_values,
736                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
737
738
739        quantity2 = Quantity(self.mesh4)
740        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
741                             location = 'vertices')
742
743
744
745        quantity3 = Quantity(self.mesh4)
746        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
747                             location = 'vertices')
748
749
750        #Negation
751        Q = -quantity1
752        assert allclose(Q.vertex_values, -quantity1.vertex_values)
753        assert allclose(Q.centroid_values, -quantity1.centroid_values)
754        assert allclose(Q.edge_values, -quantity1.edge_values)
755
756        #Addition
757        Q = quantity1 + 7
758        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
759        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
760        assert allclose(Q.edge_values, quantity1.edge_values + 7)
761
762        Q = 7 + quantity1
763        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
764        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
765        assert allclose(Q.edge_values, quantity1.edge_values + 7)
766
767        Q = quantity1 + quantity2
768        assert allclose(Q.vertex_values,
769                        quantity1.vertex_values + quantity2.vertex_values)
770        assert allclose(Q.centroid_values,
771                        quantity1.centroid_values + quantity2.centroid_values)
772        assert allclose(Q.edge_values,
773                        quantity1.edge_values + quantity2.edge_values)
774
775
776        Q = quantity1 + quantity2 - 3
777        assert allclose(Q.vertex_values,
778                        quantity1.vertex_values + quantity2.vertex_values - 3)
779
780        Q = quantity1 - quantity2
781        assert allclose(Q.vertex_values,
782                        quantity1.vertex_values - quantity2.vertex_values)
783
784        #Scaling
785        Q = quantity1*3
786        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
787        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
788        assert allclose(Q.edge_values, quantity1.edge_values*3)
789        Q = 3*quantity1
790        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
791
792        #Multiplication
793        Q = quantity1 * quantity2
794        #print Q.vertex_values
795        #print Q.centroid_values
796        #print quantity1.centroid_values
797        #print quantity2.centroid_values
798
799        assert allclose(Q.vertex_values,
800                        quantity1.vertex_values * quantity2.vertex_values)
801
802        #Linear combinations
803        Q = 4*quantity1 + 2
804        assert allclose(Q.vertex_values,
805                        4*quantity1.vertex_values + 2)
806
807        Q = quantity1*quantity2 + 2
808        assert allclose(Q.vertex_values,
809                        quantity1.vertex_values * quantity2.vertex_values + 2)
810
811        Q = quantity1*quantity2 + quantity3
812        assert allclose(Q.vertex_values,
813                        quantity1.vertex_values * quantity2.vertex_values +
814                        quantity3.vertex_values)
815        Q = quantity1*quantity2 + 3*quantity3
816        assert allclose(Q.vertex_values,
817                        quantity1.vertex_values * quantity2.vertex_values +
818                        3*quantity3.vertex_values)
819        Q = quantity1*quantity2 + 3*quantity3 + 5.0
820        assert allclose(Q.vertex_values,
821                        quantity1.vertex_values * quantity2.vertex_values +
822                        3*quantity3.vertex_values + 5)
823
824        Q = quantity1*quantity2 - quantity3
825        assert allclose(Q.vertex_values,
826                        quantity1.vertex_values * quantity2.vertex_values -
827                        quantity3.vertex_values)
828        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
829        assert allclose(Q.vertex_values,
830                        1.5*quantity1.vertex_values * quantity2.vertex_values -
831                        3*quantity3.vertex_values + 5)
832
833        #Try combining quantities and arrays and scalars
834        Q = 1.5*quantity1*quantity2.vertex_values -\
835            3*quantity3.vertex_values + 5.0
836        assert allclose(Q.vertex_values,
837                        1.5*quantity1.vertex_values * quantity2.vertex_values -
838                        3*quantity3.vertex_values + 5)
839
840
841        #Powers
842        Q = quantity1**2
843        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
844
845        Q = quantity1**2 +quantity2**2
846        assert allclose(Q.vertex_values,
847                        quantity1.vertex_values**2 + \
848                        quantity2.vertex_values**2)
849
850        Q = (quantity1**2 +quantity2**2)**0.5
851        assert allclose(Q.vertex_values,
852                        (quantity1.vertex_values**2 + \
853                        quantity2.vertex_values**2)**0.5)
854
855
856
857
858
859
860
861    def test_gradient(self):
862        quantity = Conserved_quantity(self.mesh4)
863
864        #Set up for a gradient of (2,0) at mid triangle
865        quantity.set_values([2.0, 4.0, 6.0, 2.0],
866                            location = 'centroids')
867
868
869        #Gradients
870        a, b = quantity.compute_gradients()
871
872        #print self.mesh4.centroid_coordinates
873        #print a, b
874
875        #The central triangle (1)
876        #(using standard gradient based on neigbours controid values)
877        assert allclose(a[1], 2.0)
878        assert allclose(b[1], 0.0)
879
880
881        #Left triangle (0) using two point gradient
882        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
883        #2  = 4  + a*(-2/3)  + b*(-2/3)
884        assert allclose(a[0] + b[0], 3)
885        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
886        assert allclose(a[0] - b[0], 0)
887
888
889        #Right triangle (2) using two point gradient
890        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
891        #6  = 4  + a*(4/3)  + b*(-2/3)
892        assert allclose(2*a[2] - b[2], 3)
893        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
894        assert allclose(a[2] + 2*b[2], 0)
895
896
897        #Top triangle (3) using two point gradient
898        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
899        #2  = 4  + a*(-2/3)  + b*(4/3)
900        assert allclose(a[3] - 2*b[3], 3)
901        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
902        assert allclose(2*a[3] + b[3], 0)
903
904
905
906        #print a, b
907        quantity.extrapolate_second_order()
908
909        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
910        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
911        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
912
913
914        #a = 1.2, b=-0.6
915        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
916        assert allclose(quantity.vertex_values[2,2], 8)
917
918
919    def test_second_order_extrapolation2(self):
920        quantity = Conserved_quantity(self.mesh4)
921
922        #Set up for a gradient of (3,1), f(x) = 3x+y
923        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
924                            location = 'centroids')
925
926        #Gradients
927        a, b = quantity.compute_gradients()
928
929        #print a, b
930
931        assert allclose(a[1], 3.0)
932        assert allclose(b[1], 1.0)
933
934        #Work out the others
935
936        quantity.extrapolate_second_order()
937
938        #print quantity.vertex_values
939        assert allclose(quantity.vertex_values[1,0], 2.0)
940        assert allclose(quantity.vertex_values[1,1], 6.0)
941        assert allclose(quantity.vertex_values[1,2], 8.0)
942
943
944
945    def test_first_order_extrapolator(self):
946        quantity = Conserved_quantity(self.mesh4)
947
948        #Test centroids
949        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
950        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
951
952        #Extrapolate
953        quantity.extrapolate_first_order()
954
955        #Check vertices but not edge values
956        assert allclose(quantity.vertex_values,
957                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
958
959
960    def test_second_order_extrapolator(self):
961        quantity = Conserved_quantity(self.mesh4)
962
963        #Set up for a gradient of (3,0) at mid triangle
964        quantity.set_values([2.0, 4.0, 8.0, 2.0],
965                            location = 'centroids')
966
967
968
969        quantity.extrapolate_second_order()
970        quantity.limit()
971
972
973        #Assert that central triangle is limited by neighbours
974        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
975        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
976
977        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
978        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
979
980        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
981        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
982
983
984        #Assert that quantities are conserved
985        from Numeric import sum
986        for k in range(quantity.centroid_values.shape[0]):
987            assert allclose (quantity.centroid_values[k],
988                             sum(quantity.vertex_values[k,:])/3)
989
990
991
992
993
994    def test_limiter(self):
995        quantity = Conserved_quantity(self.mesh4)
996
997        #Create a deliberate overshoot (e.g. from gradient computation)
998        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
999
1000
1001        #Limit
1002        quantity.limit()
1003
1004        #Assert that central triangle is limited by neighbours
1005        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
1006        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
1007
1008        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
1009        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
1010
1011        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
1012        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
1013
1014
1015
1016        #Assert that quantities are conserved
1017        from Numeric import sum
1018        for k in range(quantity.centroid_values.shape[0]):
1019            assert allclose (quantity.centroid_values[k],
1020                             sum(quantity.vertex_values[k,:])/3)
1021
1022
1023    def test_limiter2(self):
1024        """Taken from test_shallow_water
1025        """
1026        quantity = Conserved_quantity(self.mesh4)
1027        quantity.domain.beta_w = 0.9
1028       
1029        #Test centroids
1030        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1031        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1032
1033
1034        #Extrapolate
1035        quantity.extrapolate_second_order()
1036
1037        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1038
1039        #Limit
1040        quantity.limit()
1041
1042        # limited value for beta_w = 0.9
1043       
1044        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
1045        # limited values for beta_w = 0.5
1046        #assert allclose(quantity.vertex_values[1,:], [3.0, 4.5, 4.5])
1047
1048
1049        #Assert that quantities are conserved
1050        from Numeric import sum
1051        for k in range(quantity.centroid_values.shape[0]):
1052            assert allclose (quantity.centroid_values[k],
1053                             sum(quantity.vertex_values[k,:])/3)
1054
1055
1056
1057
1058
1059    def test_distribute_first_order(self):
1060        quantity = Conserved_quantity(self.mesh4)
1061
1062        #Test centroids
1063        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1064        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1065
1066
1067        #Extrapolate
1068        quantity.extrapolate_first_order()
1069
1070        #Interpolate
1071        quantity.interpolate_from_vertices_to_edges()
1072
1073        assert allclose(quantity.vertex_values,
1074                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
1075        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
1076                                               [3,3,3], [4, 4, 4]])
1077
1078
1079    def test_distribute_second_order(self):
1080        quantity = Conserved_quantity(self.mesh4)
1081
1082        #Test centroids
1083        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
1084        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
1085
1086
1087        #Extrapolate
1088        quantity.extrapolate_second_order()
1089
1090        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
1091
1092
1093    def test_update_explicit(self):
1094        quantity = Conserved_quantity(self.mesh4)
1095
1096        #Test centroids
1097        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1098        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1099
1100        #Set explicit_update
1101        quantity.explicit_update = array( [1.,1.,1.,1.] )
1102
1103        #Update with given timestep
1104        quantity.update(0.1)
1105
1106        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1107        assert allclose( quantity.centroid_values, x)
1108
1109    def test_update_semi_implicit(self):
1110        quantity = Conserved_quantity(self.mesh4)
1111
1112        #Test centroids
1113        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1114        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1115
1116        #Set semi implicit update
1117        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1118
1119        #Update with given timestep
1120        timestep = 0.1
1121        quantity.update(timestep)
1122
1123        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1124        denom = ones(4, Float)-timestep*sem
1125
1126        x = array([1, 2, 3, 4])/denom
1127        assert allclose( quantity.centroid_values, x)
1128
1129
1130    def test_both_updates(self):
1131        quantity = Conserved_quantity(self.mesh4)
1132
1133        #Test centroids
1134        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1135        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1136
1137        #Set explicit_update
1138        quantity.explicit_update = array( [4.,3.,2.,1.] )
1139
1140        #Set semi implicit update
1141        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1142
1143        #Update with given timestep
1144        timestep = 0.1
1145        quantity.update(0.1)
1146
1147        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1148        denom = ones(4, Float)-timestep*sem
1149
1150        x = array([1., 2., 3., 4.])
1151        x /= denom
1152        x += timestep*array( [4.0, 3.0, 2.0, 1.0] )
1153
1154        assert allclose( quantity.centroid_values, x)
1155
1156
1157
1158
1159    #Test smoothing
1160    def test_smoothing(self):
1161
1162        from mesh_factory import rectangular
1163        from shallow_water import Domain, Transmissive_boundary
1164        from Numeric import zeros, Float
1165        from anuga.utilities.numerical_tools import mean
1166
1167        #Create basic mesh
1168        points, vertices, boundary = rectangular(2, 2)
1169
1170        #Create shallow water domain
1171        domain = Domain(points, vertices, boundary)
1172        domain.default_order=2
1173        domain.reduction = mean
1174
1175
1176        #Set some field values
1177        domain.set_quantity('elevation', lambda x,y: x)
1178        domain.set_quantity('friction', 0.03)
1179
1180
1181        ######################
1182        # Boundary conditions
1183        B = Transmissive_boundary(domain)
1184        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1185
1186
1187        ######################
1188        #Initial condition - with jumps
1189
1190        bed = domain.quantities['elevation'].vertex_values
1191        stage = zeros(bed.shape, Float)
1192
1193        h = 0.03
1194        for i in range(stage.shape[0]):
1195            if i % 2 == 0:
1196                stage[i,:] = bed[i,:] + h
1197            else:
1198                stage[i,:] = bed[i,:]
1199
1200        domain.set_quantity('stage', stage)
1201
1202        stage = domain.quantities['stage']
1203
1204        #Get smoothed stage
1205        A, V = stage.get_vertex_values(xy=False, smooth=True)
1206        Q = stage.vertex_values
1207
1208
1209        assert A.shape[0] == 9
1210        assert V.shape[0] == 8
1211        assert V.shape[1] == 3
1212
1213        #First four points
1214        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1215        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1216        assert allclose(A[2], Q[3,0])
1217        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1218
1219        #Center point
1220        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1221                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1222
1223
1224        #Check V
1225        assert allclose(V[0,:], [3,4,0])
1226        assert allclose(V[1,:], [1,0,4])
1227        assert allclose(V[2,:], [4,5,1])
1228        assert allclose(V[3,:], [2,1,5])
1229        assert allclose(V[4,:], [6,7,3])
1230        assert allclose(V[5,:], [4,3,7])
1231        assert allclose(V[6,:], [7,8,4])
1232        assert allclose(V[7,:], [5,4,8])
1233
1234        #Get smoothed stage with XY
1235        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1236
1237        assert allclose(A, A1)
1238        assert allclose(V, V1)
1239
1240        #Check XY
1241        assert allclose(X[4], 0.5)
1242        assert allclose(Y[4], 0.5)
1243
1244        assert allclose(X[7], 1.0)
1245        assert allclose(Y[7], 0.5)
1246
1247
1248
1249
1250    def test_vertex_values_no_smoothing(self):
1251
1252        from mesh_factory import rectangular
1253        from shallow_water import Domain, Transmissive_boundary
1254        from Numeric import zeros, Float
1255        from anuga.utilities.numerical_tools import mean
1256
1257
1258        #Create basic mesh
1259        points, vertices, boundary = rectangular(2, 2)
1260
1261        #Create shallow water domain
1262        domain = Domain(points, vertices, boundary)
1263        domain.default_order=2
1264        domain.reduction = mean
1265
1266
1267        #Set some field values
1268        domain.set_quantity('elevation', lambda x,y: x)
1269        domain.set_quantity('friction', 0.03)
1270
1271
1272        ######################
1273        #Initial condition - with jumps
1274
1275        bed = domain.quantities['elevation'].vertex_values
1276        stage = zeros(bed.shape, Float)
1277
1278        h = 0.03
1279        for i in range(stage.shape[0]):
1280            if i % 2 == 0:
1281                stage[i,:] = bed[i,:] + h
1282            else:
1283                stage[i,:] = bed[i,:]
1284
1285        domain.set_quantity('stage', stage)
1286
1287        #Get stage
1288        stage = domain.quantities['stage']
1289        A, V = stage.get_vertex_values(xy=False, smooth=False)
1290        Q = stage.vertex_values.flat
1291
1292        for k in range(8):
1293            assert allclose(A[k], Q[k])
1294
1295
1296        for k in range(8):
1297            assert V[k, 0] == 3*k
1298            assert V[k, 1] == 3*k+1
1299            assert V[k, 2] == 3*k+2
1300
1301
1302
1303        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1304
1305
1306        assert allclose(A, A1)
1307        assert allclose(V, V1)
1308
1309        #Check XY
1310        assert allclose(X[1], 0.5)
1311        assert allclose(Y[1], 0.5)
1312        assert allclose(X[4], 0.0)
1313        assert allclose(Y[4], 0.0)
1314        assert allclose(X[12], 1.0)
1315        assert allclose(Y[12], 0.0)
1316
1317
1318
1319    def set_array_values_by_index(self):
1320
1321        from mesh_factory import rectangular
1322        from shallow_water import Domain
1323        from Numeric import zeros, Float
1324
1325        #Create basic mesh
1326        points, vertices, boundary = rectangular(1, 1)
1327
1328        #Create shallow water domain
1329        domain = Domain(points, vertices, boundary)
1330        #print "domain.number_of_elements ",domain.number_of_elements
1331        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1332        value = [7]
1333        indices = [1]
1334        quantity.set_array_values_by_index(value,
1335                                           location = 'centroids',
1336                                           indices = indices)
1337        #print "quantity.centroid_values",quantity.centroid_values
1338
1339        assert allclose(quantity.centroid_values, [1,7])
1340
1341        quantity.set_array_values([15,20,25], indices = indices)
1342        assert allclose(quantity.centroid_values, [1,20])
1343
1344        quantity.set_array_values([15,20,25], indices = indices)
1345        assert allclose(quantity.centroid_values, [1,20])
1346
1347    def test_setting_some_vertex_values(self):
1348        """
1349        set values based on triangle lists.
1350        """
1351        from mesh_factory import rectangular
1352        from shallow_water import Domain
1353        from Numeric import zeros, Float
1354
1355        #Create basic mesh
1356        points, vertices, boundary = rectangular(1, 3)
1357        #print "vertices",vertices
1358        #Create shallow water domain
1359        domain = Domain(points, vertices, boundary)
1360        #print "domain.number_of_elements ",domain.number_of_elements
1361        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1362                                    [4,4,4],[5,5,5],[6,6,6]])
1363
1364
1365        # Check that constants work
1366        value = 7
1367        indices = [1]
1368        quantity.set_values(value,
1369                            location = 'centroids',
1370                            indices = indices)
1371        #print "quantity.centroid_values",quantity.centroid_values
1372        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1373       
1374        value = [7]
1375        indices = [1]
1376        quantity.set_values(value,
1377                            location = 'centroids',
1378                            indices = indices)
1379        #print "quantity.centroid_values",quantity.centroid_values
1380        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1381
1382        value = [[15,20,25]]
1383        quantity.set_values(value, indices = indices)
1384        #print "1 quantity.vertex_values",quantity.vertex_values
1385        assert allclose(quantity.vertex_values[1], value[0])
1386
1387
1388        #print "quantity",quantity.vertex_values
1389        values = [10,100,50]
1390        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1391        #print "2 quantity.vertex_values",quantity.vertex_values
1392        assert allclose(quantity.vertex_values[0], [10,10,10])
1393        assert allclose(quantity.vertex_values[5], [50,50,50])
1394        #quantity.interpolate()
1395        #print "quantity.centroid_values",quantity.centroid_values
1396        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1397
1398
1399        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1400                                    [4,4,4],[5,5,5],[6,6,6]])
1401        values = [10,100,50]
1402        #this will be per unique vertex, indexing the vertices
1403        #print "quantity.vertex_values",quantity.vertex_values
1404        quantity.set_values(values, indices = [0,1,5])
1405        #print "quantity.vertex_values",quantity.vertex_values
1406        assert allclose(quantity.vertex_values[0], [1,50,10])
1407        assert allclose(quantity.vertex_values[5], [6,6,6])
1408        assert allclose(quantity.vertex_values[1], [100,10,50])
1409
1410        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1411                                    [4,4,4],[5,5,5],[6,6,6]])
1412        values = [[31,30,29],[400,400,400],[1000,999,998]]
1413        quantity.set_values(values, indices = [3,3,5])
1414        quantity.interpolate()
1415        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1416
1417        values = [[1,1,1],[2,2,2],[3,3,3],
1418                                    [4,4,4],[5,5,5],[6,6,6]]
1419        quantity.set_values(values)
1420
1421        # testing the standard set values by vertex
1422        # indexed by vertex_id in general_mesh.coordinates
1423        values = [0,1,2,3,4,5,6,7]
1424
1425        quantity.set_values(values)
1426        #print "1 quantity.vertex_values",quantity.vertex_values
1427        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1428                                                [ 1.,  0.,  5.],
1429                                                [ 5.,  6.,  1.],
1430                                                [ 2.,  1.,  6.],
1431                                                [ 6.,  7.,  2.],
1432                                                [ 3.,  2.,  7.]])
1433
1434    def test_setting_unique_vertex_values(self):
1435        """
1436        set values based on unique_vertex lists.
1437        """
1438        from mesh_factory import rectangular
1439        from shallow_water import Domain
1440        from Numeric import zeros, Float
1441
1442        #Create basic mesh
1443        points, vertices, boundary = rectangular(1, 3)
1444        #print "vertices",vertices
1445        #Create shallow water domain
1446        domain = Domain(points, vertices, boundary)
1447        #print "domain.number_of_elements ",domain.number_of_elements
1448        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1449                                    [4,4,4],[5,5,5]])
1450        value = 7
1451        indices = [1,5]
1452        quantity.set_values(value,
1453                            location = 'unique vertices',
1454                            indices = indices)
1455        #print "quantity.centroid_values",quantity.centroid_values
1456        assert allclose(quantity.vertex_values[0], [0,7,0])
1457        assert allclose(quantity.vertex_values[1], [7,1,7])
1458        assert allclose(quantity.vertex_values[2], [7,2,7])
1459
1460
1461    def test_get_values(self):
1462        """
1463        get values based on triangle lists.
1464        """
1465        from mesh_factory import rectangular
1466        from shallow_water import Domain
1467        from Numeric import zeros, Float
1468
1469        #Create basic mesh
1470        points, vertices, boundary = rectangular(1, 3)
1471
1472        #print "points",points
1473        #print "vertices",vertices
1474        #print "boundary",boundary
1475
1476        #Create shallow water domain
1477        domain = Domain(points, vertices, boundary)
1478        #print "domain.number_of_elements ",domain.number_of_elements
1479        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1480                                    [4,4,4],[5,5,5]])
1481
1482        #print "quantity.get_values(location = 'unique vertices')", \
1483        #      quantity.get_values(location = 'unique vertices')
1484
1485        #print "quantity.get_values(location = 'unique vertices')", \
1486        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1487        #                          location = 'unique vertices')
1488
1489        answer = [0.5,2,4,5,0,1,3,4.5]
1490        assert allclose(answer,
1491                        quantity.get_values(location = 'unique vertices'))
1492
1493        indices = [0,5,3]
1494        answer = [0.5,1,5]
1495        assert allclose(answer,
1496                        quantity.get_values(indices=indices, \
1497                                            location = 'unique vertices'))
1498        #print "quantity.centroid_values",quantity.centroid_values
1499        #print "quantity.get_values(location = 'centroids') ",\
1500        #      quantity.get_values(location = 'centroids')
1501
1502
1503
1504
1505    def test_get_values_2(self):
1506        """Different mesh (working with domain object) - also check centroids.
1507        """
1508
1509       
1510        a = [0.0, 0.0]
1511        b = [0.0, 2.0]
1512        c = [2.0,0.0]
1513        d = [0.0, 4.0]
1514        e = [2.0, 2.0]
1515        f = [4.0,0.0]
1516
1517        points = [a, b, c, d, e, f]
1518        #bac, bce, ecf, dbe
1519        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1520
1521        domain = Domain(points, vertices)
1522
1523        quantity = Quantity(domain)
1524        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1525       
1526        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
1527        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
1528
1529
1530        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
1531                                                                   [4,2,6],
1532                                                                   [6,2,4],
1533                                                                   [8,4,6]])
1534       
1535        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
1536                                                                                  [8,4,6]])
1537
1538
1539        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
1540                                                                [4,5,3],
1541                                                                [3,5,4],
1542                                                                [5,7,6]])
1543        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
1544                        [[4,5,3],
1545                         [5,7,6]])       
1546
1547        # Check averaging over vertices
1548        #a: 0
1549        #b: (4+4+4)/3
1550        #c: (2+2+2)/3
1551        #d: 8
1552        #e: (6+6+6)/3       
1553        #f: 4
1554        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])       
1555                                                                                 
1556       
1557
1558
1559
1560
1561    def test_get_interpolated_values(self):
1562
1563        from mesh_factory import rectangular
1564        from shallow_water import Domain
1565        from Numeric import zeros, Float
1566
1567        #Create basic mesh
1568        points, vertices, boundary = rectangular(1, 3)
1569        domain = Domain(points, vertices, boundary)
1570
1571        #Constant values
1572        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1573                                    [4,4,4],[5,5,5]])
1574
1575       
1576
1577        # Get interpolated values at centroids
1578        interpolation_points = domain.get_centroid_coordinates()
1579        answer = quantity.get_values(location='centroids')
1580
1581       
1582        #print quantity.get_values(points=interpolation_points)
1583        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
1584
1585
1586        #Arbitrary values
1587        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
1588                                    [1,4,-9],[2,5,0]])
1589
1590
1591        # Get interpolated values at centroids
1592        interpolation_points = domain.get_centroid_coordinates()
1593        answer = quantity.get_values(location='centroids')
1594        #print answer
1595        #print quantity.get_values(points=interpolation_points)
1596        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))       
1597                       
1598
1599        #FIXME TODO
1600        #indices = [0,5,3]
1601        #answer = [0.5,1,5]
1602        #assert allclose(answer,
1603        #                quantity.get_values(indices=indices, \
1604        #                                    location = 'unique vertices'))
1605
1606
1607
1608
1609    def test_get_interpolated_values_2(self):
1610        a = [0.0, 0.0]
1611        b = [0.0, 2.0]
1612        c = [2.0,0.0]
1613        d = [0.0, 4.0]
1614        e = [2.0, 2.0]
1615        f = [4.0,0.0]
1616
1617        points = [a, b, c, d, e, f]
1618        #bac, bce, ecf, dbe
1619        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1620
1621        domain = Domain(points, vertices)
1622
1623        quantity = Quantity(domain)
1624        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
1625
1626        #First pick one point
1627        x, y = 2.0/3, 8.0/3
1628        v = quantity.get_values(interpolation_points = [[x,y]])
1629        assert allclose(v, 6)       
1630
1631        # Then another to test that algorithm won't blindly
1632        # reuse interpolation matrix
1633        x, y = 4.0/3, 4.0/3
1634        v = quantity.get_values(interpolation_points = [[x,y]])
1635        assert allclose(v, 4)       
1636
1637
1638
1639
1640    def test_getting_some_vertex_values(self):
1641        """
1642        get values based on triangle lists.
1643        """
1644        from mesh_factory import rectangular
1645        from shallow_water import Domain
1646        from Numeric import zeros, Float
1647
1648        #Create basic mesh
1649        points, vertices, boundary = rectangular(1, 3)
1650
1651        #print "points",points
1652        #print "vertices",vertices
1653        #print "boundary",boundary
1654
1655        #Create shallow water domain
1656        domain = Domain(points, vertices, boundary)
1657        #print "domain.number_of_elements ",domain.number_of_elements
1658        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1659                                    [4,4,4],[5,5,5],[6,6,6]])
1660        value = [7]
1661        indices = [1]
1662        quantity.set_values(value,
1663                            location = 'centroids',
1664                            indices = indices)
1665        #print "quantity.centroid_values",quantity.centroid_values
1666        #print "quantity.get_values(location = 'centroids') ",\
1667        #      quantity.get_values(location = 'centroids')
1668        assert allclose(quantity.centroid_values,
1669                        quantity.get_values(location = 'centroids'))
1670
1671
1672        value = [[15,20,25]]
1673        quantity.set_values(value, indices = indices)
1674        #print "1 quantity.vertex_values",quantity.vertex_values
1675        assert allclose(quantity.vertex_values, quantity.get_values())
1676
1677        assert allclose(quantity.edge_values,
1678                        quantity.get_values(location = 'edges'))
1679
1680        # get a subset of elements
1681        subset = quantity.get_values(location='centroids', indices=[0,5])
1682        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1683        assert allclose(subset, answer)
1684
1685
1686        subset = quantity.get_values(location='edges', indices=[0,5])
1687        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1688        #print "subset",subset
1689        #print "answer",answer
1690        assert allclose(subset, answer)
1691
1692        subset = quantity.get_values( indices=[1,5])
1693        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1694        #print "subset",subset
1695        #print "answer",answer
1696        assert allclose(subset, answer)
1697
1698
1699
1700
1701#-------------------------------------------------------------
1702if __name__ == "__main__":
1703    suite = unittest.makeSuite(Test_Quantity, 'test')
1704    #print "restricted test"
1705    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_from_file_with_georef2')
1706    runner = unittest.TextTestRunner()
1707    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.