source: storm_surge/pyvolution/test_mesh.py @ 81

Last change on this file since 81 was 1, checked in by duncan, 20 years ago

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