source: inundation/pyvolution/test_quantity.py @ 2347

Last change on this file since 2347 was 2347, checked in by ole, 19 years ago

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