source: branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 15 years ago

Initial commit of numpy changes. Still a long way to go.

File size: 24.9 KB
Line 
1#!/usr/bin/env python
2
3import unittest
4from math import sqrt
5
6from anuga.abstract_2d_finite_volumes.domain import *
7from anuga.config import epsilon
8
9import numpy as num
10
11
12def add_to_verts(tag, elements, domain):
13    if tag == "mound":
14        domain.test = "Mound"
15
16
17def set_bottom_friction(tag, elements, domain):
18    if tag == "bottom":
19        #print 'bottom - indices',elements
20        domain.set_quantity('friction', 0.09, indices = elements)
21
22def set_top_friction(tag, elements, domain):
23    if tag == "top":
24        #print 'top - indices',elements
25        domain.set_quantity('friction', 1., indices = elements)
26
27
28def set_all_friction(tag, elements, domain):
29    if tag == 'all':
30        new_values = domain.get_quantity('friction').get_values(indices = elements) + 10.0
31
32        domain.set_quantity('friction', new_values, indices = elements)
33
34
35class Test_Domain(unittest.TestCase):
36    def setUp(self):
37        pass
38
39
40    def tearDown(self):
41        pass
42
43
44    def test_simple(self):
45        a = [0.0, 0.0]
46        b = [0.0, 2.0]
47        c = [2.0,0.0]
48        d = [0.0, 4.0]
49        e = [2.0, 2.0]
50        f = [4.0,0.0]
51
52        points = [a, b, c, d, e, f]
53        #bac, bce, ecf, dbe, daf, dae
54        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
55
56        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
57        other_quantities = ['elevation', 'friction']
58
59        domain = Domain(points, vertices, None,
60                        conserved_quantities, other_quantities)
61        domain.check_integrity()
62
63        for name in conserved_quantities + other_quantities:
64            assert domain.quantities.has_key(name)
65
66
67        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
68
69
70    def test_conserved_quantities(self):
71
72        a = [0.0, 0.0]
73        b = [0.0, 2.0]
74        c = [2.0,0.0]
75        d = [0.0, 4.0]
76        e = [2.0, 2.0]
77        f = [4.0,0.0]
78
79        points = [a, b, c, d, e, f]
80        #bac, bce, ecf, dbe, daf, dae
81        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
82
83        domain = Domain(points, vertices, boundary=None,
84                        conserved_quantities =\
85                        ['stage', 'xmomentum', 'ymomentum'])
86
87
88        domain.set_quantity('stage', [[1,2,3], [5,5,5],
89                                      [0,0,9], [-6, 3, 3]])
90
91        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
92                                          [0,0,9], [-6, 3, 3]])
93
94        domain.check_integrity()
95
96        #Centroids
97        q = domain.get_conserved_quantities(0)
98        assert num.allclose(q, [2., 2., 0.])
99
100        q = domain.get_conserved_quantities(1)
101        assert num.allclose(q, [5., 5., 0.])
102
103        q = domain.get_conserved_quantities(2)
104        assert num.allclose(q, [3., 3., 0.])
105
106        q = domain.get_conserved_quantities(3)
107        assert num.allclose(q, [0., 0., 0.])
108
109
110        #Edges
111        q = domain.get_conserved_quantities(0, edge=0)
112        assert num.allclose(q, [2.5, 2.5, 0.])
113        q = domain.get_conserved_quantities(0, edge=1)
114        assert num.allclose(q, [2., 2., 0.])
115        q = domain.get_conserved_quantities(0, edge=2)
116        assert num.allclose(q, [1.5, 1.5, 0.])
117
118        for i in range(3):
119            q = domain.get_conserved_quantities(1, edge=i)
120            assert num.allclose(q, [5, 5, 0.])
121
122
123        q = domain.get_conserved_quantities(2, edge=0)
124        assert num.allclose(q, [4.5, 4.5, 0.])
125        q = domain.get_conserved_quantities(2, edge=1)
126        assert num.allclose(q, [4.5, 4.5, 0.])
127        q = domain.get_conserved_quantities(2, edge=2)
128        assert num.allclose(q, [0., 0., 0.])
129
130
131        q = domain.get_conserved_quantities(3, edge=0)
132        assert num.allclose(q, [3., 3., 0.])
133        q = domain.get_conserved_quantities(3, edge=1)
134        assert num.allclose(q, [-1.5, -1.5, 0.])
135        q = domain.get_conserved_quantities(3, edge=2)
136        assert num.allclose(q, [-1.5, -1.5, 0.])
137
138
139
140    def test_create_quantity_from_expression(self):
141        """Quantity created from other quantities using arbitrary expression
142
143        """
144
145
146        a = [0.0, 0.0]
147        b = [0.0, 2.0]
148        c = [2.0,0.0]
149        d = [0.0, 4.0]
150        e = [2.0, 2.0]
151        f = [4.0,0.0]
152
153        points = [a, b, c, d, e, f]
154        #bac, bce, ecf, dbe, daf, dae
155        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
156
157        domain = Domain(points, vertices, boundary=None,
158                        conserved_quantities =\
159                        ['stage', 'xmomentum', 'ymomentum'],
160                        other_quantities = ['elevation', 'friction'])
161
162
163        domain.set_quantity('elevation', -1)
164
165
166        domain.set_quantity('stage', [[1,2,3], [5,5,5],
167                                      [0,0,9], [-6, 3, 3]])
168
169        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
170                                          [0,0,9], [-6, 3, 3]])
171
172        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
173                                          [2,4,-1], [1, 0, 1]])
174
175        domain.check_integrity()
176
177
178
179        expression = 'stage - elevation'
180        Q = domain.create_quantity_from_expression(expression)
181
182        assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
183                                              [1,1,10], [-5, 4, 4]])
184
185        expression = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'
186        Q = domain.create_quantity_from_expression(expression)
187
188        X = domain.quantities['xmomentum'].vertex_values
189        Y = domain.quantities['ymomentum'].vertex_values
190
191        assert num.allclose(Q.vertex_values, (X**2 + Y**2)**0.5)
192
193
194
195    def test_set_quanitities_to_be_monitored(self):
196        """test_set_quanitities_to_be_monitored
197        """
198
199        a = [0.0, 0.0]
200        b = [0.0, 2.0]
201        c = [2.0,0.0]
202        d = [0.0, 4.0]
203        e = [2.0, 2.0]
204        f = [4.0,0.0]
205
206        points = [a, b, c, d, e, f]
207        #bac, bce, ecf, dbe, daf, dae
208        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
209
210
211        domain = Domain(points, vertices, boundary=None,
212                        conserved_quantities =\
213                        ['stage', 'xmomentum', 'ymomentum'],
214                        other_quantities = ['elevation', 'friction', 'depth'])
215
216
217        assert domain.quantities_to_be_monitored is None
218        domain.set_quantities_to_be_monitored(['stage', 'stage-elevation'])
219        assert len(domain.quantities_to_be_monitored) == 2
220        assert domain.quantities_to_be_monitored.has_key('stage')
221        assert domain.quantities_to_be_monitored.has_key('stage-elevation')
222        for key in domain.quantities_to_be_monitored['stage'].keys():
223            assert domain.quantities_to_be_monitored['stage'][key] is None
224
225
226        # Check that invalid requests are dealt with
227        try:
228            domain.set_quantities_to_be_monitored(['yyyyy'])       
229        except:
230            pass
231        else:
232            msg = 'Should have caught illegal quantity'
233            raise Exception, msg
234
235        try:
236            domain.set_quantities_to_be_monitored(['stage-xx'])       
237        except NameError:
238            pass
239        else:
240            msg = 'Should have caught illegal quantity'
241            raise Exception, msg
242
243        try:
244            domain.set_quantities_to_be_monitored('stage', 'stage-elevation')
245        except:
246            pass
247        else:
248            msg = 'Should have caught too many arguments'
249            raise Exception, msg
250
251        try:
252            domain.set_quantities_to_be_monitored('stage', 'blablabla')
253        except:
254            pass
255        else:
256            msg = 'Should have caught polygon as a string'
257            raise Exception, msg       
258
259
260
261        # Now try with a polygon restriction
262        domain.set_quantities_to_be_monitored('xmomentum',
263                                              polygon=[[1,1], [1,3], [3,3], [3,1]],
264                                              time_interval = [0,3])
265        assert domain.monitor_indices[0] == 1
266        assert domain.monitor_time_interval[0] == 0
267        assert domain.monitor_time_interval[1] == 3       
268       
269
270    def test_set_quantity_from_expression(self):
271        """Quantity set using arbitrary expression
272
273        """
274
275
276        a = [0.0, 0.0]
277        b = [0.0, 2.0]
278        c = [2.0,0.0]
279        d = [0.0, 4.0]
280        e = [2.0, 2.0]
281        f = [4.0,0.0]
282
283        points = [a, b, c, d, e, f]
284        #bac, bce, ecf, dbe, daf, dae
285        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
286
287        domain = Domain(points, vertices, boundary=None,
288                        conserved_quantities =\
289                        ['stage', 'xmomentum', 'ymomentum'],
290                        other_quantities = ['elevation', 'friction', 'depth'])
291
292
293        domain.set_quantity('elevation', -1)
294
295
296        domain.set_quantity('stage', [[1,2,3], [5,5,5],
297                                      [0,0,9], [-6, 3, 3]])
298
299        domain.set_quantity('xmomentum', [[1,2,3], [5,5,5],
300                                          [0,0,9], [-6, 3, 3]])
301
302        domain.set_quantity('ymomentum', [[3,3,3], [4,2,1],
303                                          [2,4,-1], [1, 0, 1]])
304
305
306
307
308        domain.set_quantity('depth', expression = 'stage - elevation')
309
310        domain.check_integrity()
311
312
313
314
315        Q = domain.quantities['depth']
316
317        assert num.allclose(Q.vertex_values, [[2,3,4], [6,6,6],
318                                              [1,1,10], [-5, 4, 4]])
319
320
321
322                                     
323    def test_add_quantity(self):
324        """Test that quantities already set can be added to using
325        add_quantity
326
327        """
328
329
330        a = [0.0, 0.0]
331        b = [0.0, 2.0]
332        c = [2.0,0.0]
333        d = [0.0, 4.0]
334        e = [2.0, 2.0]
335        f = [4.0,0.0]
336
337        points = [a, b, c, d, e, f]
338        #bac, bce, ecf, dbe, daf, dae
339        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
340
341        domain = Domain(points, vertices, boundary=None,
342                        conserved_quantities =\
343                        ['stage', 'xmomentum', 'ymomentum'],
344                        other_quantities = ['elevation', 'friction', 'depth'])
345
346
347        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.float)
348        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.float)
349       
350        #print A
351        #print B
352        #print A+B
353       
354        # Shorthands
355        stage = domain.quantities['stage']
356        elevation = domain.quantities['elevation']
357        depth = domain.quantities['depth']
358       
359        # Go testing
360        domain.set_quantity('elevation', A)
361        domain.add_quantity('elevation', B)
362        assert num.allclose(elevation.vertex_values, A+B)
363       
364        domain.add_quantity('elevation', 4)
365        assert num.allclose(elevation.vertex_values, A+B+4)       
366       
367       
368        # Test using expression
369        domain.set_quantity('stage', [[1,2,3], [5,5,5],
370                                      [0,0,9], [-6, 3, 3]])       
371        domain.set_quantity('depth', 1.0)                                     
372        domain.add_quantity('depth', expression = 'stage - elevation')       
373        assert num.allclose(depth.vertex_values, stage.vertex_values-elevation.vertex_values+1)
374               
375       
376        # Check self referential expression
377        reference = 2*stage.vertex_values - depth.vertex_values
378        domain.add_quantity('stage', expression = 'stage - depth')               
379        assert num.allclose(stage.vertex_values, reference)       
380                                     
381
382        # Test using a function
383        def f(x, y):
384            return x+y
385           
386        domain.set_quantity('elevation', f)           
387        domain.set_quantity('stage', 5.0)
388        domain.set_quantity('depth', expression = 'stage - elevation')
389       
390        domain.add_quantity('depth', f)
391        assert num.allclose(stage.vertex_values, depth.vertex_values)               
392         
393           
394       
395       
396                                     
397                                     
398    def test_setting_timestepping_method(self):
399        """test_set_quanitities_to_be_monitored
400        """
401
402        a = [0.0, 0.0]
403        b = [0.0, 2.0]
404        c = [2.0,0.0]
405        d = [0.0, 4.0]
406        e = [2.0, 2.0]
407        f = [4.0,0.0]
408
409        points = [a, b, c, d, e, f]
410        #bac, bce, ecf, dbe, daf, dae
411        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
412
413
414        domain = Domain(points, vertices, boundary=None,
415                        conserved_quantities =\
416                        ['stage', 'xmomentum', 'ymomentum'],
417                        other_quantities = ['elevation', 'friction', 'depth'])
418
419
420        domain.timestepping_method = None
421
422
423        # Check that invalid requests are dealt with
424        try:
425            domain.set_timestepping_method('eee')       
426        except:
427            pass
428        else:
429            msg = 'Should have caught illegal method'
430            raise Exception, msg
431
432
433        #Should have no trouble with euler, rk2 or rk3
434        domain.set_timestepping_method('euler')
435        domain.set_timestepping_method('rk2')
436        domain.set_timestepping_method('rk3')
437
438        #test get timestepping method
439        assert domain.get_timestepping_method() == 'rk3'
440
441
442
443    def test_boundary_indices(self):
444
445        from anuga.config import default_boundary_tag
446
447
448        a = [0.0, 0.5]
449        b = [0.0, 0.0]
450        c = [0.5, 0.5]
451
452        points = [a, b, c]
453        vertices = [ [0,1,2] ]
454        domain = Domain(points, vertices)
455
456        domain.set_boundary( {default_boundary_tag: Dirichlet_boundary([5,2,1])} )
457
458
459        domain.check_integrity()
460
461        assert num.allclose(domain.neighbours, [[-1,-2,-3]])
462
463
464
465    def test_boundary_conditions(self):
466
467        a = [0.0, 0.0]
468        b = [0.0, 2.0]
469        c = [2.0,0.0]
470        d = [0.0, 4.0]
471        e = [2.0, 2.0]
472        f = [4.0,0.0]
473
474        points = [a, b, c, d, e, f]
475        #bac, bce, ecf, dbe
476        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
477        boundary = { (0, 0): 'First',
478                     (0, 2): 'First',
479                     (2, 0): 'Second',
480                     (2, 1): 'Second',
481                     (3, 1): 'Second',
482                     (3, 2): 'Second'}
483
484
485        domain = Domain(points, vertices, boundary,
486                        conserved_quantities =\
487                        ['stage', 'xmomentum', 'ymomentum'])
488        domain.check_integrity()
489
490
491
492        domain.set_quantity('stage', [[1,2,3], [5,5,5],
493                                      [0,0,9], [-6, 3, 3]])
494
495
496        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]),
497                              'Second': Transmissive_boundary(domain)} )
498
499        domain.update_boundary()
500
501        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
502        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
503        assert domain.quantities['stage'].boundary_values[2] ==\
504               domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
505        assert domain.quantities['stage'].boundary_values[3] ==\
506               domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
507        assert domain.quantities['stage'].boundary_values[4] ==\
508               domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
509        assert domain.quantities['stage'].boundary_values[5] ==\
510               domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5)
511
512        #Check enumeration
513        for k, ((vol_id, edge_id), _) in enumerate(domain.boundary_objects):
514            assert domain.neighbours[vol_id, edge_id] == -k-1
515
516
517
518
519    def test_distribute_first_order(self):
520        """Domain implements a default first order gradient limiter
521        """
522
523        a = [0.0, 0.0]
524        b = [0.0, 2.0]
525        c = [2.0,0.0]
526        d = [0.0, 4.0]
527        e = [2.0, 2.0]
528        f = [4.0,0.0]
529
530        points = [a, b, c, d, e, f]
531        #bac, bce, ecf, dbe
532        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
533        boundary = { (0, 0): 'Third',
534                     (0, 2): 'First',
535                     (2, 0): 'Second',
536                     (2, 1): 'Second',
537                     (3, 1): 'Second',
538                     (3, 2): 'Third'}
539
540
541        domain = Domain(points, vertices, boundary,
542                        conserved_quantities =\
543                        ['stage', 'xmomentum', 'ymomentum'])
544        domain.check_integrity()
545
546
547        domain.set_quantity('stage', [[1,2,3], [5,5,5],
548                                      [0,0,9], [-6, 3, 3]])
549
550        assert num.allclose( domain.quantities['stage'].centroid_values,
551                             [2,5,3,0] )
552
553        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
554                                          [3,3,3], [4, 4, 4]])
555
556        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
557                                          [30,30,30], [40, 40, 40]])
558
559
560        domain.distribute_to_vertices_and_edges()
561
562        #First order extrapolation
563        assert num.allclose( domain.quantities['stage'].vertex_values,
564                             [[ 2.,  2.,  2.],
565                              [ 5.,  5.,  5.],
566                              [ 3.,  3.,  3.],
567                              [ 0.,  0.,  0.]])
568
569
570
571
572    def test_update_conserved_quantities(self):
573        a = [0.0, 0.0]
574        b = [0.0, 2.0]
575        c = [2.0,0.0]
576        d = [0.0, 4.0]
577        e = [2.0, 2.0]
578        f = [4.0,0.0]
579
580        points = [a, b, c, d, e, f]
581        #bac, bce, ecf, dbe
582        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
583        boundary = { (0, 0): 'Third',
584                     (0, 2): 'First',
585                     (2, 0): 'Second',
586                     (2, 1): 'Second',
587                     (3, 1): 'Second',
588                     (3, 2): 'Third'}
589
590
591        domain = Domain(points, vertices, boundary,
592                        conserved_quantities =\
593                        ['stage', 'xmomentum', 'ymomentum'])
594        domain.check_integrity()
595
596
597        domain.set_quantity('stage', [1,2,3,4], location='centroids')
598        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
599        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
600
601
602        #Assign some values to update vectors
603        #Set explicit_update
604
605        for name in domain.conserved_quantities:
606            domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.])
607            domain.quantities[name].semi_implicit_update = num.array([1.,1.,1.,1.])
608
609
610        #Update with given timestep (assuming no other forcing terms)
611        domain.timestep = 0.1
612        domain.update_conserved_quantities()
613
614        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
615        denom = num.ones(4, num.float)-domain.timestep*sem
616
617#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
618#        x /= denom
619
620        x = num.array([1., 2., 3., 4.])
621        x /= denom
622        x += domain.timestep*num.array( [4,3,2,1] )
623
624        for name in domain.conserved_quantities:
625            assert num.allclose(domain.quantities[name].centroid_values, x)
626
627
628    def test_set_region(self):
629        """Set quantities for sub region
630        """
631
632        a = [0.0, 0.0]
633        b = [0.0, 2.0]
634        c = [2.0,0.0]
635        d = [0.0, 4.0]
636        e = [2.0, 2.0]
637        f = [4.0,0.0]
638
639        points = [a, b, c, d, e, f]
640        #bac, bce, ecf, dbe
641        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
642        boundary = { (0, 0): 'Third',
643                     (0, 2): 'First',
644                     (2, 0): 'Second',
645                     (2, 1): 'Second',
646                     (3, 1): 'Second',
647                     (3, 2): 'Third'}
648
649        domain = Domain(points, vertices, boundary,
650                        conserved_quantities =\
651                        ['stage', 'xmomentum', 'ymomentum'])
652        domain.check_integrity()
653
654        domain.set_quantity('stage', [[1,2,3], [5,5,5],
655                                      [0,0,9], [-6, 3, 3]])
656
657        assert num.allclose( domain.quantities['stage'].centroid_values,
658                             [2,5,3,0] )
659
660        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
661                                          [3,3,3], [4, 4, 4]])
662
663        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
664                                          [30,30,30], [40, 40, 40]])
665
666
667        domain.distribute_to_vertices_and_edges()
668
669        #First order extrapolation
670        assert num.allclose( domain.quantities['stage'].vertex_values,
671                             [[ 2.,  2.,  2.],
672                              [ 5.,  5.,  5.],
673                              [ 3.,  3.,  3.],
674                              [ 0.,  0.,  0.]])
675
676        domain.build_tagged_elements_dictionary({'mound':[0,1]})
677        domain.set_region([add_to_verts])
678
679        self.failUnless(domain.test == "Mound",
680                        'set region failed')
681
682
683    def test_track_speeds(self):
684        """
685        get values based on triangle lists.
686        """
687        from mesh_factory import rectangular
688        from shallow_water import Domain
689
690        #Create basic mesh
691        points, vertices, boundary = rectangular(1, 3)
692
693        #Create shallow water domain
694        domain = Domain(points, vertices, boundary)
695        domain.timestepping_statistics(track_speeds=True)
696
697
698
699    def test_region_tags(self):
700        """
701        get values based on triangle lists.
702        """
703        from mesh_factory import rectangular
704        from shallow_water import Domain
705
706        #Create basic mesh
707        points, vertices, boundary = rectangular(1, 3)
708
709        #Create shallow water domain
710        domain = Domain(points, vertices, boundary)
711        domain.build_tagged_elements_dictionary({'bottom':[0,1],
712                                                 'top':[4,5],
713                                                 'all':[0,1,2,3,4,5]})
714
715
716        #Set friction
717        manning = 0.07
718        domain.set_quantity('friction', manning)
719
720        domain.set_region([set_bottom_friction, set_top_friction])
721        #print domain.quantities['friction'].get_values()
722        assert num.allclose(domain.quantities['friction'].get_values(),\
723                            [[ 0.09,  0.09,  0.09],
724                             [ 0.09,  0.09,  0.09],
725                             [ 0.07,  0.07,  0.07],
726                             [ 0.07,  0.07,  0.07],
727                             [ 1.0,  1.0,  1.0],
728                             [ 1.0,  1.0,  1.0]])
729
730        domain.set_region([set_all_friction])
731        #print domain.quantities['friction'].get_values()
732        assert num.allclose(domain.quantities['friction'].get_values(),
733                            [[ 10.09, 10.09, 10.09],
734                             [ 10.09, 10.09, 10.09],
735                             [ 10.07, 10.07, 10.07],
736                             [ 10.07, 10.07, 10.07],
737                             [ 11.0,  11.0,  11.0],
738                             [ 11.0,  11.0,  11.0]])
739
740
741    def test_region_tags2(self):
742        """
743        get values based on triangle lists.
744        """
745        from mesh_factory import rectangular
746        from shallow_water import Domain
747
748        #Create basic mesh
749        points, vertices, boundary = rectangular(1, 3)
750
751        #Create shallow water domain
752        domain = Domain(points, vertices, boundary)
753        domain.build_tagged_elements_dictionary({'bottom':[0,1],
754                                                 'top':[4,5],
755                                                 'all':[0,1,2,3,4,5]})
756
757
758        #Set friction
759        manning = 0.07
760        domain.set_quantity('friction', manning)
761
762        domain.set_region('top', 'friction', 1.0)
763        domain.set_region('bottom', 'friction', 0.09)
764       
765        #print domain.quantities['friction'].get_values()
766        assert num.allclose(domain.quantities['friction'].get_values(),
767                            [[ 0.09,  0.09,  0.09],
768                             [ 0.09,  0.09,  0.09],
769                             [ 0.07,  0.07,  0.07],
770                             [ 0.07,  0.07,  0.07],
771                             [ 1.0,  1.0,  1.0],
772                             [ 1.0,  1.0,  1.0]])
773       
774        domain.set_region([set_bottom_friction, set_top_friction])
775        #print domain.quantities['friction'].get_values()
776        assert num.allclose(domain.quantities['friction'].get_values(),
777                            [[ 0.09,  0.09,  0.09],
778                             [ 0.09,  0.09,  0.09],
779                             [ 0.07,  0.07,  0.07],
780                             [ 0.07,  0.07,  0.07],
781                             [ 1.0,  1.0,  1.0],
782                             [ 1.0,  1.0,  1.0]])
783
784        domain.set_region([set_all_friction])
785        #print domain.quantities['friction'].get_values()
786        assert num.allclose(domain.quantities['friction'].get_values(),
787                            [[ 10.09, 10.09, 10.09],
788                             [ 10.09, 10.09, 10.09],
789                             [ 10.07, 10.07, 10.07],
790                             [ 10.07, 10.07, 10.07],
791                             [ 11.0,  11.0,  11.0],
792                             [ 11.0,  11.0,  11.0]])
793
794#-------------------------------------------------------------
795if __name__ == "__main__":
796    suite = unittest.makeSuite(Test_Domain,'test')
797    runner = unittest.TextTestRunner()
798    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.