source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py @ 7276

Last change on this file since 7276 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

File size: 29.6 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_setting_timestepping_method
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.set_default_order(1)
545        domain.check_integrity()
546
547
548        domain.set_quantity('stage', [[1,2,3], [5,5,5],
549                                      [0,0,9], [-6, 3, 3]])
550
551        assert num.allclose( domain.quantities['stage'].centroid_values,
552                             [2,5,3,0] )
553
554        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
555                                          [3,3,3], [4, 4, 4]])
556
557        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
558                                          [30,30,30], [40, 40, 40]])
559
560
561        domain.distribute_to_vertices_and_edges()
562
563        #First order extrapolation
564        assert num.allclose( domain.quantities['stage'].vertex_values,
565                             [[ 2.,  2.,  2.],
566                              [ 5.,  5.,  5.],
567                              [ 3.,  3.,  3.],
568                              [ 0.,  0.,  0.]])
569
570
571
572
573    def test_update_conserved_quantities(self):
574        a = [0.0, 0.0]
575        b = [0.0, 2.0]
576        c = [2.0,0.0]
577        d = [0.0, 4.0]
578        e = [2.0, 2.0]
579        f = [4.0,0.0]
580
581        points = [a, b, c, d, e, f]
582        #bac, bce, ecf, dbe
583        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
584        boundary = { (0, 0): 'Third',
585                     (0, 2): 'First',
586                     (2, 0): 'Second',
587                     (2, 1): 'Second',
588                     (3, 1): 'Second',
589                     (3, 2): 'Third'}
590
591
592        domain = Domain(points, vertices, boundary,
593                        conserved_quantities =\
594                        ['stage', 'xmomentum', 'ymomentum'])
595        domain.check_integrity()
596
597
598        domain.set_quantity('stage', [1,2,3,4], location='centroids')
599        domain.set_quantity('xmomentum', [1,2,3,4], location='centroids')
600        domain.set_quantity('ymomentum', [1,2,3,4], location='centroids')
601
602
603        #Assign some values to update vectors
604        #Set explicit_update
605
606        for name in domain.conserved_quantities:
607            domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.])
608            domain.quantities[name].semi_implicit_update = num.array([1.,1.,1.,1.])
609
610
611        #Update with given timestep (assuming no other forcing terms)
612        domain.timestep = 0.1
613        domain.update_conserved_quantities()
614
615        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
616        denom = num.ones(4, num.float) - domain.timestep*sem
617
618#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
619#        x /= denom
620
621        x = num.array([1., 2., 3., 4.])
622        x /= denom
623        x += domain.timestep*num.array( [4,3,2,1] )
624
625        for name in domain.conserved_quantities:
626            assert num.allclose(domain.quantities[name].centroid_values, x)
627
628
629    def test_set_region(self):
630        """Set quantities for sub region
631        """
632
633        a = [0.0, 0.0]
634        b = [0.0, 2.0]
635        c = [2.0,0.0]
636        d = [0.0, 4.0]
637        e = [2.0, 2.0]
638        f = [4.0,0.0]
639
640        points = [a, b, c, d, e, f]
641        #bac, bce, ecf, dbe
642        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
643        boundary = { (0, 0): 'Third',
644                     (0, 2): 'First',
645                     (2, 0): 'Second',
646                     (2, 1): 'Second',
647                     (3, 1): 'Second',
648                     (3, 2): 'Third'}
649
650        domain = Domain(points, vertices, boundary,
651                        conserved_quantities =\
652                        ['stage', 'xmomentum', 'ymomentum'])
653        domain.set_default_order(1)                       
654        domain.check_integrity()
655
656        domain.set_quantity('stage', [[1,2,3], [5,5,5],
657                                      [0,0,9], [-6, 3, 3]])
658
659        assert num.allclose( domain.quantities['stage'].centroid_values,
660                             [2,5,3,0] )
661
662        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
663                                          [3,3,3], [4, 4, 4]])
664
665        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
666                                          [30,30,30], [40, 40, 40]])
667
668
669        domain.distribute_to_vertices_and_edges()
670
671        #First order extrapolation
672        assert num.allclose( domain.quantities['stage'].vertex_values,
673                             [[ 2.,  2.,  2.],
674                              [ 5.,  5.,  5.],
675                              [ 3.,  3.,  3.],
676                              [ 0.,  0.,  0.]])
677
678        domain.build_tagged_elements_dictionary({'mound':[0,1]})
679        domain.set_region([add_to_verts])
680
681        self.failUnless(domain.test == "Mound",
682                        'set region failed')
683
684
685    def test_track_speeds(self):
686        """
687        get values based on triangle lists.
688        """
689        from mesh_factory import rectangular
690        from shallow_water import Domain
691
692        #Create basic mesh
693        points, vertices, boundary = rectangular(1, 3)
694
695        #Create shallow water domain
696        domain = Domain(points, vertices, boundary)
697        domain.timestepping_statistics(track_speeds=True)
698
699
700
701    def test_region_tags(self):
702        """
703        get values based on triangle lists.
704        """
705        from mesh_factory import rectangular
706        from shallow_water import Domain
707
708        #Create basic mesh
709        points, vertices, boundary = rectangular(1, 3)
710
711        #Create shallow water domain
712        domain = Domain(points, vertices, boundary)
713        domain.build_tagged_elements_dictionary({'bottom':[0,1],
714                                                 'top':[4,5],
715                                                 'all':[0,1,2,3,4,5]})
716
717
718        #Set friction
719        manning = 0.07
720        domain.set_quantity('friction', manning)
721
722        domain.set_region([set_bottom_friction, set_top_friction])
723        #print domain.quantities['friction'].get_values()
724        assert num.allclose(domain.quantities['friction'].get_values(),\
725                            [[ 0.09,  0.09,  0.09],
726                             [ 0.09,  0.09,  0.09],
727                             [ 0.07,  0.07,  0.07],
728                             [ 0.07,  0.07,  0.07],
729                             [ 1.0,  1.0,  1.0],
730                             [ 1.0,  1.0,  1.0]])
731
732        domain.set_region([set_all_friction])
733        #print domain.quantities['friction'].get_values()
734        assert num.allclose(domain.quantities['friction'].get_values(),
735                            [[ 10.09, 10.09, 10.09],
736                             [ 10.09, 10.09, 10.09],
737                             [ 10.07, 10.07, 10.07],
738                             [ 10.07, 10.07, 10.07],
739                             [ 11.0,  11.0,  11.0],
740                             [ 11.0,  11.0,  11.0]])
741
742
743    def test_region_tags2(self):
744        """
745        get values based on triangle lists.
746        """
747        from mesh_factory import rectangular
748        from shallow_water import Domain
749
750        #Create basic mesh
751        points, vertices, boundary = rectangular(1, 3)
752
753        #Create shallow water domain
754        domain = Domain(points, vertices, boundary)
755        domain.build_tagged_elements_dictionary({'bottom':[0,1],
756                                                 'top':[4,5],
757                                                 'all':[0,1,2,3,4,5]})
758
759
760        #Set friction
761        manning = 0.07
762        domain.set_quantity('friction', manning)
763
764        domain.set_region('top', 'friction', 1.0)
765        domain.set_region('bottom', 'friction', 0.09)
766       
767        #print domain.quantities['friction'].get_values()
768        msg = ("domain.quantities['friction'].get_values()=\n%s\n"
769               'should equal\n'
770               '[[ 0.09,  0.09,  0.09],\n'
771               ' [ 0.09,  0.09,  0.09],\n'
772               ' [ 0.07,  0.07,  0.07],\n'
773               ' [ 0.07,  0.07,  0.07],\n'
774               ' [ 1.0,  1.0,  1.0],\n'
775               ' [ 1.0,  1.0,  1.0]]'
776               % str(domain.quantities['friction'].get_values()))
777        assert num.allclose(domain.quantities['friction'].get_values(),
778                            [[ 0.09,  0.09,  0.09],
779                             [ 0.09,  0.09,  0.09],
780                             [ 0.07,  0.07,  0.07],
781                             [ 0.07,  0.07,  0.07],
782                             [ 1.0,  1.0,  1.0],
783                             [ 1.0,  1.0,  1.0]]), msg
784       
785        domain.set_region([set_bottom_friction, set_top_friction])
786        #print domain.quantities['friction'].get_values()
787        assert num.allclose(domain.quantities['friction'].get_values(),
788                            [[ 0.09,  0.09,  0.09],
789                             [ 0.09,  0.09,  0.09],
790                             [ 0.07,  0.07,  0.07],
791                             [ 0.07,  0.07,  0.07],
792                             [ 1.0,  1.0,  1.0],
793                             [ 1.0,  1.0,  1.0]])
794
795        domain.set_region([set_all_friction])
796        #print domain.quantities['friction'].get_values()
797        assert num.allclose(domain.quantities['friction'].get_values(),
798                            [[ 10.09, 10.09, 10.09],
799                             [ 10.09, 10.09, 10.09],
800                             [ 10.07, 10.07, 10.07],
801                             [ 10.07, 10.07, 10.07],
802                             [ 11.0,  11.0,  11.0],
803                             [ 11.0,  11.0,  11.0]])
804                             
805    def test_rectangular_periodic_and_ghosts(self):
806
807        from mesh_factory import rectangular_periodic
808       
809
810        M=5
811        N=2
812        points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N)
813
814        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
815        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
816
817        #print 'HERE'
818       
819        conserved_quantities = ['quant1', 'quant2']
820        domain = Domain(points, vertices, boundary, conserved_quantities,
821                        full_send_dict=full_send_dict,
822                        ghost_recv_dict=ghost_recv_dict)
823
824
825
826
827        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
828        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
829
830        def xylocation(x,y):
831            return 15*x + 9*y
832
833       
834        domain.set_quantity('quant1',xylocation,location='centroids')
835        domain.set_quantity('quant2',xylocation,location='centroids')
836
837
838        assert num.allclose(domain.quantities['quant1'].centroid_values,
839                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
840                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
841                               18.5,  19.,   23.,   23.5])
842
843
844
845        assert num.allclose(domain.quantities['quant2'].centroid_values,
846                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
847                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
848                               18.5,  19.,   23.,   23.5])
849
850        domain.update_ghosts()
851
852
853        assert num.allclose(domain.quantities['quant1'].centroid_values,
854                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
855                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
856                                3.5,   4.,    8.,    8.5])
857
858
859
860        assert num.allclose(domain.quantities['quant2'].centroid_values,
861                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
862                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
863                                3.5,   4.,    8.,    8.5])
864
865       
866        assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
867                                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])
868
869        #Test that points are arranged in a counter clock wise order
870        domain.check_integrity()
871
872
873    def test_that_mesh_methods_exist(self):
874        """test_that_mesh_methods_exist
875       
876        Test that relavent mesh methods are made available in
877        domain through composition
878        """
879        from mesh_factory import rectangular
880        from shallow_water import Domain
881
882        # Create basic mesh
883        points, vertices, boundary = rectangular(1, 3)
884
885        # Create shallow water domain
886        domain = Domain(points, vertices, boundary)                             
887       
888       
889        domain.get_centroid_coordinates()
890        domain.get_radii()
891        domain.get_areas()
892        domain.get_area()
893        domain.get_vertex_coordinates()
894        domain.get_triangles()
895        domain.get_nodes()
896        domain.get_number_of_nodes()
897        domain.get_normal(0,0)
898        domain.get_triangle_containing_point([0.4,0.5])
899        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
900        domain.get_disconnected_triangles()
901        domain.get_boundary_tags()
902        domain.get_boundary_polygon()
903        #domain.get_number_of_triangles_per_node()
904        domain.get_triangles_and_vertices_per_node()
905        domain.get_interpolation_object()
906        domain.get_tagged_elements()
907        domain.get_lone_vertices()
908        domain.get_unique_vertices()
909        g = domain.get_georeference()
910        domain.set_georeference(g)
911        domain.build_tagged_elements_dictionary()
912        domain.statistics()
913        domain.get_extent()
914
915       
916       
917
918#-------------------------------------------------------------
919
920if __name__ == "__main__":
921    suite = unittest.makeSuite(Test_Domain,'test')
922    runner = unittest.TextTestRunner()
923    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.