source: inundation/pyvolution/test_quantity.py @ 1753

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

Embedded caching functionality within quantity.set_values and modified validation example lwru2.py to illustrate the advantages that can be gained from supervised caching.

File size: 35.6 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_vertex_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        #Obsoleted bit
303        #interp = Interpolation(quantity.domain.coordinates,
304        #                       quantity.domain.triangles,
305        #                       data_points, alpha=0.0)
306        #
307        #answer = linear_function(quantity.domain.coordinates)
308        #
309        #f = interp.fit(z)
310        #
311        #print "f",f
312        #print "answer",answer
313        #assert allclose(f, answer)
314
315
316        #Use built-in least squares fit
317        quantity.set_values(points = data_points, values = z, alpha = 0)
318        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
319        #print quantity.vertex_values, answer
320        assert allclose(quantity.vertex_values.flat, answer)
321
322
323        #Now try by setting the same values directly
324        from least_squares import fit_to_mesh       
325        vertex_attributes = fit_to_mesh(quantity.domain.coordinates,
326                                        quantity.domain.triangles,
327                                        data_points,
328                                        z,
329                                        alpha = 0,
330                                        verbose=False)
331
332        #print vertex_attributes
333        quantity.set_values(vertex_attributes)
334        assert allclose(quantity.vertex_values.flat, answer)
335
336
337    def test_set_values_from_file(self):
338        quantity = Quantity(self.mesh4)
339
340        #Get (enough) datapoints
341        data_points = [[ 0.66666667, 0.66666667],
342                       [ 1.33333333, 1.33333333],
343                       [ 2.66666667, 0.66666667],
344                       [ 0.66666667, 2.66666667],
345                       [ 0.0, 1.0],
346                       [ 0.0, 3.0],
347                       [ 1.0, 0.0],
348                       [ 1.0, 1.0],
349                       [ 1.0, 2.0],
350                       [ 1.0, 3.0],
351                       [ 2.0, 1.0],
352                       [ 3.0, 0.0],
353                       [ 3.0, 1.0]]
354
355        z = linear_function(data_points)
356       
357
358        #Create pts file
359        from data_manager import write_ptsfile
360        ptsfile = 'testptsfile.pts'
361        att = 'spam_and_eggs'
362        write_ptsfile(ptsfile, data_points, z, attribute_name = att)
363
364        #Check that values can be set from file
365        quantity.set_values(filename = ptsfile,
366                            attribute_name = att, alpha = 0)
367        answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True))
368        #print quantity.vertex_values, answer
369        assert allclose(quantity.vertex_values.flat, answer)
370
371
372        #Check that values can be set from file using default attribute
373        quantity.set_values(filename = ptsfile, alpha = 0)
374        assert allclose(quantity.vertex_values.flat, answer)       
375
376        #Cleanup
377        import os
378        os.remove(ptsfile)
379       
380
381
382    def test_gradient(self):
383        quantity = Conserved_quantity(self.mesh4)
384
385        #Set up for a gradient of (2,0) at mid triangle
386        quantity.set_values([2.0, 4.0, 6.0, 2.0],
387                            location = 'centroids')
388
389
390        #Gradients
391        a, b = quantity.compute_gradients()
392
393        #print self.mesh4.centroid_coordinates
394        #print a, b
395
396        #The central triangle (1)
397        #(using standard gradient based on neigbours controid values)
398        assert allclose(a[1], 2.0)
399        assert allclose(b[1], 0.0)
400
401
402        #Left triangle (0) using two point gradient
403        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
404        #2  = 4  + a*(-2/3)  + b*(-2/3)
405        assert allclose(a[0] + b[0], 3)
406        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)       
407        assert allclose(a[0] - b[0], 0) 
408
409
410        #Right triangle (2) using two point gradient
411        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
412        #6  = 4  + a*(4/3)  + b*(-2/3)       
413        assert allclose(2*a[2] - b[2], 3)
414        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
415        assert allclose(a[2] + 2*b[2], 0) 
416
417
418        #Top triangle (3) using two point gradient
419        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
420        #2  = 4  + a*(-2/3)  + b*(4/3)       
421        assert allclose(a[3] - 2*b[3], 3)
422        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)       
423        assert allclose(2*a[3] + b[3], 0) 
424
425
426
427        #print a, b
428        quantity.extrapolate_second_order()
429
430        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
431        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
432        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
433
434
435        #a = 1.2, b=-0.6
436        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
437        assert allclose(quantity.vertex_values[2,2], 8)       
438
439
440    def test_second_order_extrapolation2(self):
441        quantity = Conserved_quantity(self.mesh4)
442
443        #Set up for a gradient of (3,1), f(x) = 3x+y
444        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
445                            location = 'centroids')
446
447        #Gradients
448        a, b = quantity.compute_gradients()
449
450        #print a, b
451
452        assert allclose(a[1], 3.0)
453        assert allclose(b[1], 1.0)
454       
455        #Work out the others
456       
457        quantity.extrapolate_second_order()
458
459        #print quantity.vertex_values
460        assert allclose(quantity.vertex_values[1,0], 2.0)
461        assert allclose(quantity.vertex_values[1,1], 6.0)
462        assert allclose(quantity.vertex_values[1,2], 8.0)
463
464
465
466    def test_first_order_extrapolator(self):
467        quantity = Conserved_quantity(self.mesh4)
468
469        #Test centroids
470        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
471        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
472
473        #Extrapolate
474        quantity.extrapolate_first_order()
475
476        #Check vertices but not edge values
477        assert allclose(quantity.vertex_values,
478                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
479
480
481    def test_second_order_extrapolator(self):
482        quantity = Conserved_quantity(self.mesh4)
483
484        #Set up for a gradient of (3,0) at mid triangle
485        quantity.set_values([2.0, 4.0, 8.0, 2.0],
486                            location = 'centroids')
487
488
489
490        quantity.extrapolate_second_order()
491        quantity.limit()
492
493
494        #Assert that central triangle is limited by neighbours
495        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
496        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
497
498        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
499        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
500
501        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
502        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
503
504
505        #Assert that quantities are conserved
506        from Numeric import sum
507        for k in range(quantity.centroid_values.shape[0]):
508            assert allclose (quantity.centroid_values[k],
509                             sum(quantity.vertex_values[k,:])/3)
510
511
512
513
514
515    def test_limiter(self):
516        quantity = Conserved_quantity(self.mesh4)
517
518        #Create a deliberate overshoot (e.g. from gradient computation)
519        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
520
521
522        #Limit
523        quantity.limit()
524
525        #Assert that central triangle is limited by neighbours
526        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
527        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
528
529        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
530        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
531
532        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
533        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
534
535
536
537        #Assert that quantities are conserved
538        from Numeric import sum
539        for k in range(quantity.centroid_values.shape[0]):
540            assert allclose (quantity.centroid_values[k],
541                             sum(quantity.vertex_values[k,:])/3)
542
543
544    def test_limiter2(self):
545        """Taken from test_shallow_water
546        """
547        quantity = Conserved_quantity(self.mesh4)
548
549        #Test centroids
550        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
551        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
552
553
554        #Extrapolate
555        quantity.extrapolate_second_order()
556
557        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
558
559        #Limit
560        quantity.limit()
561
562
563        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
564
565
566        #Assert that quantities are conserved
567        from Numeric import sum
568        for k in range(quantity.centroid_values.shape[0]):
569            assert allclose (quantity.centroid_values[k],
570                             sum(quantity.vertex_values[k,:])/3)
571
572
573
574
575
576    def test_distribute_first_order(self):
577        quantity = Conserved_quantity(self.mesh4)
578
579        #Test centroids
580        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
581        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
582
583
584        #Extrapolate
585        quantity.extrapolate_first_order()
586
587        #Interpolate
588        quantity.interpolate_from_vertices_to_edges()
589
590        assert allclose(quantity.vertex_values,
591                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
592        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
593                                               [3,3,3], [4, 4, 4]])
594
595
596    def test_distribute_second_order(self):
597        quantity = Conserved_quantity(self.mesh4)
598
599        #Test centroids
600        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
601        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
602
603
604        #Extrapolate
605        quantity.extrapolate_second_order()
606
607        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
608
609
610    def test_update_explicit(self):
611        quantity = Conserved_quantity(self.mesh4)
612
613        #Test centroids
614        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
615        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
616
617        #Set explicit_update
618        quantity.explicit_update = array( [1.,1.,1.,1.] )
619
620        #Update with given timestep
621        quantity.update(0.1)
622
623        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
624        assert allclose( quantity.centroid_values, x)
625
626    def test_update_semi_implicit(self):
627        quantity = Conserved_quantity(self.mesh4)
628
629        #Test centroids
630        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
631        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
632
633        #Set semi implicit update
634        quantity.semi_implicit_update = array([1.,1.,1.,1.])
635
636        #Update with given timestep
637        timestep = 0.1
638        quantity.update(timestep)
639
640        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
641        denom = ones(4, Float)-timestep*sem
642
643        x = array([1, 2, 3, 4])/denom
644        assert allclose( quantity.centroid_values, x)
645
646
647    def test_both_updates(self):
648        quantity = Conserved_quantity(self.mesh4)
649
650        #Test centroids
651        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
652        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
653
654        #Set explicit_update
655        quantity.explicit_update = array( [4.,3.,2.,1.] )
656
657        #Set semi implicit update
658        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
659
660        #Update with given timestep
661        timestep = 0.1
662        quantity.update(0.1)
663
664        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
665        denom = ones(4, Float)-timestep*sem
666
667        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
668        x /= denom
669        assert allclose( quantity.centroid_values, x)
670
671
672
673
674    #Test smoothing
675    def test_smoothing(self):
676
677        from mesh_factory import rectangular
678        from shallow_water import Domain, Transmissive_boundary
679        from Numeric import zeros, Float
680        from util import mean
681
682        #Create basic mesh
683        points, vertices, boundary = rectangular(2, 2)
684
685        #Create shallow water domain
686        domain = Domain(points, vertices, boundary)
687        domain.default_order=2
688        domain.reduction = mean
689
690
691        #Set some field values
692        domain.set_quantity('elevation', lambda x,y: x)
693        domain.set_quantity('friction', 0.03)
694
695
696        ######################
697        # Boundary conditions
698        B = Transmissive_boundary(domain)
699        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
700
701
702        ######################
703        #Initial condition - with jumps
704
705        bed = domain.quantities['elevation'].vertex_values
706        stage = zeros(bed.shape, Float)
707
708        h = 0.03
709        for i in range(stage.shape[0]):
710            if i % 2 == 0:
711                stage[i,:] = bed[i,:] + h
712            else:
713                stage[i,:] = bed[i,:]
714
715        domain.set_quantity('stage', stage)
716
717        stage = domain.quantities['stage']
718
719        #Get smoothed stage
720        A, V = stage.get_vertex_values(xy=False, smooth=True)
721        Q = stage.vertex_values
722
723
724        assert A.shape[0] == 9
725        assert V.shape[0] == 8
726        assert V.shape[1] == 3
727
728        #First four points
729        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
730        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
731        assert allclose(A[2], Q[3,0])
732        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
733
734        #Center point
735        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
736                               Q[5,0] + Q[6,2] + Q[7,1])/6)
737
738
739        #Check V
740        assert allclose(V[0,:], [3,4,0])
741        assert allclose(V[1,:], [1,0,4])
742        assert allclose(V[2,:], [4,5,1])
743        assert allclose(V[3,:], [2,1,5])
744        assert allclose(V[4,:], [6,7,3])
745        assert allclose(V[5,:], [4,3,7])
746        assert allclose(V[6,:], [7,8,4])
747        assert allclose(V[7,:], [5,4,8])
748
749        #Get smoothed stage with XY
750        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
751
752        assert allclose(A, A1)
753        assert allclose(V, V1)
754
755        #Check XY
756        assert allclose(X[4], 0.5)
757        assert allclose(Y[4], 0.5)
758
759        assert allclose(X[7], 1.0)
760        assert allclose(Y[7], 0.5)
761
762
763
764
765    def test_vertex_values_no_smoothing(self):
766
767        from mesh_factory import rectangular
768        from shallow_water import Domain, Transmissive_boundary
769        from Numeric import zeros, Float
770        from util import mean
771
772
773        #Create basic mesh
774        points, vertices, boundary = rectangular(2, 2)
775
776        #Create shallow water domain
777        domain = Domain(points, vertices, boundary)
778        domain.default_order=2
779        domain.reduction = mean
780
781
782        #Set some field values
783        domain.set_quantity('elevation', lambda x,y: x)
784        domain.set_quantity('friction', 0.03)
785
786
787        ######################
788        #Initial condition - with jumps
789
790        bed = domain.quantities['elevation'].vertex_values
791        stage = zeros(bed.shape, Float)
792
793        h = 0.03
794        for i in range(stage.shape[0]):
795            if i % 2 == 0:
796                stage[i,:] = bed[i,:] + h
797            else:
798                stage[i,:] = bed[i,:]
799
800        domain.set_quantity('stage', stage)
801
802        #Get stage
803        stage = domain.quantities['stage']
804        A, V = stage.get_vertex_values(xy=False, smooth=False)
805        Q = stage.vertex_values.flat
806
807        for k in range(8):
808            assert allclose(A[k], Q[k])
809
810
811        for k in range(8):
812            assert V[k, 0] == 3*k
813            assert V[k, 1] == 3*k+1
814            assert V[k, 2] == 3*k+2
815
816
817
818        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
819
820
821        assert allclose(A, A1)
822        assert allclose(V, V1)
823
824        #Check XY
825        assert allclose(X[1], 0.5)
826        assert allclose(Y[1], 0.5)
827        assert allclose(X[4], 0.0)
828        assert allclose(Y[4], 0.0)
829        assert allclose(X[12], 1.0)
830        assert allclose(Y[12], 0.0)
831
832
833
834    def set_array_values_by_index(self):
835
836        from mesh_factory import rectangular
837        from shallow_water import Domain
838        from Numeric import zeros, Float
839
840        #Create basic mesh
841        points, vertices, boundary = rectangular(1, 1)
842
843        #Create shallow water domain
844        domain = Domain(points, vertices, boundary)
845        #print "domain.number_of_elements ",domain.number_of_elements
846        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
847        value = [7]
848        indices = [1]
849        quantity.set_array_values_by_index(value,
850                                           location = 'centroids',
851                                           indices = indices)
852        #print "quantity.centroid_values",quantity.centroid_values
853
854        assert allclose(quantity.centroid_values, [1,7])
855
856        quantity.set_array_values([15,20,25], indices = indices)
857        assert allclose(quantity.centroid_values, [1,20])
858
859        quantity.set_array_values([15,20,25], indices = indices)
860        assert allclose(quantity.centroid_values, [1,20])
861
862    def test_setting_some_vertex_values(self):
863        """
864        set values based on triangle lists.
865        """
866        from mesh_factory import rectangular
867        from shallow_water import Domain
868        from Numeric import zeros, Float
869
870        #Create basic mesh
871        points, vertices, boundary = rectangular(1, 3)
872        #print "vertices",vertices
873        #Create shallow water domain
874        domain = Domain(points, vertices, boundary)
875        #print "domain.number_of_elements ",domain.number_of_elements
876        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
877                                    [4,4,4],[5,5,5],[6,6,6]])
878        value = [7]
879        indices = [1]
880        quantity.set_values(value,
881                            location = 'centroids',
882                            indices = indices)
883        #print "quantity.centroid_values",quantity.centroid_values
884        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
885
886        value = [[15,20,25]]
887        quantity.set_values(value, indices = indices)
888        #print "1 quantity.vertex_values",quantity.vertex_values
889        assert allclose(quantity.vertex_values[1], value[0])
890
891
892        #print "quantity",quantity.vertex_values
893        values = [10,100,50]
894        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
895        #print "2 quantity.vertex_values",quantity.vertex_values
896        assert allclose(quantity.vertex_values[0], [10,10,10])
897        assert allclose(quantity.vertex_values[5], [50,50,50])
898        #quantity.interpolate()
899        #print "quantity.centroid_values",quantity.centroid_values
900        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
901
902
903        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
904                                    [4,4,4],[5,5,5],[6,6,6]])
905        values = [10,100,50]
906        #this will be per unique vertex, indexing the vertices
907        #print "quantity.vertex_values",quantity.vertex_values
908        quantity.set_values(values, indices = [0,1,5])
909        #print "quantity.vertex_values",quantity.vertex_values
910        assert allclose(quantity.vertex_values[0], [1,50,10])
911        assert allclose(quantity.vertex_values[5], [6,6,6])
912        assert allclose(quantity.vertex_values[1], [100,10,50])
913
914        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
915                                    [4,4,4],[5,5,5],[6,6,6]])
916        values = [[31,30,29],[400,400,400],[1000,999,998]]
917        quantity.set_values(values, indices = [3,3,5])
918        quantity.interpolate()
919        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
920
921        values = [[1,1,1],[2,2,2],[3,3,3],
922                                    [4,4,4],[5,5,5],[6,6,6]]
923        quantity.set_values(values)
924
925        # testing the standard set values by vertex
926        # indexed by vertex_id in general_mesh.coordinates
927        values = [0,1,2,3,4,5,6,7]
928
929        quantity.set_values(values)
930        #print "1 quantity.vertex_values",quantity.vertex_values
931        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
932                                                [ 1.,  0.,  5.],
933                                                [ 5.,  6.,  1.],
934                                                [ 2.,  1.,  6.],
935                                                [ 6.,  7.,  2.],
936                                                [ 3.,  2.,  7.]])
937
938    def test_setting_unique_vertex_values(self):
939        """
940        set values based on unique_vertex lists.
941        """
942        from mesh_factory import rectangular
943        from shallow_water import Domain
944        from Numeric import zeros, Float
945
946        #Create basic mesh
947        points, vertices, boundary = rectangular(1, 3)
948        #print "vertices",vertices
949        #Create shallow water domain
950        domain = Domain(points, vertices, boundary)
951        #print "domain.number_of_elements ",domain.number_of_elements
952        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
953                                    [4,4,4],[5,5,5]])
954        value = 7
955        indices = [1,5]
956        quantity.set_values(value,
957                            location = 'unique vertices',
958                            indices = indices)
959        #print "quantity.centroid_values",quantity.centroid_values
960        assert allclose(quantity.vertex_values[0], [0,7,0])
961        assert allclose(quantity.vertex_values[1], [7,1,7])
962        assert allclose(quantity.vertex_values[2], [7,2,7])
963
964
965    def test_get_values(self):
966        """
967        get values based on triangle lists.
968        """
969        from mesh_factory import rectangular
970        from shallow_water import Domain
971        from Numeric import zeros, Float
972
973        #Create basic mesh
974        points, vertices, boundary = rectangular(1, 3)
975
976        #print "points",points
977        #print "vertices",vertices
978        #print "boundary",boundary
979
980        #Create shallow water domain
981        domain = Domain(points, vertices, boundary)
982        #print "domain.number_of_elements ",domain.number_of_elements
983        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
984                                    [4,4,4],[5,5,5]])
985
986        #print "quantity.get_values(location = 'unique vertices')", \
987        #      quantity.get_values(location = 'unique vertices')
988
989        #print "quantity.get_values(location = 'unique vertices')", \
990        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
991        #                          location = 'unique vertices')
992
993        answer = [0.5,2,4,5,0,1,3,4.5]
994        assert allclose(answer,
995                        quantity.get_values(location = 'unique vertices'))
996
997        indices = [0,5,3]
998        answer = [0.5,1,5]
999        assert allclose(answer,
1000                        quantity.get_values(indices=indices, \
1001                                            location = 'unique vertices'))
1002        #print "quantity.centroid_values",quantity.centroid_values
1003        #print "quantity.get_values(location = 'centroids') ",\
1004        #      quantity.get_values(location = 'centroids')
1005
1006    def test_getting_some_vertex_values(self):
1007        """
1008        get values based on triangle lists.
1009        """
1010        from mesh_factory import rectangular
1011        from shallow_water import Domain
1012        from Numeric import zeros, Float
1013
1014        #Create basic mesh
1015        points, vertices, boundary = rectangular(1, 3)
1016
1017        #print "points",points
1018        #print "vertices",vertices
1019        #print "boundary",boundary
1020
1021        #Create shallow water domain
1022        domain = Domain(points, vertices, boundary)
1023        #print "domain.number_of_elements ",domain.number_of_elements
1024        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1025                                    [4,4,4],[5,5,5],[6,6,6]])
1026        value = [7]
1027        indices = [1]
1028        quantity.set_values(value,
1029                            location = 'centroids',
1030                            indices = indices)
1031        #print "quantity.centroid_values",quantity.centroid_values
1032        #print "quantity.get_values(location = 'centroids') ",\
1033        #      quantity.get_values(location = 'centroids')
1034        assert allclose(quantity.centroid_values,
1035                        quantity.get_values(location = 'centroids'))
1036
1037
1038        value = [[15,20,25]]
1039        quantity.set_values(value, indices = indices)
1040        #print "1 quantity.vertex_values",quantity.vertex_values
1041        assert allclose(quantity.vertex_values, quantity.get_values())
1042
1043        assert allclose(quantity.edge_values,
1044                        quantity.get_values(location = 'edges'))
1045
1046        # get a subset of elements
1047        subset = quantity.get_values(location='centroids', indices=[0,5])
1048        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1049        assert allclose(subset, answer)
1050
1051
1052        subset = quantity.get_values(location='edges', indices=[0,5])
1053        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1054        #print "subset",subset
1055        #print "answer",answer
1056        assert allclose(subset, answer)
1057
1058        subset = quantity.get_values( indices=[1,5])
1059        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1060        #print "subset",subset
1061        #print "answer",answer
1062        assert allclose(subset, answer)
1063
1064
1065
1066#-------------------------------------------------------------
1067if __name__ == "__main__":
1068    suite = unittest.makeSuite(Test_Quantity,'test')
1069    #print "restricted test"
1070    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1071    runner = unittest.TextTestRunner()
1072    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.