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