source: inundation/ga/storm_surge/pyvolution/test_quantity.py @ 1108

Last change on this file since 1108 was 1018, checked in by steve, 20 years ago

Cleaned up test function
File size: 30.9 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
12class Test_Quantity(unittest.TestCase):
13    def setUp(self):
14        from domain import Domain
15
16        a = [0.0, 0.0]
17        b = [0.0, 2.0]
18        c = [2.0, 0.0]
19        d = [0.0, 4.0]
20        e = [2.0, 2.0]
21        f = [4.0, 0.0]
22
23        points = [a, b, c, d, e, f]
24
25        #bac, bce, ecf, dbe
26        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
27
28        self.mesh1 = Domain(points[:3], [elements[0]])
29        self.mesh1.check_integrity()
30
31        self.mesh4 = Domain(points, elements)
32        self.mesh4.check_integrity()
33
34    def tearDown(self):
35        pass
36        #print "  Tearing down"
37
38
39    def test_creation(self):
40
41        quantity = Quantity(self.mesh1, [[1,2,3]])
42        assert allclose(quantity.vertex_values, [[1.,2.,3.]])
43
44        try:
45            quantity = Quantity()
46        except:
47            pass
48        else:
49            raise 'Should have raised empty quantity exception'
50
51
52        try:
53            quantity = Quantity([1,2,3])
54        except AssertionError:
55            pass
56        except:
57            raise 'Should have raised "mising mesh object" error'
58
59
60    def test_creation_zeros(self):
61
62        quantity = Quantity(self.mesh1)
63        assert allclose(quantity.vertex_values, [[0.,0.,0.]])
64
65
66        quantity = Quantity(self.mesh4)
67        assert allclose(quantity.vertex_values, [[0.,0.,0.], [0.,0.,0.],
68                                                 [0.,0.,0.], [0.,0.,0.]])
69
70
71    def test_interpolation(self):
72        quantity = Quantity(self.mesh1, [[1,2,3]])
73        assert allclose(quantity.centroid_values, [2.0]) #Centroid
74
75        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5]])
76
77
78    def test_interpolation2(self):
79        quantity = Conserved_quantity(self.mesh4,
80                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
81        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
82
83
84        quantity.extrapolate_second_order()
85
86        assert allclose(quantity.vertex_values, [[2., 2., 2.],
87                                               [3.+2./3, 6.+2./3, 4.+2./3],
88                                               [7.5, 0.5, 1.],
89                                               [-5, -2.5, 7.5]])
90
91        #print quantity.edge_values
92        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
93                                               [5., 5., 5.],
94                                               [4.5, 4.5, 0.],
95                                               [3.0, -1.5, -1.5]])
96
97    def test_boundary_allocation(self):
98        quantity = Conserved_quantity(self.mesh4,
99                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
100
101        assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)
102
103
104    def test_set_values(self):
105        quantity = Quantity(self.mesh4)
106
107
108        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
109                            location = 'vertices')
110        assert allclose(quantity.vertex_values,
111                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
112        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
113        assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
114                                               [5., 5., 5.],
115                                               [4.5, 4.5, 0.],
116                                               [3.0, -1.5, -1.5]])
117
118
119        #Test default
120        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
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        #Test centroids
130        quantity.set_values([1,2,3,4], location = 'centroids')
131        assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
132
133        #Test edges
134        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
135                            location = 'edges')
136        assert allclose(quantity.edge_values,
137                        [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
138
139        #Test exceptions
140        try:
141            quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
142                                location = 'bas kamel tuba')
143        except:
144            pass
145
146
147        try:
148            quantity.set_values([[1,2,3], [0,0,9]])
149        except AssertionError:
150            pass
151        except:
152            raise 'should have raised Assertionerror'
153
154
155
156    def test_set_values_const(self):
157        quantity = Quantity(self.mesh4)
158
159        quantity.set_values(1.0, location = 'vertices')
160        assert allclose(quantity.vertex_values,
161                        [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
162        assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
163        assert allclose(quantity.edge_values, [[1, 1, 1],
164                                               [1, 1, 1],
165                                               [1, 1, 1],
166                                               [1, 1, 1]])
167
168
169        quantity.set_values(2.0, location = 'centroids')
170        assert allclose(quantity.centroid_values, [2, 2, 2, 2])
171
172        quantity.set_values(3.0, location = 'edges')
173        assert allclose(quantity.edge_values, [[3, 3, 3],
174                                               [3, 3, 3],
175                                               [3, 3, 3],
176                                               [3, 3, 3]])
177
178
179    def test_set_values_func(self):
180        quantity = Quantity(self.mesh4)
181
182        def f(x, y):
183            return x+y
184
185        quantity.set_values(f, location = 'vertices')
186        #print "quantity.vertex_values",quantity.vertex_values
187        assert allclose(quantity.vertex_values,
188                        [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])
189        assert allclose(quantity.centroid_values,
190                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
191        assert allclose(quantity.edge_values,
192                        [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])
193
194
195        quantity.set_values(f, location = 'centroids')
196        assert allclose(quantity.centroid_values,
197                        [4.0/3, 8.0/3, 10.0/3, 10.0/3])
198
199
200    def test_integral(self):
201        quantity = Quantity(self.mesh4)
202
203        #Try constants first
204        const = 5
205        quantity.set_values(const, location = 'vertices')
206        #print 'Q', quantity.get_integral()
207
208        assert allclose(quantity.get_integral(),
209                        self.mesh4.get_area() * const)
210
211        #Try with a linear function
212        def f(x, y):
213            return x+y
214
215        quantity.set_values(f, location = 'vertices')
216
217
218        ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2
219
220        assert allclose (quantity.get_integral(), ref_integral)
221
222
223
224    def test_set_vertex_values(self):
225        quantity = Quantity(self.mesh4)
226        quantity.set_vertex_values([0,1,2,3,4,5])
227
228        assert allclose(quantity.vertex_values,
229                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
230        assert allclose(quantity.centroid_values,
231                        [1., 7./3, 11./3, 8./3]) #Centroid
232        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
233                                               [3., 2.5, 1.5],
234                                               [3.5, 4.5, 3.],
235                                               [2.5, 3.5, 2]])
236
237
238    def test_set_vertex_values_subset(self):
239        quantity = Quantity(self.mesh4)
240        quantity.set_vertex_values([0,1,2,3,4,5])
241        quantity.set_vertex_values([0,20,30,50], indexes = [0,2,3,5])
242
243        assert allclose(quantity.vertex_values,
244                        [[1,0,20], [1,20,4], [4,20,50], [30,1,4]])
245
246    def test_set_vertex_values_using_general_interface(self):
247        quantity = Quantity(self.mesh4)
248
249
250        quantity.set_values([0,1,2,3,4,5])
251
252
253        assert allclose(quantity.vertex_values,
254                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
255
256        #Centroid
257        assert allclose(quantity.centroid_values, [1., 7./3, 11./3, 8./3])
258
259        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
260                                               [3., 2.5, 1.5],
261                                               [3.5, 4.5, 3.],
262                                               [2.5, 3.5, 2]])
263
264
265
266
267    def test_gradient(self):
268        quantity = Conserved_quantity(self.mesh4)
269
270        #Set up for a gradient of (3,0) at mid triangle
271        quantity.set_values([2.0, 4.0, 8.0, 2.0],
272                            location = 'centroids')
273
274
275
276        #Gradients
277        a, b = quantity.compute_gradients()
278
279        #gradient bewteen t0 and t1 is undefined as det==0
280        assert a[0] == 0.0
281        assert b[0] == 0.0
282        #The others are OK
283        for i in range(1,4):
284            assert a[i] == 3.0
285            assert b[i] == 0.0
286
287
288        quantity.extrapolate_second_order()
289
290        #print quantity.vertex_values
291        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
292                                                 [0., 6.,  6.],
293                                                 [6., 6., 12.],
294                                                 [0., 0.,  6.]])
295
296
297    def test_second_order_extrapolation2(self):
298        quantity = Conserved_quantity(self.mesh4)
299
300        #Set up for a gradient of (3,1), f(x) = 3x+y
301        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
302                            location = 'centroids')
303
304        #Gradients
305        a, b = quantity.compute_gradients()
306
307        #gradient bewteen t0 and t1 is undefined as det==0
308        assert a[0] == 0.0
309        assert b[0] == 0.0
310        #The others are OK
311        for i in range(1,4):
312            assert allclose(a[i], 3.0)
313            assert allclose(b[i], 1.0)
314
315        quantity.extrapolate_second_order()
316
317        #print quantity.vertex_values
318        assert allclose(quantity.vertex_values[1,0], 2.0)
319        assert allclose(quantity.vertex_values[1,1], 6.0)
320        assert allclose(quantity.vertex_values[1,2], 8.0)
321
322
323
324    def test_first_order_extrapolator(self):
325        quantity = Conserved_quantity(self.mesh4)
326
327        #Test centroids
328        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
329        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
330
331        #Extrapolate
332        quantity.extrapolate_first_order()
333
334        #Check vertices but not edge values
335        assert allclose(quantity.vertex_values,
336                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
337
338
339    def test_second_order_extrapolator(self):
340        quantity = Conserved_quantity(self.mesh4)
341
342        #Set up for a gradient of (3,0) at mid triangle
343        quantity.set_values([2.0, 4.0, 8.0, 2.0],
344                            location = 'centroids')
345
346
347
348        quantity.extrapolate_second_order()
349        quantity.limit()
350
351
352        #Assert that central triangle is limited by neighbours
353        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
354        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
355
356        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
357        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
358
359        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
360        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
361
362
363        #Assert that quantities are conserved
364        from Numeric import sum
365        for k in range(quantity.centroid_values.shape[0]):
366            assert allclose (quantity.centroid_values[k],
367                             sum(quantity.vertex_values[k,:])/3)
368
369
370
371
372
373    def test_limiter(self):
374        quantity = Conserved_quantity(self.mesh4)
375
376        #Create a deliberate overshoot (e.g. from gradient computation)
377        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
378
379
380        #Limit
381        quantity.limit()
382
383        #Assert that central triangle is limited by neighbours
384        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
385        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
386
387        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
388        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
389
390        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
391        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
392
393
394
395        #Assert that quantities are conserved
396        from Numeric import sum
397        for k in range(quantity.centroid_values.shape[0]):
398            assert allclose (quantity.centroid_values[k],
399                             sum(quantity.vertex_values[k,:])/3)
400
401
402    def test_limiter2(self):
403        """Taken from test_shallow_water
404        """
405        quantity = Conserved_quantity(self.mesh4)
406
407        #Test centroids
408        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
409        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
410
411
412        #Extrapolate
413        quantity.extrapolate_second_order()
414
415        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
416
417        #Limit
418        quantity.limit()
419
420
421        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
422
423
424        #Assert that quantities are conserved
425        from Numeric import sum
426        for k in range(quantity.centroid_values.shape[0]):
427            assert allclose (quantity.centroid_values[k],
428                             sum(quantity.vertex_values[k,:])/3)
429
430
431
432
433
434    def test_distribute_first_order(self):
435        quantity = Conserved_quantity(self.mesh4)
436
437        #Test centroids
438        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
439        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
440
441
442        #Extrapolate
443        quantity.extrapolate_first_order()
444
445        #Interpolate
446        quantity.interpolate_from_vertices_to_edges()
447
448        assert allclose(quantity.vertex_values,
449                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
450        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
451                                               [3,3,3], [4, 4, 4]])
452
453
454    def test_distribute_second_order(self):
455        quantity = Conserved_quantity(self.mesh4)
456
457        #Test centroids
458        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
459        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
460
461
462        #Extrapolate
463        quantity.extrapolate_second_order()
464
465        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
466
467
468    def test_update_explicit(self):
469        quantity = Conserved_quantity(self.mesh4)
470
471        #Test centroids
472        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
473        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
474
475        #Set explicit_update
476        quantity.explicit_update = array( [1.,1.,1.,1.] )
477
478        #Update with given timestep
479        quantity.update(0.1)
480
481        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
482        assert allclose( quantity.centroid_values, x)
483
484    def test_update_semi_implicit(self):
485        quantity = Conserved_quantity(self.mesh4)
486
487        #Test centroids
488        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
489        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
490
491        #Set semi implicit update
492        quantity.semi_implicit_update = array([1.,1.,1.,1.])
493
494        #Update with given timestep
495        timestep = 0.1
496        quantity.update(timestep)
497
498        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
499        denom = ones(4, Float)-timestep*sem
500
501        x = array([1, 2, 3, 4])/denom
502        assert allclose( quantity.centroid_values, x)
503
504
505    def test_both_updates(self):
506        quantity = Conserved_quantity(self.mesh4)
507
508        #Test centroids
509        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
510        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
511
512        #Set explicit_update
513        quantity.explicit_update = array( [4.,3.,2.,1.] )
514
515        #Set semi implicit update
516        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
517
518        #Update with given timestep
519        timestep = 0.1
520        quantity.update(0.1)
521
522        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
523        denom = ones(4, Float)-timestep*sem
524
525        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
526        x /= denom
527        assert allclose( quantity.centroid_values, x)
528
529
530
531
532    #Test smoothing
533    def test_smoothing(self):
534
535        from mesh_factory import rectangular
536        from shallow_water import Domain, Transmissive_boundary
537        from Numeric import zeros, Float
538        from util import mean
539
540        #Create basic mesh
541        points, vertices, boundary = rectangular(2, 2)
542
543        #Create shallow water domain
544        domain = Domain(points, vertices, boundary)
545        domain.default_order=2
546        domain.reduction = mean
547
548
549        #Set some field values
550        domain.set_quantity('elevation', lambda x,y: x)
551        domain.set_quantity('friction', 0.03)
552
553
554        ######################
555        # Boundary conditions
556        B = Transmissive_boundary(domain)
557        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
558
559
560        ######################
561        #Initial condition - with jumps
562
563        bed = domain.quantities['elevation'].vertex_values
564        stage = zeros(bed.shape, Float)
565
566        h = 0.03
567        for i in range(stage.shape[0]):
568            if i % 2 == 0:
569                stage[i,:] = bed[i,:] + h
570            else:
571                stage[i,:] = bed[i,:]
572
573        domain.set_quantity('stage', stage)
574
575        stage = domain.quantities['stage']
576
577        #Get smoothed stage
578        A, V = stage.get_vertex_values(xy=False, smooth=True)
579        Q = stage.vertex_values
580
581
582        assert A.shape[0] == 9
583        assert V.shape[0] == 8
584        assert V.shape[1] == 3
585
586        #First four points
587        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
588        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
589        assert allclose(A[2], Q[3,0])
590        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
591
592        #Center point
593        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
594                               Q[5,0] + Q[6,2] + Q[7,1])/6)
595
596
597        #Check V
598        assert allclose(V[0,:], [3,4,0])
599        assert allclose(V[1,:], [1,0,4])
600        assert allclose(V[2,:], [4,5,1])
601        assert allclose(V[3,:], [2,1,5])
602        assert allclose(V[4,:], [6,7,3])
603        assert allclose(V[5,:], [4,3,7])
604        assert allclose(V[6,:], [7,8,4])
605        assert allclose(V[7,:], [5,4,8])
606
607        #Get smoothed stage with XY
608        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
609
610        assert allclose(A, A1)
611        assert allclose(V, V1)
612
613        #Check XY
614        assert allclose(X[4], 0.5)
615        assert allclose(Y[4], 0.5)
616
617        assert allclose(X[7], 1.0)
618        assert allclose(Y[7], 0.5)
619
620
621
622
623    def test_vertex_values_no_smoothing(self):
624
625        from mesh_factory import rectangular
626        from shallow_water import Domain, Transmissive_boundary
627        from Numeric import zeros, Float
628        from util import mean
629
630
631        #Create basic mesh
632        points, vertices, boundary = rectangular(2, 2)
633
634        #Create shallow water domain
635        domain = Domain(points, vertices, boundary)
636        domain.default_order=2
637        domain.reduction = mean
638
639
640        #Set some field values
641        domain.set_quantity('elevation', lambda x,y: x)
642        domain.set_quantity('friction', 0.03)
643
644
645        ######################
646        #Initial condition - with jumps
647
648        bed = domain.quantities['elevation'].vertex_values
649        stage = zeros(bed.shape, Float)
650
651        h = 0.03
652        for i in range(stage.shape[0]):
653            if i % 2 == 0:
654                stage[i,:] = bed[i,:] + h
655            else:
656                stage[i,:] = bed[i,:]
657
658        domain.set_quantity('stage', stage)
659
660        #Get stage
661        stage = domain.quantities['stage']
662        A, V = stage.get_vertex_values(xy=False, smooth=False)
663        Q = stage.vertex_values.flat
664
665        for k in range(8):
666            assert allclose(A[k], Q[k])
667
668
669        for k in range(8):
670            assert V[k, 0] == 3*k
671            assert V[k, 1] == 3*k+1
672            assert V[k, 2] == 3*k+2
673
674
675
676        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
677
678
679        assert allclose(A, A1)
680        assert allclose(V, V1)
681
682        #Check XY
683        assert allclose(X[1], 0.5)
684        assert allclose(Y[1], 0.5)
685        assert allclose(X[4], 0.0)
686        assert allclose(Y[4], 0.0)
687        assert allclose(X[12], 1.0)
688        assert allclose(Y[12], 0.0)
689
690
691
692    def set_array_values_by_index(self):
693
694        from mesh_factory import rectangular
695        from shallow_water import Domain
696        from Numeric import zeros, Float
697
698        #Create basic mesh
699        points, vertices, boundary = rectangular(1, 1)
700
701        #Create shallow water domain
702        domain = Domain(points, vertices, boundary)
703        #print "domain.number_of_elements ",domain.number_of_elements
704        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
705        value = [7]
706        indexes = [1]
707        quantity.set_array_values_by_index(value,
708                                           location = 'centroids',
709                                           indexes = indexes)
710        #print "quantity.centroid_values",quantity.centroid_values
711
712        assert allclose(quantity.centroid_values, [1,7])
713
714        quantity.set_array_values([15,20,25], indexes = indexes)
715        assert allclose(quantity.centroid_values, [1,20])
716
717        quantity.set_array_values([15,20,25], indexes = indexes)
718        assert allclose(quantity.centroid_values, [1,20])
719
720    def test_setting_some_vertex_values(self):
721        """
722        set values based on triangle lists.
723        """
724        from mesh_factory import rectangular
725        from shallow_water import Domain
726        from Numeric import zeros, Float
727
728        #Create basic mesh
729        points, vertices, boundary = rectangular(1, 3)
730        #print "vertices",vertices
731        #Create shallow water domain
732        domain = Domain(points, vertices, boundary)
733        #print "domain.number_of_elements ",domain.number_of_elements
734        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
735                                    [4,4,4],[5,5,5],[6,6,6]])
736        value = [7]
737        indexes = [1]
738        quantity.set_values(value,
739                                  location = 'centroids',
740                                  indexes = indexes)
741        #print "quantity.centroid_values",quantity.centroid_values
742        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
743
744        value = [[15,20,25]]
745        quantity.set_values(value, indexes = indexes)
746        #print "1 quantity.vertex_values",quantity.vertex_values
747        assert allclose(quantity.vertex_values[1], value[0])
748
749
750        #print "quantity",quantity.vertex_values
751        values = [10,100,50]
752        quantity.set_values(values, indexes = [0,1,5], location = 'centroids')
753        #print "2 quantity.vertex_values",quantity.vertex_values
754        assert allclose(quantity.vertex_values[0], [10,10,10])
755        assert allclose(quantity.vertex_values[5], [50,50,50])
756        #quantity.interpolate()
757        #print "quantity.centroid_values",quantity.centroid_values
758        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
759
760
761        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
762                                    [4,4,4],[5,5,5],[6,6,6]])
763        values = [10,100,50]
764        #this will be per unique vertex, indexing the vertices
765        #print "quantity.vertex_values",quantity.vertex_values
766        quantity.set_values(values, indexes = [0,1,5])
767        #print "quantity.vertex_values",quantity.vertex_values
768        assert allclose(quantity.vertex_values[0], [1,50,10])
769        assert allclose(quantity.vertex_values[5], [6,6,6])
770        assert allclose(quantity.vertex_values[1], [100,10,50])
771
772        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
773                                    [4,4,4],[5,5,5],[6,6,6]])
774        values = [[31,30,29],[400,400,400],[1000,999,998]]
775        quantity.set_values(values, indexes = [3,3,5])
776        quantity.interpolate()
777        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
778
779        values = [[1,1,1],[2,2,2],[3,3,3],
780                                    [4,4,4],[5,5,5],[6,6,6]]
781        quantity.set_values(values)
782
783        # testing the standard set values by vertex
784        # indexed by vertex_id in general_mesh.coordinates
785        values = [0,1,2,3,4,5,6,7]
786
787        quantity.set_values(values)
788        #print "1 quantity.vertex_values",quantity.vertex_values
789        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
790                                                [ 1.,  0.,  5.],
791                                                [ 5.,  6.,  1.],
792                                                [ 2.,  1.,  6.],
793                                                [ 6.,  7.,  2.],
794                                                [ 3.,  2.,  7.]])
795
796    def test_setting_unique_vertex_values(self):
797        """
798        set values based on unique_vertex lists.
799        """
800        from mesh_factory import rectangular
801        from shallow_water import Domain
802        from Numeric import zeros, Float
803
804        #Create basic mesh
805        points, vertices, boundary = rectangular(1, 3)
806        #print "vertices",vertices
807        #Create shallow water domain
808        domain = Domain(points, vertices, boundary)
809        #print "domain.number_of_elements ",domain.number_of_elements
810        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
811                                    [4,4,4],[5,5,5]])
812        value = 7
813        indexes = [1,5]
814        quantity.set_values(value,
815                                  location = 'unique vertices',
816                                  indexes = indexes)
817        #print "quantity.centroid_values",quantity.centroid_values
818        assert allclose(quantity.vertex_values[0], [0,7,0])
819        assert allclose(quantity.vertex_values[1], [7,1,7])
820        assert allclose(quantity.vertex_values[2], [7,2,7])
821
822
823    def test_get_values(self):
824        """
825        get values based on triangle lists.
826        """
827        from mesh_factory import rectangular
828        from shallow_water import Domain
829        from Numeric import zeros, Float
830
831        #Create basic mesh
832        points, vertices, boundary = rectangular(1, 3)
833
834        #print "points",points
835        #print "vertices",vertices
836        #print "boundary",boundary
837
838        #Create shallow water domain
839        domain = Domain(points, vertices, boundary)
840        #print "domain.number_of_elements ",domain.number_of_elements
841        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
842                                    [4,4,4],[5,5,5]])
843
844        #print "quantity.get_values(location = 'unique vertices')", \
845        #      quantity.get_values(location = 'unique vertices')
846
847        #print "quantity.get_values(location = 'unique vertices')", \
848        #      quantity.get_values(indexes=[0,1,2,3,4,5,6,7], \
849        #                          location = 'unique vertices')
850
851        answer = [0.5,2,4,5,0,1,3,4.5]
852        assert allclose(answer,
853                        quantity.get_values(location = 'unique vertices'))
854
855        indexes = [0,5,3]
856        answer = [0.5,1,5]
857        assert allclose(answer,
858                        quantity.get_values(indexes=indexes, \
859                                            location = 'unique vertices'))
860        #print "quantity.centroid_values",quantity.centroid_values
861        #print "quantity.get_values(location = 'centroids') ",\
862        #      quantity.get_values(location = 'centroids')
863
864    def test_getting_some_vertex_values(self):
865        """
866        get values based on triangle lists.
867        """
868        from mesh_factory import rectangular
869        from shallow_water import Domain
870        from Numeric import zeros, Float
871
872        #Create basic mesh
873        points, vertices, boundary = rectangular(1, 3)
874
875        #print "points",points
876        #print "vertices",vertices
877        #print "boundary",boundary
878
879        #Create shallow water domain
880        domain = Domain(points, vertices, boundary)
881        #print "domain.number_of_elements ",domain.number_of_elements
882        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
883                                    [4,4,4],[5,5,5],[6,6,6]])
884        value = [7]
885        indexes = [1]
886        quantity.set_values(value,
887                                  location = 'centroids',
888                                  indexes = indexes)
889        #print "quantity.centroid_values",quantity.centroid_values
890        #print "quantity.get_values(location = 'centroids') ",\
891        #      quantity.get_values(location = 'centroids')
892        assert allclose(quantity.centroid_values,
893                        quantity.get_values(location = 'centroids'))
894
895
896        value = [[15,20,25]]
897        quantity.set_values(value, indexes = indexes)
898        #print "1 quantity.vertex_values",quantity.vertex_values
899        assert allclose(quantity.vertex_values, quantity.get_values())
900
901        assert allclose(quantity.edge_values,
902                        quantity.get_values(location = 'edges'))
903
904        # get a subset of elements
905        subset = quantity.get_values(location='centroids', indexes=[0,5])
906        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
907        assert allclose(subset, answer)
908
909
910        subset = quantity.get_values(location='edges', indexes=[0,5])
911        answer = [quantity.edge_values[0],quantity.edge_values[5]]
912        #print "subset",subset
913        #print "answer",answer
914        assert allclose(subset, answer)
915
916        subset = quantity.get_values( indexes=[1,5])
917        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
918        #print "subset",subset
919        #print "answer",answer
920        assert allclose(subset, answer)
921
922
923
924#-------------------------------------------------------------
925if __name__ == "__main__":
926    suite = unittest.makeSuite(Test_Quantity,'test')
927    #print "restricted test"
928    #suite = unittest.makeSuite(Test_Quantity,'test_set_vertex_values_subset')
929    runner = unittest.TextTestRunner()
930    runner.run(suite)
931
Note: See TracBrowser for help on using the repository browser.