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

Last change on this file since 7737 was 7484, checked in by steve, 15 years ago

Added edge_midpoint_coordinates together with get_edge_midpoint_coordinates
and get_edge_midpoint_coordinate into the general mesh. This is so that
boundary classes can get access to the location of the edge midpoint.

I will use this to allow Dirichlet_boundary to take a function of
time and space.

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