source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py @ 5571

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

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

File size: 9.6 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5
6from Numeric import allclose, array
7
8
9#from anuga.pyvolution.pmesh2domain import *
10from pmesh2domain import *
11
12from anuga.shallow_water import Domain,\
13     Reflective_boundary, Dirichlet_boundary,\
14     Transmissive_boundary
15
16from anuga.coordinate_transforms.geo_reference import Geo_reference
17from anuga.pmesh.mesh import importMeshFromFile
18
19class Test_pmesh2domain(unittest.TestCase):
20
21    def setUp(self):
22        pass
23
24    def tearDown(self):
25        pass
26
27    def test_pmesh2Domain(self):
28         import os
29         import tempfile
30
31         fileName = tempfile.mktemp(".tsh")
32         file = open(fileName,"w")
33         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
340 0.0 0.0 0.0 0.0 0.01 \n \
351 1.0 0.0 10.0 10.0 0.02  \n \
362 0.0 1.0 0.0 10.0 0.03  \n \
373 0.5 0.25 8.0 12.0 0.04  \n \
38# Vert att title  \n \
39elevation  \n \
40stage  \n \
41friction  \n \
422 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
430 0 3 2 -1  -1  1 dsg\n\
441 0 1 3 -1  0 -1   ole nielsen\n\
454 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
460 1 0 2 \n\
471 0 2 3 \n\
482 2 3 \n\
493 3 1 1 \n\
503 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
510 216.0 -86.0 \n \
521 160.0 -167.0 \n \
532 114.0 -91.0 \n \
543 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
550 0 1 0 \n \
561 1 2 0 \n \
572 2 0 0 \n \
580 # <x> <y> ...Mesh Holes... \n \
590 # <x> <y> <attribute>...Mesh Regions... \n \
600 # <x> <y> <attribute>...Mesh Regions, area... \n\
61#Geo reference \n \
6256 \n \
63140 \n \
64120 \n")
65         file.close()
66
67         tags = {}
68         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
69         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
70         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
71         tags["1"] = b1
72         tags["2"] = b2
73         tags["3"] = b3
74
75         domain = pmesh_to_domain_instance(fileName, Domain)
76         os.remove(fileName)
77         #print "domain.tagged_elements", domain.tagged_elements
78         ## check the quantities
79         #print domain.quantities['elevation'].vertex_values
80         answer = [[0., 8., 0.],
81                   [0., 10., 8.]]
82         assert allclose(domain.quantities['elevation'].vertex_values,
83                        answer)
84
85         #print domain.quantities['stage'].vertex_values
86         answer = [[0., 12., 10.],
87                   [0., 10., 12.]]
88         assert allclose(domain.quantities['stage'].vertex_values,
89                        answer)
90
91         #print domain.quantities['friction'].vertex_values
92         answer = [[0.01, 0.04, 0.03],
93                   [0.01, 0.02, 0.04]]
94         assert allclose(domain.quantities['friction'].vertex_values,
95                        answer)
96
97         #print domain.quantities['friction'].vertex_values
98         assert allclose(domain.tagged_elements['dsg'][0],0)
99         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
100
101         self.failUnless( domain.boundary[(1, 0)]  == '1',
102                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
103         self.failUnless( domain.boundary[(1, 2)]  == '2',
104                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
105         self.failUnless( domain.boundary[(0, 1)]  == '3',
106                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
107         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
108                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
109         #print "domain.boundary",domain.boundary
110         self.failUnless( len(domain.boundary)  == 4,
111                          "test_pmesh2Domain Too many boundaries")
112         #FIXME change to use get_xllcorner
113         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
114         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
115                          "bad geo_referece")
116    #************
117   
118    def test_pmesh2Domain_instance(self):
119         import os
120         import tempfile
121
122         fileName = tempfile.mktemp(".tsh")
123         file = open(fileName,"w")
124         file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
1250 0.0 0.0 0.0 0.0 0.01 \n \
1261 1.0 0.0 10.0 10.0 0.02  \n \
1272 0.0 1.0 0.0 10.0 0.03  \n \
1283 0.5 0.25 8.0 12.0 0.04  \n \
129# Vert att title  \n \
130elevation  \n \
131stage  \n \
132friction  \n \
1332 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
1340 0 3 2 -1  -1  1 dsg\n\
1351 0 1 3 -1  0 -1   ole nielsen\n\
1364 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
1370 1 0 2 \n\
1381 0 2 3 \n\
1392 2 3 \n\
1403 3 1 1 \n\
1413 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
1420 216.0 -86.0 \n \
1431 160.0 -167.0 \n \
1442 114.0 -91.0 \n \
1453 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
1460 0 1 0 \n \
1471 1 2 0 \n \
1482 2 0 0 \n \
1490 # <x> <y> ...Mesh Holes... \n \
1500 # <x> <y> <attribute>...Mesh Regions... \n \
1510 # <x> <y> <attribute>...Mesh Regions, area... \n\
152#Geo reference \n \
15356 \n \
154140 \n \
155120 \n")
156         file.close()
157
158         mesh_instance = importMeshFromFile(fileName)
159       
160         tags = {}
161         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
162         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
163         b3 =  Dirichlet_boundary(conserved_quantities = array([2.0]))
164         tags["1"] = b1
165         tags["2"] = b2
166         tags["3"] = b3
167
168         domain = pmesh_instance_to_domain_instance(mesh_instance, Domain)
169
170         os.remove(fileName)
171         #print "domain.tagged_elements", domain.tagged_elements
172         ## check the quantities
173         #print domain.quantities['elevation'].vertex_values
174         answer = [[0., 8., 0.],
175                   [0., 10., 8.]]
176         assert allclose(domain.quantities['elevation'].vertex_values,
177                        answer)
178
179         #print domain.quantities['stage'].vertex_values
180         answer = [[0., 12., 10.],
181                   [0., 10., 12.]]
182         assert allclose(domain.quantities['stage'].vertex_values,
183                        answer)
184
185         #print domain.quantities['friction'].vertex_values
186         answer = [[0.01, 0.04, 0.03],
187                   [0.01, 0.02, 0.04]]
188         assert allclose(domain.quantities['friction'].vertex_values,
189                        answer)
190
191         #print domain.quantities['friction'].vertex_values
192         assert allclose(domain.tagged_elements['dsg'][0],0)
193         assert allclose(domain.tagged_elements['ole nielsen'][0],1)
194
195         self.failUnless( domain.boundary[(1, 0)]  == '1',
196                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
197         self.failUnless( domain.boundary[(1, 2)]  == '2',
198                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
199         self.failUnless( domain.boundary[(0, 1)]  == '3',
200                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
201         self.failUnless( domain.boundary[(0, 0)]  == 'exterior',
202                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
203         #print "domain.boundary",domain.boundary
204         self.failUnless( len(domain.boundary)  == 4,
205                          "test_pmesh2Domain Too many boundaries")
206         #FIXME change to use get_xllcorner
207         #print "d.geo_reference.xllcorner",domain.geo_reference.xllcorner
208         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
209                          "bad geo_referece")
210         
211    #***********
212    def old_test_tags_to_boundaries (self):
213         meshDict = {}
214         p0 = [0.0, 0.0]
215         p1 = [1.0, 0.0]
216         p2 = [0.0, 1.0]
217         p3 = [0.646446609407, 0.353553390593]
218         meshDict['vertices'] = [p0,p1,p2,p3]
219         meshDict['vertex_attributes'] = [[0.0, 0.0,7.0],[10.0, 0.0,7.0],[0.0, 10.0,7.0],[6.46446609407, 3.53553390593,7.0]]
220         meshDict['triangles'] = [[0,3,2],[0,1,3]]
221         meshDict['triangle_tags'] = [6.6,6.6]
222         meshDict['triangle_neighbors'] = [[-1,-1,1],[-1,0,-1]]
223         meshDict['segments'] = [[1,0],[0,2],[2,3],[3,1]]
224         meshDict['segment_tags'] = [2,3,1,1]
225
226         domain = Domain.pmesh_dictionary_to_domain(meshDict)
227
228         #domain.set_tag_dict(tag_dict)
229         #Boundary tests
230         b1 =  Dirichlet_boundary(conserved_quantities = array([0.0]))
231         b2 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
232         b3 =  Dirichlet_boundary(conserved_quantities = array([1.0]))
233         #test adding a boundary
234         tags = {}
235         tags[1] = b1
236         tags[2] = b2
237         tags[3] = b3
238         domain.set_boundary(tags)
239         inverted_id = Volume.instances[0].neighbours[0]
240         id = -(inverted_id+1)
241         boundary_value = Boundary_value.instances[id]
242         boundary_obj = boundary_value.boundary_object
243
244         self.failUnless( boundary_obj  == b1,
245                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
246
247         inverted_id = Volume.instances[0].neighbours[1]
248         id = -(inverted_id+1)
249         boundary_value = Boundary_value.instances[id]
250         boundary_obj = boundary_value.boundary_object
251         self.failUnless( boundary_obj  == b3,
252                          "test_tags_to_boundaries  failed. Single boundary wasn't added.")
253
254#-------------------------------------------------------------
255if __name__ == "__main__":
256    suite = unittest.makeSuite(Test_pmesh2domain, 'test')
257    runner = unittest.TextTestRunner()
258    runner.run(suite)
259
260
261
262
Note: See TracBrowser for help on using the repository browser.