source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py @ 7928

Last change on this file since 7928 was 7780, checked in by hudson, 14 years ago

Almost all failing tests fixed.

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