source: inundation-numpy-branch/pyvolution/test_quantity.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 47.5 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 numpy import allclose, array, zeros, ones, Float, sum
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 anuga.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 anuga.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 anuga.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 numpy 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        for k in range(quantity.centroid_values.shape[0]):
926            assert allclose (quantity.centroid_values[k],
927                             sum(quantity.vertex_values[k,:])/3)
928
929
930    def test_limiter2(self):
931        """Taken from test_shallow_water
932        """
933        quantity = Conserved_quantity(self.mesh4)
934
935        #Test centroids
936        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
937        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
938
939
940        #Extrapolate
941        quantity.extrapolate_second_order()
942
943        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
944
945        #Limit
946        quantity.limit()
947
948
949        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
950
951
952        #Assert that quantities are conserved
953        for k in range(quantity.centroid_values.shape[0]):
954            assert allclose (quantity.centroid_values[k],
955                             sum(quantity.vertex_values[k,:])/3)
956
957
958
959
960
961    def test_distribute_first_order(self):
962        quantity = Conserved_quantity(self.mesh4)
963
964        #Test centroids
965        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
966        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
967
968
969        #Extrapolate
970        quantity.extrapolate_first_order()
971
972        #Interpolate
973        quantity.interpolate_from_vertices_to_edges()
974
975        assert allclose(quantity.vertex_values,
976                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
977        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
978                                               [3,3,3], [4, 4, 4]])
979
980
981    def test_distribute_second_order(self):
982        quantity = Conserved_quantity(self.mesh4)
983
984        #Test centroids
985        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
986        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
987
988
989        #Extrapolate
990        quantity.extrapolate_second_order()
991
992        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
993
994
995    def test_update_explicit(self):
996        quantity = Conserved_quantity(self.mesh4)
997
998        #Test centroids
999        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1000        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1001
1002        #Set explicit_update
1003        quantity.explicit_update = array( [1.,1.,1.,1.] )
1004
1005        #Update with given timestep
1006        quantity.update(0.1)
1007
1008        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
1009        assert allclose( quantity.centroid_values, x)
1010
1011    def test_update_semi_implicit(self):
1012        quantity = Conserved_quantity(self.mesh4)
1013
1014        #Test centroids
1015        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1016        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1017
1018        #Set semi implicit update
1019        quantity.semi_implicit_update = array([1.,1.,1.,1.])
1020
1021        #Update with given timestep
1022        timestep = 0.1
1023        quantity.update(timestep)
1024
1025        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1026        denom = ones(4, Float)-timestep*sem
1027
1028        x = array([1, 2, 3, 4])/denom
1029        assert allclose( quantity.centroid_values, x)
1030
1031
1032    def test_both_updates(self):
1033        quantity = Conserved_quantity(self.mesh4)
1034
1035        #Test centroids
1036        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
1037        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
1038
1039        #Set explicit_update
1040        quantity.explicit_update = array( [4.,3.,2.,1.] )
1041
1042        #Set semi implicit update
1043        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
1044
1045        #Update with given timestep
1046        timestep = 0.1
1047        quantity.update(0.1)
1048
1049        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
1050        denom = ones(4, Float)-timestep*sem
1051
1052        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
1053        x /= denom
1054        assert allclose( quantity.centroid_values, x)
1055
1056
1057
1058
1059    #Test smoothing
1060    def test_smoothing(self):
1061
1062        from mesh_factory import rectangular
1063        from shallow_water import Domain, Transmissive_boundary
1064        from anuga.utilities.numerical_tools import mean
1065
1066        #Create basic mesh
1067        points, vertices, boundary = rectangular(2, 2)
1068
1069        #Create shallow water domain
1070        domain = Domain(points, vertices, boundary)
1071        domain.default_order=2
1072        domain.reduction = mean
1073
1074
1075        #Set some field values
1076        domain.set_quantity('elevation', lambda x,y: x)
1077        domain.set_quantity('friction', 0.03)
1078
1079
1080        ######################
1081        # Boundary conditions
1082        B = Transmissive_boundary(domain)
1083        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1084
1085
1086        ######################
1087        #Initial condition - with jumps
1088
1089        bed = domain.quantities['elevation'].vertex_values
1090        stage = zeros(bed.shape, Float)
1091
1092        h = 0.03
1093        for i in range(stage.shape[0]):
1094            if i % 2 == 0:
1095                stage[i,:] = bed[i,:] + h
1096            else:
1097                stage[i,:] = bed[i,:]
1098
1099        domain.set_quantity('stage', stage)
1100
1101        stage = domain.quantities['stage']
1102
1103        #Get smoothed stage
1104        A, V = stage.get_vertex_values(xy=False, smooth=True)
1105        Q = stage.vertex_values
1106
1107
1108        assert A.shape[0] == 9
1109        assert V.shape[0] == 8
1110        assert V.shape[1] == 3
1111
1112        #First four points
1113        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
1114        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
1115        assert allclose(A[2], Q[3,0])
1116        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
1117
1118        #Center point
1119        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
1120                               Q[5,0] + Q[6,2] + Q[7,1])/6)
1121
1122
1123        #Check V
1124        assert allclose(V[0,:], [3,4,0])
1125        assert allclose(V[1,:], [1,0,4])
1126        assert allclose(V[2,:], [4,5,1])
1127        assert allclose(V[3,:], [2,1,5])
1128        assert allclose(V[4,:], [6,7,3])
1129        assert allclose(V[5,:], [4,3,7])
1130        assert allclose(V[6,:], [7,8,4])
1131        assert allclose(V[7,:], [5,4,8])
1132
1133        #Get smoothed stage with XY
1134        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
1135
1136        assert allclose(A, A1)
1137        assert allclose(V, V1)
1138
1139        #Check XY
1140        assert allclose(X[4], 0.5)
1141        assert allclose(Y[4], 0.5)
1142
1143        assert allclose(X[7], 1.0)
1144        assert allclose(Y[7], 0.5)
1145
1146
1147
1148
1149    def test_vertex_values_no_smoothing(self):
1150
1151        from mesh_factory import rectangular
1152        from shallow_water import Domain, Transmissive_boundary
1153        from anuga.utilities.numerical_tools import mean
1154
1155
1156        #Create basic mesh
1157        points, vertices, boundary = rectangular(2, 2)
1158
1159        #Create shallow water domain
1160        domain = Domain(points, vertices, boundary)
1161        domain.default_order=2
1162        domain.reduction = mean
1163
1164
1165        #Set some field values
1166        domain.set_quantity('elevation', lambda x,y: x)
1167        domain.set_quantity('friction', 0.03)
1168
1169
1170        ######################
1171        #Initial condition - with jumps
1172
1173        bed = domain.quantities['elevation'].vertex_values
1174        stage = zeros(bed.shape, Float)
1175
1176        h = 0.03
1177        for i in range(stage.shape[0]):
1178            if i % 2 == 0:
1179                stage[i,:] = bed[i,:] + h
1180            else:
1181                stage[i,:] = bed[i,:]
1182
1183        domain.set_quantity('stage', stage)
1184
1185        #Get stage
1186        stage = domain.quantities['stage']
1187        A, V = stage.get_vertex_values(xy=False, smooth=False)
1188        Q = stage.vertex_values.flat
1189
1190        for k in range(8):
1191            assert allclose(A[k], Q[k])
1192
1193
1194        for k in range(8):
1195            assert V[k, 0] == 3*k
1196            assert V[k, 1] == 3*k+1
1197            assert V[k, 2] == 3*k+2
1198
1199
1200
1201        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
1202
1203
1204        assert allclose(A, A1)
1205        assert allclose(V, V1)
1206
1207        #Check XY
1208        assert allclose(X[1], 0.5)
1209        assert allclose(Y[1], 0.5)
1210        assert allclose(X[4], 0.0)
1211        assert allclose(Y[4], 0.0)
1212        assert allclose(X[12], 1.0)
1213        assert allclose(Y[12], 0.0)
1214
1215
1216
1217    def set_array_values_by_index(self):
1218
1219        from mesh_factory import rectangular
1220        from shallow_water import Domain
1221
1222        #Create basic mesh
1223        points, vertices, boundary = rectangular(1, 1)
1224
1225        #Create shallow water domain
1226        domain = Domain(points, vertices, boundary)
1227        #print "domain.number_of_elements ",domain.number_of_elements
1228        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1229        value = [7]
1230        indices = [1]
1231        quantity.set_array_values_by_index(value,
1232                                           location = 'centroids',
1233                                           indices = indices)
1234        #print "quantity.centroid_values",quantity.centroid_values
1235
1236        assert allclose(quantity.centroid_values, [1,7])
1237
1238        quantity.set_array_values([15,20,25], indices = indices)
1239        assert allclose(quantity.centroid_values, [1,20])
1240
1241        quantity.set_array_values([15,20,25], indices = indices)
1242        assert allclose(quantity.centroid_values, [1,20])
1243
1244    def test_setting_some_vertex_values(self):
1245        """
1246        set values based on triangle lists.
1247        """
1248        from mesh_factory import rectangular
1249        from shallow_water import Domain
1250
1251        #Create basic mesh
1252        points, vertices, boundary = rectangular(1, 3)
1253        #print "vertices",vertices
1254        #Create shallow water domain
1255        domain = Domain(points, vertices, boundary)
1256        #print "domain.number_of_elements ",domain.number_of_elements
1257        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1258                                    [4,4,4],[5,5,5],[6,6,6]])
1259        value = [7]
1260        indices = [1]
1261        quantity.set_values(value,
1262                            location = 'centroids',
1263                            indices = indices)
1264        #print "quantity.centroid_values",quantity.centroid_values
1265        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1266
1267        value = [[15,20,25]]
1268        quantity.set_values(value, indices = indices)
1269        #print "1 quantity.vertex_values",quantity.vertex_values
1270        assert allclose(quantity.vertex_values[1], value[0])
1271
1272
1273        #print "quantity",quantity.vertex_values
1274        values = [10,100,50]
1275        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1276        #print "2 quantity.vertex_values",quantity.vertex_values
1277        assert allclose(quantity.vertex_values[0], [10,10,10])
1278        assert allclose(quantity.vertex_values[5], [50,50,50])
1279        #quantity.interpolate()
1280        #print "quantity.centroid_values",quantity.centroid_values
1281        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1282
1283
1284        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1285                                    [4,4,4],[5,5,5],[6,6,6]])
1286        values = [10,100,50]
1287        #this will be per unique vertex, indexing the vertices
1288        #print "quantity.vertex_values",quantity.vertex_values
1289        quantity.set_values(values, indices = [0,1,5])
1290        #print "quantity.vertex_values",quantity.vertex_values
1291        assert allclose(quantity.vertex_values[0], [1,50,10])
1292        assert allclose(quantity.vertex_values[5], [6,6,6])
1293        assert allclose(quantity.vertex_values[1], [100,10,50])
1294
1295        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1296                                    [4,4,4],[5,5,5],[6,6,6]])
1297        values = [[31,30,29],[400,400,400],[1000,999,998]]
1298        quantity.set_values(values, indices = [3,3,5])
1299        quantity.interpolate()
1300        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1301
1302        values = [[1,1,1],[2,2,2],[3,3,3],
1303                                    [4,4,4],[5,5,5],[6,6,6]]
1304        quantity.set_values(values)
1305
1306        # testing the standard set values by vertex
1307        # indexed by vertex_id in general_mesh.coordinates
1308        values = [0,1,2,3,4,5,6,7]
1309
1310        quantity.set_values(values)
1311        #print "1 quantity.vertex_values",quantity.vertex_values
1312        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1313                                                [ 1.,  0.,  5.],
1314                                                [ 5.,  6.,  1.],
1315                                                [ 2.,  1.,  6.],
1316                                                [ 6.,  7.,  2.],
1317                                                [ 3.,  2.,  7.]])
1318
1319    def test_setting_unique_vertex_values(self):
1320        """
1321        set values based on unique_vertex lists.
1322        """
1323        from mesh_factory import rectangular
1324        from shallow_water import Domain
1325
1326        #Create basic mesh
1327        points, vertices, boundary = rectangular(1, 3)
1328        #print "vertices",vertices
1329        #Create shallow water domain
1330        domain = Domain(points, vertices, boundary)
1331        #print "domain.number_of_elements ",domain.number_of_elements
1332        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1333                                    [4,4,4],[5,5,5]])
1334        value = 7
1335        indices = [1,5]
1336        quantity.set_values(value,
1337                            location = 'unique vertices',
1338                            indices = indices)
1339        #print "quantity.centroid_values",quantity.centroid_values
1340        assert allclose(quantity.vertex_values[0], [0,7,0])
1341        assert allclose(quantity.vertex_values[1], [7,1,7])
1342        assert allclose(quantity.vertex_values[2], [7,2,7])
1343
1344
1345    def test_get_values(self):
1346        """
1347        get values based on triangle lists.
1348        """
1349        from mesh_factory import rectangular
1350        from shallow_water import Domain
1351
1352        #Create basic mesh
1353        points, vertices, boundary = rectangular(1, 3)
1354
1355        #print "points",points
1356        #print "vertices",vertices
1357        #print "boundary",boundary
1358
1359        #Create shallow water domain
1360        domain = Domain(points, vertices, boundary)
1361        #print "domain.number_of_elements ",domain.number_of_elements
1362        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1363                                    [4,4,4],[5,5,5]])
1364
1365        #print "quantity.get_values(location = 'unique vertices')", \
1366        #      quantity.get_values(location = 'unique vertices')
1367
1368        #print "quantity.get_values(location = 'unique vertices')", \
1369        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1370        #                          location = 'unique vertices')
1371
1372        answer = [0.5,2,4,5,0,1,3,4.5]
1373        assert allclose(answer,
1374                        quantity.get_values(location = 'unique vertices'))
1375
1376        indices = [0,5,3]
1377        answer = [0.5,1,5]
1378        assert allclose(answer,
1379                        quantity.get_values(indices=indices, \
1380                                            location = 'unique vertices'))
1381        #print "quantity.centroid_values",quantity.centroid_values
1382        #print "quantity.get_values(location = 'centroids') ",\
1383        #      quantity.get_values(location = 'centroids')
1384
1385    def test_getting_some_vertex_values(self):
1386        """
1387        get values based on triangle lists.
1388        """
1389        from mesh_factory import rectangular
1390        from shallow_water import Domain
1391
1392        #Create basic mesh
1393        points, vertices, boundary = rectangular(1, 3)
1394
1395        #print "points",points
1396        #print "vertices",vertices
1397        #print "boundary",boundary
1398
1399        #Create shallow water domain
1400        domain = Domain(points, vertices, boundary)
1401        #print "domain.number_of_elements ",domain.number_of_elements
1402        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1403                                    [4,4,4],[5,5,5],[6,6,6]])
1404        value = [7]
1405        indices = [1]
1406        quantity.set_values(value,
1407                            location = 'centroids',
1408                            indices = indices)
1409        #print "quantity.centroid_values",quantity.centroid_values
1410        #print "quantity.get_values(location = 'centroids') ",\
1411        #      quantity.get_values(location = 'centroids')
1412        assert allclose(quantity.centroid_values,
1413                        quantity.get_values(location = 'centroids'))
1414
1415
1416        value = [[15,20,25]]
1417        quantity.set_values(value, indices = indices)
1418        #print "1 quantity.vertex_values",quantity.vertex_values
1419        assert allclose(quantity.vertex_values, quantity.get_values())
1420
1421        assert allclose(quantity.edge_values,
1422                        quantity.get_values(location = 'edges'))
1423
1424        # get a subset of elements
1425        subset = quantity.get_values(location='centroids', indices=[0,5])
1426        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1427        assert allclose(subset, answer)
1428
1429
1430        subset = quantity.get_values(location='edges', indices=[0,5])
1431        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1432        #print "subset",subset
1433        #print "answer",answer
1434        assert allclose(subset, answer)
1435
1436        subset = quantity.get_values( indices=[1,5])
1437        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1438        #print "subset",subset
1439        #print "answer",answer
1440        assert allclose(subset, answer)
1441
1442
1443
1444
1445#-------------------------------------------------------------
1446if __name__ == "__main__":
1447    suite = unittest.makeSuite(Test_Quantity, 'test')
1448    #print "restricted test"
1449    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1450    runner = unittest.TextTestRunner()
1451    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.