source: inundation/pyvolution/test_quantity.py @ 2891

Last change on this file since 2891 was 2754, checked in by duncan, 19 years ago

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