source: inundation/pyvolution/test_quantity.py @ 1752

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

Refactored quantity.set_values

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