source: branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py @ 6883

Last change on this file since 6883 was 6689, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk to numpy branch.

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