source: inundation/pyvolution/test_quantity.py @ 1916

Last change on this file since 1916 was 1916, checked in by ole, 18 years ago

Implemented operator overloading for (pow) in class Quantity and tested.
Wrote create_quantity_from_expression and test.

File size: 41.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
280    def test_set_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_set_values_from_quantity(self):
383
384        quantity1 = Quantity(self.mesh4)
385        quantity1.set_vertex_values([0,1,2,3,4,5])
386
387        assert allclose(quantity1.vertex_values,
388                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
389
390       
391        quantity2 = Quantity(self.mesh4)
392        quantity2.set_values(quantity = quantity1)
393        assert allclose(quantity2.vertex_values,
394                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
395
396        quantity2.set_values(quantity = 2*quantity1)
397        assert allclose(quantity2.vertex_values,
398                        [[2,0,4], [2,4,8], [8,4,10], [6,2,8]])
399       
400        quantity2.set_values(quantity = 2*quantity1 + 3)
401        assert allclose(quantity2.vertex_values,
402                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])
403
404
405        #Check detection of quantity as first orgument
406        quantity2.set_values(2*quantity1 + 3)
407        assert allclose(quantity2.vertex_values,
408                        [[5,3,7], [5,7,11], [11,7,13], [9,5,11]])             
409
410
411
412
413
414    def test_overloading(self):
415
416        quantity1 = Quantity(self.mesh4)
417        quantity1.set_vertex_values([0,1,2,3,4,5])
418
419        assert allclose(quantity1.vertex_values,
420                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
421
422       
423        quantity2 = Quantity(self.mesh4)
424        quantity2.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
425                             location = 'vertices')
426
427
428
429        quantity3 = Quantity(self.mesh4)
430        quantity3.set_values([[2,2,2], [7,8,9], [7,6,3], [3, 8, -8]],
431                             location = 'vertices')       
432
433
434        #Negation
435        Q = -quantity1
436        assert allclose(Q.vertex_values, -quantity1.vertex_values)
437        assert allclose(Q.centroid_values, -quantity1.centroid_values)
438        assert allclose(Q.edge_values, -quantity1.edge_values)               
439       
440        #Addition
441        Q = quantity1 + 7
442        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)
443        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
444        assert allclose(Q.edge_values, quantity1.edge_values + 7)       
445       
446        Q = 7 + quantity1
447        assert allclose(Q.vertex_values, quantity1.vertex_values + 7)       
448        assert allclose(Q.centroid_values, quantity1.centroid_values + 7)
449        assert allclose(Q.edge_values, quantity1.edge_values + 7)       
450       
451        Q = quantity1 + quantity2
452        assert allclose(Q.vertex_values,
453                        quantity1.vertex_values + quantity2.vertex_values)
454        assert allclose(Q.centroid_values,
455                        quantity1.centroid_values + quantity2.centroid_values)
456        assert allclose(Q.edge_values,
457                        quantity1.edge_values + quantity2.edge_values)       
458       
459
460        Q = quantity1 + quantity2 - 3
461        assert allclose(Q.vertex_values,
462                        quantity1.vertex_values + quantity2.vertex_values - 3)
463
464        Q = quantity1 - quantity2
465        assert allclose(Q.vertex_values,
466                        quantity1.vertex_values - quantity2.vertex_values)   
467
468        #Scaling
469        Q = quantity1*3
470        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
471        assert allclose(Q.centroid_values, quantity1.centroid_values*3)
472        assert allclose(Q.edge_values, quantity1.edge_values*3)               
473        Q = 3*quantity1
474        assert allclose(Q.vertex_values, quantity1.vertex_values*3)
475
476        #Multiplication
477        Q = quantity1 * quantity2
478        #print Q.vertex_values
479        #print Q.centroid_values
480        #print quantity1.centroid_values
481        #print quantity2.centroid_values
482       
483        assert allclose(Q.vertex_values,
484                        quantity1.vertex_values * quantity2.vertex_values)
485
486        #Linear combinations
487        Q = 4*quantity1 + 2
488        assert allclose(Q.vertex_values,
489                        4*quantity1.vertex_values + 2)
490       
491        Q = quantity1*quantity2 + 2
492        assert allclose(Q.vertex_values,
493                        quantity1.vertex_values * quantity2.vertex_values + 2)
494       
495        Q = quantity1*quantity2 + quantity3
496        assert allclose(Q.vertex_values,
497                        quantity1.vertex_values * quantity2.vertex_values +
498                        quantity3.vertex_values)       
499        Q = quantity1*quantity2 + 3*quantity3
500        assert allclose(Q.vertex_values,
501                        quantity1.vertex_values * quantity2.vertex_values +
502                        3*quantity3.vertex_values)               
503        Q = quantity1*quantity2 + 3*quantity3 + 5.0
504        assert allclose(Q.vertex_values,
505                        quantity1.vertex_values * quantity2.vertex_values +
506                        3*quantity3.vertex_values + 5)                       
507       
508        Q = quantity1*quantity2 - quantity3
509        assert allclose(Q.vertex_values,
510                        quantity1.vertex_values * quantity2.vertex_values -
511                        quantity3.vertex_values)                               
512        Q = 1.5*quantity1*quantity2 - 3*quantity3 + 5.0
513        assert allclose(Q.vertex_values,
514                        1.5*quantity1.vertex_values * quantity2.vertex_values -
515                        3*quantity3.vertex_values + 5)                               
516
517        #Try combining quantities and arrays and scalars
518        Q = 1.5*quantity1*quantity2.vertex_values -\
519            3*quantity3.vertex_values + 5.0 
520        assert allclose(Q.vertex_values,
521                        1.5*quantity1.vertex_values * quantity2.vertex_values -
522                        3*quantity3.vertex_values + 5)                                       
523       
524
525        #Powers
526        Q = quantity1**2
527        assert allclose(Q.vertex_values, quantity1.vertex_values**2)
528
529        Q = quantity1**2 +quantity2**2 
530        assert allclose(Q.vertex_values,
531                        quantity1.vertex_values**2 + \
532                        quantity2.vertex_values**2)
533
534        Q = (quantity1**2 +quantity2**2)**0.5 
535        assert allclose(Q.vertex_values,
536                        (quantity1.vertex_values**2 + \
537                        quantity2.vertex_values**2)**0.5)
538       
539
540
541       
542       
543       
544
545    def test_gradient(self):
546        quantity = Conserved_quantity(self.mesh4)
547
548        #Set up for a gradient of (2,0) at mid triangle
549        quantity.set_values([2.0, 4.0, 6.0, 2.0],
550                            location = 'centroids')
551
552
553        #Gradients
554        a, b = quantity.compute_gradients()
555
556        #print self.mesh4.centroid_coordinates
557        #print a, b
558
559        #The central triangle (1)
560        #(using standard gradient based on neigbours controid values)
561        assert allclose(a[1], 2.0)
562        assert allclose(b[1], 0.0)
563
564
565        #Left triangle (0) using two point gradient
566        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
567        #2  = 4  + a*(-2/3)  + b*(-2/3)
568        assert allclose(a[0] + b[0], 3)
569        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)       
570        assert allclose(a[0] - b[0], 0) 
571
572
573        #Right triangle (2) using two point gradient
574        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
575        #6  = 4  + a*(4/3)  + b*(-2/3)       
576        assert allclose(2*a[2] - b[2], 3)
577        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
578        assert allclose(a[2] + 2*b[2], 0) 
579
580
581        #Top triangle (3) using two point gradient
582        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
583        #2  = 4  + a*(-2/3)  + b*(4/3)       
584        assert allclose(a[3] - 2*b[3], 3)
585        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)       
586        assert allclose(2*a[3] + b[3], 0) 
587
588
589
590        #print a, b
591        quantity.extrapolate_second_order()
592
593        #Apply q(x,y) = qc + a*(x-xc) + b*(y-yc)
594        assert allclose(quantity.vertex_values[0,:], [3., 0.,  3.])
595        assert allclose(quantity.vertex_values[1,:], [4./3, 16./3,  16./3])
596
597
598        #a = 1.2, b=-0.6
599        #q(4,0) = 6 + a*(4 - 8/3) + b*(-2/3)
600        assert allclose(quantity.vertex_values[2,2], 8)       
601
602
603    def test_second_order_extrapolation2(self):
604        quantity = Conserved_quantity(self.mesh4)
605
606        #Set up for a gradient of (3,1), f(x) = 3x+y
607        quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
608                            location = 'centroids')
609
610        #Gradients
611        a, b = quantity.compute_gradients()
612
613        #print a, b
614
615        assert allclose(a[1], 3.0)
616        assert allclose(b[1], 1.0)
617       
618        #Work out the others
619       
620        quantity.extrapolate_second_order()
621
622        #print quantity.vertex_values
623        assert allclose(quantity.vertex_values[1,0], 2.0)
624        assert allclose(quantity.vertex_values[1,1], 6.0)
625        assert allclose(quantity.vertex_values[1,2], 8.0)
626
627
628
629    def test_first_order_extrapolator(self):
630        quantity = Conserved_quantity(self.mesh4)
631
632        #Test centroids
633        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
634        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
635
636        #Extrapolate
637        quantity.extrapolate_first_order()
638
639        #Check vertices but not edge values
640        assert allclose(quantity.vertex_values,
641                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
642
643
644    def test_second_order_extrapolator(self):
645        quantity = Conserved_quantity(self.mesh4)
646
647        #Set up for a gradient of (3,0) at mid triangle
648        quantity.set_values([2.0, 4.0, 8.0, 2.0],
649                            location = 'centroids')
650
651
652
653        quantity.extrapolate_second_order()
654        quantity.limit()
655
656
657        #Assert that central triangle is limited by neighbours
658        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
659        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
660
661        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
662        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
663
664        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
665        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
666
667
668        #Assert that quantities are conserved
669        from Numeric import sum
670        for k in range(quantity.centroid_values.shape[0]):
671            assert allclose (quantity.centroid_values[k],
672                             sum(quantity.vertex_values[k,:])/3)
673
674
675
676
677
678    def test_limiter(self):
679        quantity = Conserved_quantity(self.mesh4)
680
681        #Create a deliberate overshoot (e.g. from gradient computation)
682        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
683
684
685        #Limit
686        quantity.limit()
687
688        #Assert that central triangle is limited by neighbours
689        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
690        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
691
692        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
693        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
694
695        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
696        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
697
698
699
700        #Assert that quantities are conserved
701        from Numeric import sum
702        for k in range(quantity.centroid_values.shape[0]):
703            assert allclose (quantity.centroid_values[k],
704                             sum(quantity.vertex_values[k,:])/3)
705
706
707    def test_limiter2(self):
708        """Taken from test_shallow_water
709        """
710        quantity = Conserved_quantity(self.mesh4)
711
712        #Test centroids
713        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
714        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
715
716
717        #Extrapolate
718        quantity.extrapolate_second_order()
719
720        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
721
722        #Limit
723        quantity.limit()
724
725
726        assert allclose(quantity.vertex_values[1,:], [2.2, 4.9, 4.9])
727
728
729        #Assert that quantities are conserved
730        from Numeric import sum
731        for k in range(quantity.centroid_values.shape[0]):
732            assert allclose (quantity.centroid_values[k],
733                             sum(quantity.vertex_values[k,:])/3)
734
735
736
737
738
739    def test_distribute_first_order(self):
740        quantity = Conserved_quantity(self.mesh4)
741
742        #Test centroids
743        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
744        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
745
746
747        #Extrapolate
748        quantity.extrapolate_first_order()
749
750        #Interpolate
751        quantity.interpolate_from_vertices_to_edges()
752
753        assert allclose(quantity.vertex_values,
754                        [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
755        assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
756                                               [3,3,3], [4, 4, 4]])
757
758
759    def test_distribute_second_order(self):
760        quantity = Conserved_quantity(self.mesh4)
761
762        #Test centroids
763        quantity.set_values([2.,4.,8.,2.], location = 'centroids')
764        assert allclose(quantity.centroid_values, [2, 4, 8, 2]) #Centroid
765
766
767        #Extrapolate
768        quantity.extrapolate_second_order()
769
770        assert allclose(quantity.vertex_values[1,:], [0.0, 6, 6])
771
772
773    def test_update_explicit(self):
774        quantity = Conserved_quantity(self.mesh4)
775
776        #Test centroids
777        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
778        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
779
780        #Set explicit_update
781        quantity.explicit_update = array( [1.,1.,1.,1.] )
782
783        #Update with given timestep
784        quantity.update(0.1)
785
786        x = array([1, 2, 3, 4]) + array( [.1,.1,.1,.1] )
787        assert allclose( quantity.centroid_values, x)
788
789    def test_update_semi_implicit(self):
790        quantity = Conserved_quantity(self.mesh4)
791
792        #Test centroids
793        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
794        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
795
796        #Set semi implicit update
797        quantity.semi_implicit_update = array([1.,1.,1.,1.])
798
799        #Update with given timestep
800        timestep = 0.1
801        quantity.update(timestep)
802
803        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
804        denom = ones(4, Float)-timestep*sem
805
806        x = array([1, 2, 3, 4])/denom
807        assert allclose( quantity.centroid_values, x)
808
809
810    def test_both_updates(self):
811        quantity = Conserved_quantity(self.mesh4)
812
813        #Test centroids
814        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
815        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
816
817        #Set explicit_update
818        quantity.explicit_update = array( [4.,3.,2.,1.] )
819
820        #Set semi implicit update
821        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
822
823        #Update with given timestep
824        timestep = 0.1
825        quantity.update(0.1)
826
827        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
828        denom = ones(4, Float)-timestep*sem
829
830        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
831        x /= denom
832        assert allclose( quantity.centroid_values, x)
833
834
835
836
837    #Test smoothing
838    def test_smoothing(self):
839
840        from mesh_factory import rectangular
841        from shallow_water import Domain, Transmissive_boundary
842        from Numeric import zeros, Float
843        from util import mean
844
845        #Create basic mesh
846        points, vertices, boundary = rectangular(2, 2)
847
848        #Create shallow water domain
849        domain = Domain(points, vertices, boundary)
850        domain.default_order=2
851        domain.reduction = mean
852
853
854        #Set some field values
855        domain.set_quantity('elevation', lambda x,y: x)
856        domain.set_quantity('friction', 0.03)
857
858
859        ######################
860        # Boundary conditions
861        B = Transmissive_boundary(domain)
862        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
863
864
865        ######################
866        #Initial condition - with jumps
867
868        bed = domain.quantities['elevation'].vertex_values
869        stage = zeros(bed.shape, Float)
870
871        h = 0.03
872        for i in range(stage.shape[0]):
873            if i % 2 == 0:
874                stage[i,:] = bed[i,:] + h
875            else:
876                stage[i,:] = bed[i,:]
877
878        domain.set_quantity('stage', stage)
879
880        stage = domain.quantities['stage']
881
882        #Get smoothed stage
883        A, V = stage.get_vertex_values(xy=False, smooth=True)
884        Q = stage.vertex_values
885
886
887        assert A.shape[0] == 9
888        assert V.shape[0] == 8
889        assert V.shape[1] == 3
890
891        #First four points
892        assert allclose(A[0], (Q[0,2] + Q[1,1])/2)
893        assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)
894        assert allclose(A[2], Q[3,0])
895        assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)
896
897        #Center point
898        assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\
899                               Q[5,0] + Q[6,2] + Q[7,1])/6)
900
901
902        #Check V
903        assert allclose(V[0,:], [3,4,0])
904        assert allclose(V[1,:], [1,0,4])
905        assert allclose(V[2,:], [4,5,1])
906        assert allclose(V[3,:], [2,1,5])
907        assert allclose(V[4,:], [6,7,3])
908        assert allclose(V[5,:], [4,3,7])
909        assert allclose(V[6,:], [7,8,4])
910        assert allclose(V[7,:], [5,4,8])
911
912        #Get smoothed stage with XY
913        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)
914
915        assert allclose(A, A1)
916        assert allclose(V, V1)
917
918        #Check XY
919        assert allclose(X[4], 0.5)
920        assert allclose(Y[4], 0.5)
921
922        assert allclose(X[7], 1.0)
923        assert allclose(Y[7], 0.5)
924
925
926
927
928    def test_vertex_values_no_smoothing(self):
929
930        from mesh_factory import rectangular
931        from shallow_water import Domain, Transmissive_boundary
932        from Numeric import zeros, Float
933        from util import mean
934
935
936        #Create basic mesh
937        points, vertices, boundary = rectangular(2, 2)
938
939        #Create shallow water domain
940        domain = Domain(points, vertices, boundary)
941        domain.default_order=2
942        domain.reduction = mean
943
944
945        #Set some field values
946        domain.set_quantity('elevation', lambda x,y: x)
947        domain.set_quantity('friction', 0.03)
948
949
950        ######################
951        #Initial condition - with jumps
952
953        bed = domain.quantities['elevation'].vertex_values
954        stage = zeros(bed.shape, Float)
955
956        h = 0.03
957        for i in range(stage.shape[0]):
958            if i % 2 == 0:
959                stage[i,:] = bed[i,:] + h
960            else:
961                stage[i,:] = bed[i,:]
962
963        domain.set_quantity('stage', stage)
964
965        #Get stage
966        stage = domain.quantities['stage']
967        A, V = stage.get_vertex_values(xy=False, smooth=False)
968        Q = stage.vertex_values.flat
969
970        for k in range(8):
971            assert allclose(A[k], Q[k])
972
973
974        for k in range(8):
975            assert V[k, 0] == 3*k
976            assert V[k, 1] == 3*k+1
977            assert V[k, 2] == 3*k+2
978
979
980
981        X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)
982
983
984        assert allclose(A, A1)
985        assert allclose(V, V1)
986
987        #Check XY
988        assert allclose(X[1], 0.5)
989        assert allclose(Y[1], 0.5)
990        assert allclose(X[4], 0.0)
991        assert allclose(Y[4], 0.0)
992        assert allclose(X[12], 1.0)
993        assert allclose(Y[12], 0.0)
994
995
996
997    def set_array_values_by_index(self):
998
999        from mesh_factory import rectangular
1000        from shallow_water import Domain
1001        from Numeric import zeros, Float
1002
1003        #Create basic mesh
1004        points, vertices, boundary = rectangular(1, 1)
1005
1006        #Create shallow water domain
1007        domain = Domain(points, vertices, boundary)
1008        #print "domain.number_of_elements ",domain.number_of_elements
1009        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
1010        value = [7]
1011        indices = [1]
1012        quantity.set_array_values_by_index(value,
1013                                           location = 'centroids',
1014                                           indices = indices)
1015        #print "quantity.centroid_values",quantity.centroid_values
1016
1017        assert allclose(quantity.centroid_values, [1,7])
1018
1019        quantity.set_array_values([15,20,25], indices = indices)
1020        assert allclose(quantity.centroid_values, [1,20])
1021
1022        quantity.set_array_values([15,20,25], indices = indices)
1023        assert allclose(quantity.centroid_values, [1,20])
1024
1025    def test_setting_some_vertex_values(self):
1026        """
1027        set values based on triangle lists.
1028        """
1029        from mesh_factory import rectangular
1030        from shallow_water import Domain
1031        from Numeric import zeros, Float
1032
1033        #Create basic mesh
1034        points, vertices, boundary = rectangular(1, 3)
1035        #print "vertices",vertices
1036        #Create shallow water domain
1037        domain = Domain(points, vertices, boundary)
1038        #print "domain.number_of_elements ",domain.number_of_elements
1039        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1040                                    [4,4,4],[5,5,5],[6,6,6]])
1041        value = [7]
1042        indices = [1]
1043        quantity.set_values(value,
1044                            location = 'centroids',
1045                            indices = indices)
1046        #print "quantity.centroid_values",quantity.centroid_values
1047        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
1048
1049        value = [[15,20,25]]
1050        quantity.set_values(value, indices = indices)
1051        #print "1 quantity.vertex_values",quantity.vertex_values
1052        assert allclose(quantity.vertex_values[1], value[0])
1053
1054
1055        #print "quantity",quantity.vertex_values
1056        values = [10,100,50]
1057        quantity.set_values(values, indices = [0,1,5], location = 'centroids')
1058        #print "2 quantity.vertex_values",quantity.vertex_values
1059        assert allclose(quantity.vertex_values[0], [10,10,10])
1060        assert allclose(quantity.vertex_values[5], [50,50,50])
1061        #quantity.interpolate()
1062        #print "quantity.centroid_values",quantity.centroid_values
1063        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
1064
1065
1066        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1067                                    [4,4,4],[5,5,5],[6,6,6]])
1068        values = [10,100,50]
1069        #this will be per unique vertex, indexing the vertices
1070        #print "quantity.vertex_values",quantity.vertex_values
1071        quantity.set_values(values, indices = [0,1,5])
1072        #print "quantity.vertex_values",quantity.vertex_values
1073        assert allclose(quantity.vertex_values[0], [1,50,10])
1074        assert allclose(quantity.vertex_values[5], [6,6,6])
1075        assert allclose(quantity.vertex_values[1], [100,10,50])
1076
1077        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1078                                    [4,4,4],[5,5,5],[6,6,6]])
1079        values = [[31,30,29],[400,400,400],[1000,999,998]]
1080        quantity.set_values(values, indices = [3,3,5])
1081        quantity.interpolate()
1082        assert allclose(quantity.centroid_values, [1,2,3,400,5,999])
1083
1084        values = [[1,1,1],[2,2,2],[3,3,3],
1085                                    [4,4,4],[5,5,5],[6,6,6]]
1086        quantity.set_values(values)
1087
1088        # testing the standard set values by vertex
1089        # indexed by vertex_id in general_mesh.coordinates
1090        values = [0,1,2,3,4,5,6,7]
1091
1092        quantity.set_values(values)
1093        #print "1 quantity.vertex_values",quantity.vertex_values
1094        assert allclose(quantity.vertex_values,[[ 4.,  5.,  0.],
1095                                                [ 1.,  0.,  5.],
1096                                                [ 5.,  6.,  1.],
1097                                                [ 2.,  1.,  6.],
1098                                                [ 6.,  7.,  2.],
1099                                                [ 3.,  2.,  7.]])
1100
1101    def test_setting_unique_vertex_values(self):
1102        """
1103        set values based on unique_vertex lists.
1104        """
1105        from mesh_factory import rectangular
1106        from shallow_water import Domain
1107        from Numeric import zeros, Float
1108
1109        #Create basic mesh
1110        points, vertices, boundary = rectangular(1, 3)
1111        #print "vertices",vertices
1112        #Create shallow water domain
1113        domain = Domain(points, vertices, boundary)
1114        #print "domain.number_of_elements ",domain.number_of_elements
1115        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1116                                    [4,4,4],[5,5,5]])
1117        value = 7
1118        indices = [1,5]
1119        quantity.set_values(value,
1120                            location = 'unique vertices',
1121                            indices = indices)
1122        #print "quantity.centroid_values",quantity.centroid_values
1123        assert allclose(quantity.vertex_values[0], [0,7,0])
1124        assert allclose(quantity.vertex_values[1], [7,1,7])
1125        assert allclose(quantity.vertex_values[2], [7,2,7])
1126
1127
1128    def test_get_values(self):
1129        """
1130        get values based on triangle lists.
1131        """
1132        from mesh_factory import rectangular
1133        from shallow_water import Domain
1134        from Numeric import zeros, Float
1135
1136        #Create basic mesh
1137        points, vertices, boundary = rectangular(1, 3)
1138
1139        #print "points",points
1140        #print "vertices",vertices
1141        #print "boundary",boundary
1142
1143        #Create shallow water domain
1144        domain = Domain(points, vertices, boundary)
1145        #print "domain.number_of_elements ",domain.number_of_elements
1146        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
1147                                    [4,4,4],[5,5,5]])
1148
1149        #print "quantity.get_values(location = 'unique vertices')", \
1150        #      quantity.get_values(location = 'unique vertices')
1151
1152        #print "quantity.get_values(location = 'unique vertices')", \
1153        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
1154        #                          location = 'unique vertices')
1155
1156        answer = [0.5,2,4,5,0,1,3,4.5]
1157        assert allclose(answer,
1158                        quantity.get_values(location = 'unique vertices'))
1159
1160        indices = [0,5,3]
1161        answer = [0.5,1,5]
1162        assert allclose(answer,
1163                        quantity.get_values(indices=indices, \
1164                                            location = 'unique vertices'))
1165        #print "quantity.centroid_values",quantity.centroid_values
1166        #print "quantity.get_values(location = 'centroids') ",\
1167        #      quantity.get_values(location = 'centroids')
1168
1169    def test_getting_some_vertex_values(self):
1170        """
1171        get values based on triangle lists.
1172        """
1173        from mesh_factory import rectangular
1174        from shallow_water import Domain
1175        from Numeric import zeros, Float
1176
1177        #Create basic mesh
1178        points, vertices, boundary = rectangular(1, 3)
1179
1180        #print "points",points
1181        #print "vertices",vertices
1182        #print "boundary",boundary
1183
1184        #Create shallow water domain
1185        domain = Domain(points, vertices, boundary)
1186        #print "domain.number_of_elements ",domain.number_of_elements
1187        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
1188                                    [4,4,4],[5,5,5],[6,6,6]])
1189        value = [7]
1190        indices = [1]
1191        quantity.set_values(value,
1192                            location = 'centroids',
1193                            indices = indices)
1194        #print "quantity.centroid_values",quantity.centroid_values
1195        #print "quantity.get_values(location = 'centroids') ",\
1196        #      quantity.get_values(location = 'centroids')
1197        assert allclose(quantity.centroid_values,
1198                        quantity.get_values(location = 'centroids'))
1199
1200
1201        value = [[15,20,25]]
1202        quantity.set_values(value, indices = indices)
1203        #print "1 quantity.vertex_values",quantity.vertex_values
1204        assert allclose(quantity.vertex_values, quantity.get_values())
1205
1206        assert allclose(quantity.edge_values,
1207                        quantity.get_values(location = 'edges'))
1208
1209        # get a subset of elements
1210        subset = quantity.get_values(location='centroids', indices=[0,5])
1211        answer = [quantity.centroid_values[0],quantity.centroid_values[5]]
1212        assert allclose(subset, answer)
1213
1214
1215        subset = quantity.get_values(location='edges', indices=[0,5])
1216        answer = [quantity.edge_values[0],quantity.edge_values[5]]
1217        #print "subset",subset
1218        #print "answer",answer
1219        assert allclose(subset, answer)
1220
1221        subset = quantity.get_values( indices=[1,5])
1222        answer = [quantity.vertex_values[1],quantity.vertex_values[5]]
1223        #print "subset",subset
1224        #print "answer",answer
1225        assert allclose(subset, answer)
1226
1227
1228
1229
1230#-------------------------------------------------------------
1231if __name__ == "__main__":
1232    suite = unittest.makeSuite(Test_Quantity,'test')
1233    #print "restricted test"
1234    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
1235    runner = unittest.TextTestRunner()
1236    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.