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

Last change on this file since 6463 was 6463, checked in by rwilson, 15 years ago

Modified General_mesh.get_node(absolute=True) to return a copy if modifying the object node data.

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