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

Last change on this file since 707 was 648, checked in by ole, 20 years ago

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