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

Last change on this file since 3678 was 3678, checked in by steve, 18 years ago

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