source: inundation/pyvolution/test_quantity.py @ 2267

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

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