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

Last change on this file since 1632 was 1504, checked in by ole, 20 years ago

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