source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py @ 8050

Last change on this file since 8050 was 7780, checked in by hudson, 14 years ago

Almost all failing tests fixed.

File size: 15.4 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
11from mesh_factory import rectangular
12from anuga.shallow_water.shallow_water_domain import Domain
13
14import numpy as num
15
16
17class Test_General_Mesh(unittest.TestCase):
18    def setUp(self):
19        pass
20
21    def tearDown(self):
22        pass
23
24
25    def test_get_vertex_coordinates(self):
26        from mesh_factory import rectangular
27
28        #Create basic mesh
29        nodes, triangles, _ = rectangular(1, 3)
30        domain = General_mesh(nodes, triangles)
31
32
33        assert num.allclose(domain.get_nodes(), nodes)
34
35
36        M = domain.number_of_triangles       
37
38        V = domain.get_vertex_coordinates()
39        assert V.shape[0] == 3*M
40
41        for i in range(M):
42            for j in range(3):
43                k = triangles[i,j]  #Index of vertex j in triangle i
44                assert num.allclose(V[3*i+j,:], nodes[k])
45
46    def test_get_vertex_coordinates_with_geo_ref(self):
47        x0 = 314036.58727982
48        y0 = 6224951.2960092
49        geo = Geo_reference(56, x0, y0)
50       
51        a = [0.0, 0.0]
52        b = [0.0, 2.0]
53        c = [2.0, 0.0]
54        d = [0.0, 4.0]
55        e = [2.0, 2.0]
56        f = [4.0, 0.0]
57        nodes = num.array([a, b, c, d, e, f])
58
59        nodes_absolute = geo.get_absolute(nodes)
60       
61        #                        bac,     bce,     ecf,     dbe
62        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
63
64        domain = General_mesh(nodes, triangles, geo_reference=geo)
65
66        verts = domain.get_vertex_coordinates(triangle_id=0)    # bac
67        msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s"
68               % (str(num.array([b,a,c])), str(verts)))
69        self.failUnless(num.allclose(num.array([b,a,c]), verts), msg)
70
71        verts = domain.get_vertex_coordinates(triangle_id=0)       
72        msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s"
73               % (str(num.array([b,a,c])), str(verts)))
74        self.assert_(num.allclose(num.array([b,a,c]), verts), msg)
75
76        verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True)
77        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
78               % (str(num.array([nodes_absolute[1],
79                                 nodes_absolute[0],
80                                 nodes_absolute[2]])),
81                  str(verts)))
82        self.assert_(num.allclose(num.array([nodes_absolute[1],
83                                             nodes_absolute[0],
84                                             nodes_absolute[2]]),
85                                  verts), msg)
86
87        verts = domain.get_vertex_coordinates(triangle_id=0,
88                                              absolute=True)       
89        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
90               % (str(num.array([nodes_absolute[1],
91                                 nodes_absolute[0],
92                                 nodes_absolute[2]])),
93                  str(verts)))
94        self.assert_(num.allclose(num.array([nodes_absolute[1],
95                                             nodes_absolute[0],
96                                             nodes_absolute[2]]),
97                                  verts), msg)
98
99    def test_get_vertex_coordinates_triangle_id(self):
100        """test_get_vertex_coordinates_triangle_id
101        Test that vertices for one triangle can be returned.
102        """
103        from mesh_factory import rectangular
104
105        #Create basic mesh
106        nodes, triangles, _ = rectangular(1, 3)
107        domain = General_mesh(nodes, triangles)
108
109
110        assert num.allclose(domain.get_nodes(), nodes)
111
112
113        M = domain.number_of_triangles       
114
115        for i in range(M):
116            V = domain.get_vertex_coordinates(triangle_id=i)
117            assert V.shape[0] == 3
118
119            for j in range(3):
120                k = triangles[i,j]  #Index of vertex j in triangle i
121                assert num.allclose(V[j,:], nodes[k])
122
123
124    def test_get_edge_midpoint_coordinates(self):
125        from mesh_factory import rectangular
126
127        #Create basic mesh
128        nodes, triangles, _ = rectangular(1, 3)
129        domain = General_mesh(nodes, triangles)
130
131
132        assert num.allclose(domain.get_nodes(), nodes)
133
134
135        M = domain.number_of_triangles       
136
137        E = domain.get_edge_midpoint_coordinates()
138        assert E.shape[0] == 3*M
139
140        for i in range(M):
141            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
142            k1 = triangles[i,1]  #Index of vertex 1 in triangle i
143            k2 = triangles[i,2]  #Index of vertex 2 in triangle i
144
145            assert num.allclose(E[3*i+0,:], 0.5*(nodes[k1]+nodes[k2]) )
146            assert num.allclose(E[3*i+1,:], 0.5*(nodes[k0]+nodes[k2]) )
147            assert num.allclose(E[3*i+2,:], 0.5*(nodes[k1]+nodes[k0]) )
148
149    def test_get_edge_midpoint_coordinates_with_geo_ref(self):
150        x0 = 314036.58727982
151        y0 = 6224951.2960092
152        geo = Geo_reference(56, x0, y0)
153       
154        a = num.array([0.0, 0.0])
155        b = num.array([0.0, 2.0])
156        c = num.array([2.0, 0.0])
157        d = num.array([0.0, 4.0])
158        e = num.array([2.0, 2.0])
159        f = num.array([4.0, 0.0])
160        nodes = num.array([a, b, c, d, e, f])
161
162        nodes_absolute = geo.get_absolute(nodes)
163       
164        #                        bac,     bce,     ecf,     dbe
165        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
166
167        domain = General_mesh(nodes, triangles, geo_reference=geo)
168
169        verts = domain.get_edge_midpoint_coordinates(triangle_id=0)    # bac
170        msg = ("num.array(1/2[a+c,b+c,a+b])=\n%s\nshould be close to 'verts'=\n%s"
171               % (str(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)])), str(verts)))
172        self.failUnless(num.allclose(num.array([0.5*(a+c),0.5*(b+c),0.5*(a+b)]), verts), msg)
173
174
175        verts = domain.get_edge_midpoint_coordinates(triangle_id=0, absolute=True)
176        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
177               % (str(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
178                                     nodes_absolute[1]+nodes_absolute[2],
179                                     nodes_absolute[1]+nodes_absolute[0]])),
180                  str(verts)))
181        self.assert_(num.allclose(0.5*num.array([nodes_absolute[0]+nodes_absolute[2],
182                                                 nodes_absolute[1]+nodes_absolute[2],
183                                                 nodes_absolute[1]+nodes_absolute[0]]),
184                                  verts), msg)
185
186    def test_get_edge_midpoint_coordinates_triangle_id(self):
187        """test_get_vertex_coordinates_triangle_id
188        Test that vertices for one triangle can be returned.
189        """
190        from mesh_factory import rectangular
191
192        #Create basic mesh
193        nodes, triangles, _ = rectangular(1, 3)
194        domain = General_mesh(nodes, triangles)
195
196
197        assert num.allclose(domain.get_nodes(), nodes)
198
199
200        M = domain.number_of_triangles       
201
202        for i in range(M):
203            E = domain.get_edge_midpoint_coordinates(triangle_id=i)
204            assert E.shape[0] == 3
205
206
207            k0 = triangles[i,0]  #Index of vertex 0 in triangle i
208            k1 = triangles[i,1]  #Index of vertex 0 in triangle i
209            k2 = triangles[i,2]  #Index of vertex 0 in triangle i
210
211            assert num.allclose(E[0,:], 0.5*(nodes[k1]+nodes[k2]))
212            assert num.allclose(E[1,:], 0.5*(nodes[k0]+nodes[k2]))
213            assert num.allclose(E[2,:], 0.5*(nodes[k1]+nodes[k0]))
214
215            E0 = domain.get_edge_midpoint_coordinate(i, 0 )
216            E1 = domain.get_edge_midpoint_coordinate(i, 1 )
217            E2 = domain.get_edge_midpoint_coordinate(i, 2 )
218
219            assert num.allclose(E0, 0.5*(nodes[k1]+nodes[k2]))
220            assert num.allclose(E1, 0.5*(nodes[k0]+nodes[k2]))
221            assert num.allclose(E2, 0.5*(nodes[k1]+nodes[k0]))
222           
223
224
225       
226       
227
228    def test_get_vertex_values(self):
229        """Get connectivity based on triangle lists.
230        """
231        from mesh_factory import rectangular
232
233        #Create basic mesh
234        nodes, triangles, _ = rectangular(1, 3)
235        domain = General_mesh(nodes, triangles)
236
237        msg = ("domain.get_triangles()=\n%s\nshould be the same as "
238               "'triangles'=\n%s"
239               % (str(domain.get_triangles()), str(triangles)))
240        assert num.allclose(domain.get_triangles(), triangles), msg
241        msg = ("domain.get_triangles([0,4])=\n%s\nshould be the same as "
242               "'[triangles[0], triangles[4]]' which is\n%s"
243               % (str(domain.get_triangles([0,4])),
244                  str([triangles[0], triangles[4]])))
245        assert num.allclose(domain.get_triangles([0,4]),
246                            [triangles[0], triangles[4]]), msg
247       
248
249    def test_vertex_value_indices(self):
250        """Check that structures are correct.
251        """
252        from mesh_factory import rectangular
253
254        a = [0.0, 0.0]
255        b = [0.0, 2.0]
256        c = [2.0, 0.0]
257        d = [0.0, 4.0]
258        e = [2.0, 2.0]
259        f = [4.0, 0.0]
260
261        nodes = num.array([a, b, c, d, e, f])
262        #bac, bce, ecf, dbe, daf, dae
263        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
264
265        domain1 = General_mesh(nodes, triangles)
266       
267        #Create larger mesh
268        nodes, triangles, _ = rectangular(3, 6)
269        domain2 = General_mesh(nodes, triangles)
270
271        # Test both meshes
272        for domain in [domain1, domain2]:
273            assert sum(domain.number_of_triangles_per_node) ==\
274                   len(domain.vertex_value_indices)
275
276            # Check number of triangles per node
277            count = [0]*domain.number_of_nodes
278            for triangle in domain.triangles:
279                for i in triangle:
280                    count[i] += 1
281
282            #print count
283            #
284            assert num.allclose(count, domain.number_of_triangles_per_node)
285           
286            # Check indices
287            current_node = 0
288            k = 0 # Track triangles touching on node
289            for index in domain.vertex_value_indices:
290                k += 1
291               
292                triangle = index / 3
293                vertex = index % 3
294
295                assert domain.triangles[triangle, vertex] == current_node
296
297                if domain.number_of_triangles_per_node[current_node] == k:
298                    # Move on to next node
299                    k = 0
300                    current_node += 1
301               
302
303    def test_get_triangles_and_vertices_per_node(self):
304        """test_get_triangles_and_vertices_per_node -
305
306        Test that tuples of triangle, vertex can be extracted
307        from inverted triangles structure
308
309        """
310
311        a = [0.0, 0.0]
312        b = [0.0, 2.0]
313        c = [2.0, 0.0]
314        d = [0.0, 4.0]
315        e = [2.0, 2.0]
316        f = [4.0, 0.0]
317
318        nodes = num.array([a, b, c, d, e, f])
319        #bac, bce, ecf, dbe, daf, dae
320        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
321
322        domain = General_mesh(nodes, triangles)
323
324        # One node
325        L = domain.get_triangles_and_vertices_per_node(node=2)
326        assert num.allclose(L[0], [0, 2])
327        assert num.allclose(L[1], [1, 1])
328        assert num.allclose(L[2], [2, 1])
329
330        # All nodes
331        ALL = domain.get_triangles_and_vertices_per_node()
332        assert len(ALL) == 6
333        for i, Lref in enumerate(ALL):
334            L = domain.get_triangles_and_vertices_per_node(node=i)
335            assert num.allclose(L, Lref)
336           
337
338       
339
340       
341
342
343    def test_areas(self):
344        #Create basic mesh
345        points, vertices, boundary = rectangular(1, 3)
346        domain = General_mesh(points, vertices)       
347
348        assert domain.get_area() == 1.0
349
350
351    def test_get_unique_vertex_values(self):
352        """
353        get unique_vertex based on triangle lists.
354        """
355
356        #Create basic mesh
357        points, vertices, boundary = rectangular(1, 3)
358        domain = General_mesh(points, vertices)               
359
360        assert  domain.get_unique_vertices() == [0,1,2,3,4,5,6,7]
361        unique_vertices = domain.get_unique_vertices([0,1,4])
362        unique_vertices.sort()
363        assert unique_vertices == [0,1,2,4,5,6,7]
364
365        unique_vertices = domain.get_unique_vertices([0,4])
366        unique_vertices.sort()
367        assert unique_vertices == [0,2,4,5,6,7]
368
369    def test_get_node(self):
370        """test_get_triangles_and_vertices_per_node -
371
372        Test that tuples of triangle, vertex can be extracted
373        from inverted triangles structure
374
375        """
376       
377        x0 = 314036.58727982
378        y0 = 6224951.2960092
379        geo = Geo_reference(56, x0, y0)
380       
381        a = [0.0, 0.0]
382        b = [0.0, 2.0]
383        c = [2.0, 0.0]
384        d = [0.0, 4.0]
385        e = [2.0, 2.0]
386        f = [4.0, 0.0]
387
388        nodes = num.array([a, b, c, d, e, f])
389
390        nodes_absolute = geo.get_absolute(nodes)
391       
392        #                        bac,     bce,     ecf,     dbe
393        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
394
395        domain = General_mesh(nodes, triangles, geo_reference = geo)
396        node = domain.get_node(2)       
397        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
398        self.failUnless(num.alltrue(c == node), msg)
399
400        # repeat get_node(), see if result same
401        node = domain.get_node(2)       
402        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
403        self.failUnless(num.alltrue(c == node), msg)
404       
405        node = domain.get_node(2, absolute=True)     
406        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
407               % (str(nodes_absolute[2]), str(node)))
408        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
409       
410        # repeat get_node(2, absolute=True), see if result same
411        node = domain.get_node(2, absolute=True)     
412        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
413               % (str(nodes_absolute[2]), str(node)))
414        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
415       
416
417    def test_assert_index_in_nodes(self):
418        """test_assert_index_in_nodes -
419
420        Test that node indices in triangles are within nodes array.
421
422        """
423       
424        x0 = 314036.58727982
425        y0 = 6224951.2960092
426        geo = Geo_reference(56, x0, y0)
427       
428        a = [0.0, 0.0]
429        b = [0.0, 2.0]
430        c = [2.0, 0.0]
431        d = [0.0, 4.0]
432        e = [2.0, 2.0]
433        f = [4.0, 0.0]
434
435        nodes = num.array([a, b, c, d, e, f])
436
437        nodes_absolute = geo.get_absolute(nodes)
438       
439        # max index is 5, use 5, expect success
440        triangles = num.array([[1,5,2], [1,2,4], [4,2,5], [3,1,4]])
441        General_mesh(nodes, triangles, geo_reference=geo)
442       
443        # max index is 5, use 6, expect assert failure
444        triangles = num.array([[1,6,2], [1,2,4], [4,2,5], [3,1,4]])
445        self.failUnlessRaises(AssertionError, General_mesh,
446                              nodes, triangles, geo_reference=geo)
447       
448        # max index is 5, use 10, expect assert failure
449        triangles = num.array([[1,10,2], [1,2,4], [4,2,5], [3,1,4]])
450        self.failUnlessRaises(AssertionError, General_mesh,
451                              nodes, triangles, geo_reference=geo)
452
453################################################################################
454
455if __name__ == "__main__":
456    #suite = unittest.makeSuite(Test_General_Mesh, 'test')
457    suite = unittest.makeSuite(Test_General_Mesh, 'test_get_edge_midpoint_coordinates')     
458    runner = unittest.TextTestRunner()
459    runner.run(suite)
460
Note: See TracBrowser for help on using the repository browser.