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

Last change on this file since 1697 was 1657, checked in by ole, 20 years ago

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