source: inundation/ga/storm_surge/pyvolution/test_mesh.py @ 333

Last change on this file since 333 was 322, checked in by duncan, 20 years ago

mesh inherits from general_mesh

File size: 17.1 KB
Line 
1#!/usr/bin/env python
2
3 #FIXME: Seperate the tests for mesh and general_mesh
4
5import unittest
6from math import sqrt
7
8from mesh import *
9from config import epsilon
10from Numeric import allclose, array
11
12def distance(x, y):
13    return sqrt( sum( (array(x)-array(y))**2 ))         
14       
15class TestCase(unittest.TestCase):
16    def setUp(self):
17        pass
18       
19    def tearDown(self):
20        pass
21
22    def test_triangle_inputs(self):
23        points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
24        vertices = [0,1,2] #Wrong
25
26        try:
27            mesh = Mesh(points, vertices)
28        except:
29            pass
30        else:
31            msg = 'Should have raised exception'
32            raise msg
33       
34       
35    def test_basic_triangle(self):
36
37        a = [0.0, 0.0]
38        b = [4.0, 0.0]
39        c = [0.0, 3.0]
40       
41        points = [a, b, c]
42        vertices = [[0,1,2]]
43        mesh = Mesh(points, vertices) 
44
45        #Centroid
46        centroid = mesh.centroid_coordinates[0]
47        assert centroid[0] == 4.0/3
48        assert centroid[1] == 1.0
49
50        #Area
51        assert mesh.areas[0] == 6.0,\
52               'Area was %f, should have been 6.0' %mesh.areas[0]
53
54        #Normals
55        normals = mesh.get_normals()
56        assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
57        assert allclose(normals[0, 2:4], [-1.0, 0.0])
58        assert allclose(normals[0, 4:6], [0.0, -1.0])       
59
60        assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
61        assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
62        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])               
63       
64        #Edge lengths
65        assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
66
67
68        #Vertex coordinates
69        V = mesh.get_vertex_coordinates()
70        assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
71
72        V0 = mesh.get_vertex_coordinate(0, 0)
73        assert allclose(V0, [0.0, 0.0])
74
75        V1 = mesh.get_vertex_coordinate(0, 1)
76        assert allclose(V1, [4.0, 0.0])
77
78        V2 = mesh.get_vertex_coordinate(0, 2)
79        assert allclose(V2, [0.0, 3.0])               
80
81
82        #General tests:
83
84        #Test that points are arranged in a counter clock wise order etc
85        mesh.check_integrity()
86
87       
88        #Test that the centroid is located 2/3 of the way
89        #from each vertex to the midpoint of the opposite side
90
91        V = mesh.get_vertex_coordinates()
92       
93        x0 = V[0,0]
94        y0 = V[0,1]
95        x1 = V[0,2]
96        y1 = V[0,3]
97        x2 = V[0,4]
98        y2 = V[0,5]
99       
100        m0 = [(x1 + x2)/2, (y1 + y2)/2]
101        m1 = [(x0 + x2)/2, (y0 + y2)/2]
102        m2 = [(x1 + x0)/2, (y1 + y0)/2]
103
104        d0 = distance(centroid, [x0, y0])
105        d1 = distance(m0, [x0, y0])
106        assert d0 == 2*d1/3
107        #
108        d0 = distance(centroid, [x1, y1])
109        d1 = distance(m1, [x1, y1])       
110        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
111
112        d0 = distance(centroid, [x2, y2])
113        d1 = distance(m2, [x2, y2])
114        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)       
115       
116        #Radius
117        d0 = distance(centroid, m0)
118        assert d0 == 5.0/6
119       
120        d1 = distance(centroid, m1)
121        assert d1 == sqrt(73.0/36)       
122
123        d2 = distance(centroid, m2)       
124        assert d2 == sqrt(13.0/9)       
125       
126        assert mesh.radii[0] == min(d0, d1, d2)
127        assert mesh.radii[0] == 5.0/6
128
129
130
131        #Let x be the centroid of triangle abc.
132        #Test that areas of the three triangles axc, cxb, and bxa are equal.
133
134        points = [a, b, c, centroid]
135        vertices = [[0,3,2], [2,3,1], [1,3,0]]
136        new_mesh = Mesh(points, vertices)
137       
138        assert new_mesh.areas[0] == new_mesh.areas[1]
139        assert new_mesh.areas[1] == new_mesh.areas[2]       
140        assert new_mesh.areas[1] == new_mesh.areas[2]       
141
142        assert new_mesh.areas[1] == mesh.areas[0]/3       
143
144
145
146    def test_general_triangle(self):
147        a = [2.0, 1.0]
148        b = [6.0, 2.0]
149        c = [1.0, 3.0]
150
151        points = [a, b, c]
152        vertices = [[0,1,2]]
153
154        mesh = Mesh(points, vertices) 
155        centroid = mesh.centroid_coordinates[0]
156
157       
158        #Test that the centroid is located 2/3 of the way
159        #from each vertex to the midpoint of the opposite side
160
161        V = mesh.get_vertex_coordinates()
162       
163        x0 = V[0,0]
164        y0 = V[0,1]
165        x1 = V[0,2]
166        y1 = V[0,3]
167        x2 = V[0,4]
168        y2 = V[0,5]
169       
170        m0 = [(x1 + x2)/2, (y1 + y2)/2]
171        m1 = [(x0 + x2)/2, (y0 + y2)/2]
172        m2 = [(x1 + x0)/2, (y1 + y0)/2]
173
174        d0 = distance(centroid, [x0, y0])
175        d1 = distance(m0, [x0, y0])
176        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
177        #
178        d0 = distance(centroid, [x1, y1])
179        d1 = distance(m1, [x1, y1])       
180        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
181
182        d0 = distance(centroid, [x2, y2])
183        d1 = distance(m2, [x2, y2])
184        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)       
185       
186        #Radius
187        d0 = distance(centroid, m0)
188        d1 = distance(centroid, m1)
189        d2 = distance(centroid, m2)       
190        assert mesh.radii[0] == min(d0, d1, d2)
191
192
193
194        #Let x be the centroid of triangle abc.
195        #Test that areas of the three triangles axc, cxb, and bxa are equal.
196
197        points = [a, b, c, centroid]
198        vertices = [[0,3,2], [2,3,1], [1,3,0]]
199        new_mesh = Mesh(points, vertices)
200       
201        assert new_mesh.areas[0] == new_mesh.areas[1]
202        assert new_mesh.areas[1] == new_mesh.areas[2]       
203        assert new_mesh.areas[1] == new_mesh.areas[2]       
204
205        assert new_mesh.areas[1] == mesh.areas[0]/3       
206       
207
208        #Test that points are arranged in a counter clock wise order
209        mesh.check_integrity()
210
211
212    def test_boundary_indices(self):
213        a = [0.0, 0.5]
214        b = [0.0, 0.0]               
215        c = [0.5, 0.5]
216
217        points = [a, b, c]
218        vertices = [ [0,1,2] ]
219        mesh = Mesh(points, vertices)       
220        mesh.check_integrity()
221
222        assert allclose(mesh.neighbours, [[-1,-2,-3]]) 
223               
224
225
226    def test_two_triangles(self):
227        a = [0.0, 0.0]
228        b = [0.0, 2.0]
229        c = [2.0,0.0]
230        e = [2.0, 2.0]
231        points = [a, b, c, e]
232        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
233        mesh = Mesh(points, vertices)
234       
235        assert mesh.areas[0] == 2.0
236
237        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])       
238
239
240        #Test that points are arranged in a counter clock wise order
241        mesh.check_integrity()       
242
243
244           
245    def test_more_triangles(self):
246       
247        a = [0.0, 0.0]
248        b = [0.0, 2.0]
249        c = [2.0,0.0]
250        d = [0.0, 4.0]
251        e = [2.0, 2.0]
252        f = [4.0,0.0]
253
254        points = [a, b, c, d, e, f]
255        #bac, bce, ecf, dbe, daf, dae
256        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
257        mesh = Mesh(points, vertices)       
258
259        #Test that points are arranged in a counter clock wise order
260        mesh.check_integrity()       
261
262        assert mesh.areas[0] == 2.0
263        assert mesh.areas[1] == 2.0
264        assert mesh.areas[2] == 2.0
265        assert mesh.areas[3] == 2.0
266        assert mesh.areas[4] == 8.0
267        assert mesh.areas[5] == 4.0                               
268       
269        assert mesh.edgelengths[1,0] == 2.0
270        assert mesh.edgelengths[1,1] == 2.0
271        assert mesh.edgelengths[1,2] == sqrt(8.0)
272
273        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
274        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3]) 
275        assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
276        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])         
277
278
279    def test_mesh_and_neighbours(self):
280        a = [0.0, 0.0]
281        b = [0.0, 2.0]
282        c = [2.0,0.0]
283        d = [0.0, 4.0]
284        e = [2.0, 2.0]
285        f = [4.0,0.0]
286
287
288        points = [a, b, c, d, e, f]
289
290        #bac, bce, ecf, dbe
291        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
292        mesh = Mesh(points, vertices)
293       
294        mesh.check_integrity()
295
296
297        T = mesh
298        tid = 0
299        assert T.number_of_boundaries[tid] == 2   
300        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
301        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
302        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)       
303
304        tid = 1
305        assert T.number_of_boundaries[tid] == 0           
306        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
307        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
308        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
309
310        tid = 2
311        assert T.number_of_boundaries[tid] == 2           
312        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
313        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
314        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)       
315
316        tid = 3
317        assert T.number_of_boundaries[tid] == 2           
318        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
319        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
320        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
321
322        #Neighbouring edges
323        tid = 0
324        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
325        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
326        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)       
327
328        tid = 1       
329        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
330        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
331        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
332
333        tid = 2       
334        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
335        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
336        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)       
337
338        tid = 3       
339        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
340        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
341        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
342
343       
344    def test_rectangular_mesh_basic(self):
345        M=1
346        N=1
347        mesh = rectangular_mesh(M, N)
348
349        #Test that points are arranged in a counter clock wise order       
350        mesh.check_integrity()
351
352        M=2
353        N=2
354        mesh = rectangular_mesh(M, N)
355       
356        #Test that points are arranged in a counter clock wise order       
357        mesh.check_integrity()
358
359        #assert mesh.boundary[(7,1)] == 2 # top
360        assert mesh.boundary[(7,1)] == 'top' # top
361        assert mesh.boundary[(3,1)] == 'top' # top               
362
363       
364    def test_boundary_tags(self):
365       
366
367        mesh = rectangular_mesh(4, 4)
368
369        #Test that points are arranged in a counter clock wise order       
370        mesh.check_integrity()
371
372        #print mesh.get_boundary_tags()
373        #print mesh.boundary
374
375        for k in [1, 3, 5, 7]:
376            assert mesh.boundary[(k,2)] == 'left'
377
378        for k in [24, 26, 28, 30]:
379            assert mesh.boundary[(k,2)] == 'right'           
380           
381        for k in [7, 15, 23, 31]:
382            assert mesh.boundary[(k,1)] == 'top'
383        for k in [0, 8, 16, 24]:
384            assert mesh.boundary[(k,1)] == 'bottom'
385
386
387       
388    def test_rectangular_mesh(self):
389        M=4
390        N=16
391        len1 = 100.0
392        len2 = 17.0
393        mesh = rectangular_mesh(M, N, len1, len2)
394
395        assert len(mesh) == 2*M*N       
396
397        for i in range(len(mesh)):
398            assert mesh.areas[i] == len1*len2/(2*M*N)
399
400            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
401            assert mesh.edgelengths[i, 0] == hypo
402            assert mesh.edgelengths[i, 1] == len1/M #x direction
403            assert mesh.edgelengths[i, 2] == len2/N #y direction
404           
405        #Test that points are arranged in a counter clock wise order     
406        mesh.check_integrity()   
407
408
409    def test_rectangular_mesh2(self):
410        #Check that integers don't cause trouble
411        N = 16
412        mesh = rectangular_mesh(2*N, N, len1=10, len2=10)
413
414
415    def test_surrogate_neighbours(self):
416        a = [0.0, 0.0]
417        b = [0.0, 2.0]
418        c = [2.0,0.0]
419        d = [0.0, 4.0]
420        e = [2.0, 2.0]
421        f = [4.0,0.0]
422
423        points = [a, b, c, d, e, f]
424
425        #bac, bce, ecf, dbe
426        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
427        mesh = Mesh(points, vertices)       
428        mesh.check_integrity()
429
430
431        T = mesh
432        tid = 0
433        assert T.number_of_boundaries[tid] == 2   
434        assert T.surrogate_neighbours[tid, 0] == tid
435        assert T.surrogate_neighbours[tid, 1] == 1   
436        assert T.surrogate_neighbours[tid, 2] == tid
437
438        tid = 1
439        assert T.number_of_boundaries[tid] == 0           
440        assert T.surrogate_neighbours[tid, 0] == 2
441        assert T.surrogate_neighbours[tid, 1] == 3
442        assert T.surrogate_neighbours[tid, 2] == 0
443
444        tid = 2
445        assert T.number_of_boundaries[tid] == 2           
446        assert T.surrogate_neighbours[tid, 0] == tid
447        assert T.surrogate_neighbours[tid, 1] == tid
448        assert T.surrogate_neighbours[tid, 2] == 1
449
450        tid = 3
451        assert T.number_of_boundaries[tid] == 2           
452        assert T.surrogate_neighbours[tid, 0] == 1 
453        assert T.surrogate_neighbours[tid, 1] == tid
454        assert T.surrogate_neighbours[tid, 2] == tid
455
456
457    def test_boundary_inputs(self):
458        a = [0.0, 0.0]
459        b = [0.0, 2.0]
460        c = [2.0,0.0]
461        d = [0.0, 4.0]
462        e = [2.0, 2.0]
463        f = [4.0,0.0]
464
465        points = [a, b, c, d, e, f]
466
467        #bac, bce, ecf, dbe
468        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
469
470        boundary = { (0, 0): 'First',
471                     (0, 2): 'Second',
472                     (2, 0): 'Third',
473                     (2, 1): 'Fourth',
474                     (3, 1): 'Fifth',
475                     (3, 2): 'Sixth'}                                         
476                     
477       
478        mesh = Mesh(points, vertices, boundary)       
479        mesh.check_integrity()
480
481
482        #Check enumeration
483        for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):       
484            b = -k-1
485            assert mesh.neighbours[vol_id, edge_id] == b
486
487
488
489    def test_boundary_inputs_using_defaults(self):
490        a = [0.0, 0.0]
491        b = [0.0, 2.0]
492        c = [2.0,0.0]
493        d = [0.0, 4.0]
494        e = [2.0, 2.0]
495        f = [4.0,0.0]
496
497        points = [a, b, c, d, e, f]
498
499        #bac, bce, ecf, dbe
500        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
501
502        boundary = { (0, 0): 'First',
503                     (0, 2): 'Second',
504                     (2, 0): 'Third',
505                     (2, 1): 'Fourth',
506                     #(3, 1): 'Fifth',  #Skip this
507                     (3, 2): 'Sixth'}                                         
508                     
509       
510        mesh = Mesh(points, vertices, boundary)       
511        mesh.check_integrity()
512
513        from config import default_boundary_tag
514        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
515
516
517        #Check enumeration
518        for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):       
519            b = -k-1
520            assert mesh.neighbours[vol_id, edge_id] == b
521
522
523
524
525
526
527    def test_inputs(self):
528        a = [0.0, 0.0]
529        b = [0.0, 2.0]
530        c = [2.0,0.0]
531        d = [0.0, 4.0]
532        e = [2.0, 2.0]
533        f = [4.0,0.0]
534
535        points = [a, b, c, d, e, f]
536
537        #bac, bce, ecf, dbe
538        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
539
540        #Too few points
541        try:
542            mesh = Mesh([points[0]], vertices)       
543        except AssertionError:
544            pass
545        else:
546            raise 'Should have raised an exception'
547
548        #Too few points - 1 element
549        try:
550            mesh = Mesh([points[0]], [vertices[0]])       
551        except AssertionError:
552            pass
553        else:
554            raise 'Should have raised an exception'       
555       
556        #Wrong dimension of vertices
557        try:
558            mesh = Mesh(points, vertices[0])       
559        except AssertionError:
560            pass
561        else:
562            raise 'Should have raised an exception'
563
564        #Unsubscriptable coordinates object raises exception
565        try:
566            mesh = Mesh(points[0], [vertices[0]])       
567        except AssertionError:
568            pass
569        else:
570            raise 'Should have raised an exception'
571
572        #Not specifying all boundary tags
573        try:
574            mesh = Mesh(points, vertices, {(3,0): 'x'})       
575        except AssertionError:
576            pass
577        else:
578            raise 'Should have raised an exception'
579
580        #Specifying wrong non existing segment
581        try:
582            mesh = Mesh(points, vertices, {(5,0): 'x'})       
583        except AssertionError:
584            pass
585        else:
586            raise 'Should have raised an exception'
587
588
589
590       
591       
592#-------------------------------------------------------------
593if __name__ == "__main__":
594    suite = unittest.makeSuite(TestCase,'test')
595    runner = unittest.TextTestRunner()
596    runner.run(suite)
597
598   
599   
600
601
602
Note: See TracBrowser for help on using the repository browser.