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

Last change on this file since 5571 was 4839, checked in by duncan, 16 years ago

speeding up slow functions in general mesh. Fix for ticket 201. Also checking in GA validation test for profiling fitting

File size: 8.8 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from math import sqrt, pi
6
7
8from anuga.config import epsilon
9from Numeric import allclose, array, ones, Float
10from general_mesh import General_mesh
11from anuga.coordinate_transforms.geo_reference import Geo_reference
12
13
14
15class Test_General_Mesh(unittest.TestCase):
16    def setUp(self):
17        pass
18
19    def tearDown(self):
20        pass
21
22
23    def test_get_vertex_coordinates(self):
24        from mesh_factory import rectangular
25        from Numeric import zeros, Float
26
27        #Create basic mesh
28        nodes, triangles, _ = rectangular(1, 3)
29        domain = General_mesh(nodes, triangles)
30
31
32        assert allclose(domain.get_nodes(), nodes)
33
34
35        M = domain.number_of_triangles       
36
37        V = domain.get_vertex_coordinates()
38        assert V.shape[0] == 3*M
39
40        for i in range(M):
41            for j in range(3):
42                k = triangles[i,j]  #Index of vertex j in triangle i
43                assert allclose(V[3*i+j,:], nodes[k])
44
45    def test_get_vertex_coordinates_with_geo_ref(self):
46        x0 = 314036.58727982
47        y0 = 6224951.2960092
48        geo = Geo_reference(56, x0, y0)
49       
50        a = [0.0, 0.0]
51        b = [0.0, 2.0]
52        c = [2.0, 0.0]
53        d = [0.0, 4.0]
54        e = [2.0, 2.0]
55        f = [4.0, 0.0]
56
57        nodes = array([a, b, c, d, e, f])
58
59        nodes_absolute = geo.get_absolute(nodes)
60       
61        #bac, bce, ecf, dbe, daf, dae
62        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
63
64        domain = General_mesh(nodes, triangles,
65                       geo_reference = geo)
66        verts = domain.get_vertex_coordinates(triangle_id=0)       
67        self.assert_(allclose(array([b,a,c]), verts)) 
68        verts = domain.get_vertex_coordinates(triangle_id=0)       
69        self.assert_(allclose(array([b,a,c]), verts))
70        verts = domain.get_vertex_coordinates(triangle_id=0,
71                                              absolute=True)       
72        self.assert_(allclose(array([nodes_absolute[1],
73                                     nodes_absolute[0],
74                                     nodes_absolute[2]]), verts))
75        verts = domain.get_vertex_coordinates(triangle_id=0,
76                                              absolute=True)       
77        self.assert_(allclose(array([nodes_absolute[1],
78                                     nodes_absolute[0],
79                                     nodes_absolute[2]]), verts))
80       
81       
82
83    def test_get_vertex_coordinates_triangle_id(self):
84        """test_get_vertex_coordinates_triangle_id
85        Test that vertices for one triangle can be returned.
86        """
87        from mesh_factory import rectangular
88        from Numeric import zeros, Float
89
90        #Create basic mesh
91        nodes, triangles, _ = rectangular(1, 3)
92        domain = General_mesh(nodes, triangles)
93
94
95        assert allclose(domain.get_nodes(), nodes)
96
97
98        M = domain.number_of_triangles       
99
100        for i in range(M):
101            V = domain.get_vertex_coordinates(triangle_id=i)
102            assert V.shape[0] == 3
103
104            for j in range(3):
105                k = triangles[i,j]  #Index of vertex j in triangle i
106                assert allclose(V[j,:], nodes[k])
107
108
109       
110       
111
112    def test_get_vertex_values(self):
113        """Get connectivity based on triangle lists.
114        """
115        from mesh_factory import rectangular
116        from Numeric import zeros, Float
117
118        #Create basic mesh
119        nodes, triangles, _ = rectangular(1, 3)
120        domain = General_mesh(nodes, triangles)
121
122        value = [7]
123        assert allclose(domain.get_triangles(), triangles)
124        assert allclose(domain.get_triangles([0,4]),
125                        [triangles[0], triangles[4]])
126       
127
128    def test_vertex_value_indices(self):
129        """Check that structures are correct.
130        """
131        from mesh_factory import rectangular
132        from Numeric import zeros, Float, array
133
134        a = [0.0, 0.0]
135        b = [0.0, 2.0]
136        c = [2.0, 0.0]
137        d = [0.0, 4.0]
138        e = [2.0, 2.0]
139        f = [4.0, 0.0]
140
141        nodes = array([a, b, c, d, e, f])
142        #bac, bce, ecf, dbe, daf, dae
143        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
144
145        domain1 = General_mesh(nodes, triangles)
146       
147        #Create larger mesh
148        nodes, triangles, _ = rectangular(3, 6)
149        domain2 = General_mesh(nodes, triangles)
150
151        # Test both meshes
152        for domain in [domain1, domain2]:
153            assert sum(domain.number_of_triangles_per_node) ==\
154                   len(domain.vertex_value_indices)
155
156            # Check number of triangles per node
157            count = [0]*domain.number_of_nodes
158            for triangle in domain.triangles:
159                for i in triangle:
160                    count[i] += 1
161
162            #print count
163            #
164            assert allclose(count, domain.number_of_triangles_per_node)
165           
166            # Check indices
167            current_node = 0
168            k = 0 # Track triangles touching on node
169            for index in domain.vertex_value_indices:
170                k += 1
171               
172                triangle = index / 3
173                vertex = index % 3
174
175                assert domain.triangles[triangle, vertex] == current_node
176
177                if domain.number_of_triangles_per_node[current_node] == k:
178                    # Move on to next node
179                    k = 0
180                    current_node += 1
181               
182
183    def test_get_triangles_and_vertices_per_node(self):
184        """test_get_triangles_and_vertices_per_node -
185
186        Test that tuples of triangle, vertex can be extracted
187        from inverted triangles structure
188
189        """
190
191        a = [0.0, 0.0]
192        b = [0.0, 2.0]
193        c = [2.0, 0.0]
194        d = [0.0, 4.0]
195        e = [2.0, 2.0]
196        f = [4.0, 0.0]
197
198        nodes = array([a, b, c, d, e, f])
199        #bac, bce, ecf, dbe, daf, dae
200        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
201
202        domain = General_mesh(nodes, triangles)
203
204        # One node
205        L = domain.get_triangles_and_vertices_per_node(node=2)
206        assert allclose(L[0], [0, 2])
207        assert allclose(L[1], [1, 1])
208        assert allclose(L[2], [2, 1])
209
210        # All nodes
211        ALL = domain.get_triangles_and_vertices_per_node()
212        assert len(ALL) == 6
213        for i, Lref in enumerate(ALL):
214            L = domain.get_triangles_and_vertices_per_node(node=i)
215            assert allclose(L, Lref)
216           
217
218       
219
220       
221
222
223    def test_areas(self):
224        from mesh_factory import rectangular
225        from shallow_water import Domain
226        from Numeric import zeros, Float
227
228        #Create basic mesh
229        points, vertices, boundary = rectangular(1, 3)
230        domain = General_mesh(points, vertices)       
231
232        assert domain.get_area() == 1.0
233
234
235    def test_get_unique_vertex_values(self):
236        """
237        get unique_vertex based on triangle lists.
238        """
239        from mesh_factory import rectangular
240        from shallow_water import Domain
241        from Numeric import zeros, Float
242
243        #Create basic mesh
244        points, vertices, boundary = rectangular(1, 3)
245        domain = General_mesh(points, vertices)               
246
247        assert  domain.get_unique_vertices() == [0,1,2,3,4,5,6,7]
248        unique_vertices = domain.get_unique_vertices([0,1,4])
249        unique_vertices.sort()
250        assert unique_vertices == [0,1,2,4,5,6,7]
251
252        unique_vertices = domain.get_unique_vertices([0,4])
253        unique_vertices.sort()
254        assert unique_vertices == [0,2,4,5,6,7]
255
256    def test_get_node(self):
257        """test_get_triangles_and_vertices_per_node -
258
259        Test that tuples of triangle, vertex can be extracted
260        from inverted triangles structure
261
262        """
263       
264        x0 = 314036.58727982
265        y0 = 6224951.2960092
266        geo = Geo_reference(56, x0, y0)
267       
268        a = [0.0, 0.0]
269        b = [0.0, 2.0]
270        c = [2.0, 0.0]
271        d = [0.0, 4.0]
272        e = [2.0, 2.0]
273        f = [4.0, 0.0]
274
275        nodes = array([a, b, c, d, e, f])
276
277        nodes_absolute = geo.get_absolute(nodes)
278       
279        #bac, bce, ecf, dbe, daf, dae
280        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
281
282        domain = General_mesh(nodes, triangles,
283                       geo_reference = geo)
284        node = domain.get_node(2)       
285        self.assertEqual(c, node)
286       
287        node = domain.get_node(2, absolute=True)     
288        self.assertEqual(nodes_absolute[2], node)
289       
290        node = domain.get_node(2, absolute=True)     
291        self.assertEqual(nodes_absolute[2], node)
292       
293
294
295#-------------------------------------------------------------
296if __name__ == "__main__":
297    suite = unittest.makeSuite(Test_General_Mesh,'test') 
298    #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')   
299    runner = unittest.TextTestRunner()
300    runner.run(suite)
301
Note: See TracBrowser for help on using the repository browser.