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

Last change on this file since 1453 was 1392, checked in by steve, 20 years ago

Reapplied Matt's inscribed circle correction

File size: 23.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 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 Test_Mesh(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    def test_inscribed_circle_equilateral(self):
211        """test that the radius is calculated correctly by mesh in the case of an equilateral triangle"""
212        a = [0.0, 0.0]
213        b = [2.0, 0.0]
214        c = [1.0, sqrt(3.0)]
215
216        points = [a, b, c]
217        vertices = [[0,1,2]]
218
219        mesh = Mesh(points, vertices,use_inscribed_circle=False)
220        assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
221
222        mesh = Mesh(points, vertices,use_inscribed_circle=True)
223        assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
224
225    def test_inscribed_circle_rightangle_triangle(self):
226        """test that the radius is calculated correctly by mesh in the case of a right-angled triangle"""
227        a = [0.0, 0.0]
228        b = [4.0, 0.0]
229        c = [0.0, 3.0]
230
231        points = [a, b, c]
232        vertices = [[0,1,2]]
233
234        mesh = Mesh(points, vertices,use_inscribed_circle=False)
235        assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
236
237        mesh = Mesh(points, vertices,use_inscribed_circle=True)
238        assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
239
240
241    def test_two_triangles(self):
242        a = [0.0, 0.0]
243        b = [0.0, 2.0]
244        c = [2.0,0.0]
245        e = [2.0, 2.0]
246        points = [a, b, c, e]
247        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
248        mesh = Mesh(points, vertices)
249
250        assert mesh.areas[0] == 2.0
251
252        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
253
254
255        #Test that points are arranged in a counter clock wise order
256        mesh.check_integrity()
257
258
259
260    def test_more_triangles(self):
261
262        a = [0.0, 0.0]
263        b = [0.0, 2.0]
264        c = [2.0, 0.0]
265        d = [0.0, 4.0]
266        e = [2.0, 2.0]
267        f = [4.0, 0.0]
268
269        points = [a, b, c, d, e, f]
270        #bac, bce, ecf, dbe, daf, dae
271        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
272        mesh = Mesh(points, vertices)
273
274        #Test that points are arranged in a counter clock wise order
275        mesh.check_integrity()
276
277        assert mesh.areas[0] == 2.0
278        assert mesh.areas[1] == 2.0
279        assert mesh.areas[2] == 2.0
280        assert mesh.areas[3] == 2.0
281
282        assert mesh.edgelengths[1,0] == 2.0
283        assert mesh.edgelengths[1,1] == 2.0
284        assert mesh.edgelengths[1,2] == sqrt(8.0)
285
286        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
287        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
288        assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
289        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
290
291    def test_mesh_and_neighbours(self):
292        a = [0.0, 0.0]
293        b = [0.0, 2.0]
294        c = [2.0,0.0]
295        d = [0.0, 4.0]
296        e = [2.0, 2.0]
297        f = [4.0,0.0]
298
299
300        points = [a, b, c, d, e, f]
301
302        #bac, bce, ecf, dbe
303        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
304        mesh = Mesh(points, vertices)
305
306        mesh.check_integrity()
307
308
309        T = mesh
310        tid = 0
311        assert T.number_of_boundaries[tid] == 2
312        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
313        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
314        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)
315
316        tid = 1
317        assert T.number_of_boundaries[tid] == 0
318        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
319        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
320        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
321
322        tid = 2
323        assert T.number_of_boundaries[tid] == 2
324        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
325        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
326        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)
327
328        tid = 3
329        assert T.number_of_boundaries[tid] == 2
330        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
331        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
332        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
333
334        #Neighbouring edges
335        tid = 0
336        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
337        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
338        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)
339
340        tid = 1
341        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
342        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
343        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
344
345        tid = 2
346        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
347        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
348        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
349
350        tid = 3
351        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
352        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
353        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
354
355
356    def test_build_neighbour_structure_duplicates(self):
357        p0 = [-66.0, 14.0]
358        p1 = [14.0, -66.0]
359        p2 = [14.0, 14.0]
360        p3 = [60.0, 20.0]
361        p4 = [10.0, 60.0]
362        p5 = [60.0, 60.0]
363
364        points = [p0, p1, p2, p3, p4, p5]
365        triangles = [ [0, 1, 2],
366                      [3, 2, 1],
367                      [0, 2, 4],
368                      [0, 2, 4],
369                      [4, 2, 5],
370                      [5, 2, 3]]
371        try:
372            mesh = Mesh(points, triangles)
373        except:
374            pass
375        else:
376            raise "triangle edge duplicates not caught"
377
378    def test_rectangular_mesh_basic(self):
379        M=1
380        N=1
381
382        points, vertices, boundary = rectangular(M, N)
383        mesh = Mesh(points, vertices, boundary)
384
385        #Test that points are arranged in a counter clock wise order
386        mesh.check_integrity()
387
388        M=2
389        N=2
390        points, vertices, boundary = rectangular(M, N)
391        mesh = Mesh(points, vertices, boundary)
392
393        #Test that points are arranged in a counter clock wise order
394        mesh.check_integrity()
395
396        #assert mesh.boundary[(7,1)] == 2 # top
397        assert mesh.boundary[(7,1)] == 'top' # top
398        assert mesh.boundary[(3,1)] == 'top' # top
399
400
401    def test_boundary_tags(self):
402
403
404        points, vertices, boundary = rectangular(4, 4)
405        mesh = Mesh(points, vertices, boundary)
406
407
408        #Test that points are arranged in a counter clock wise order
409        mesh.check_integrity()
410
411        #print mesh.get_boundary_tags()
412        #print mesh.boundary
413
414        for k in [1, 3, 5, 7]:
415            assert mesh.boundary[(k,2)] == 'left'
416
417        for k in [24, 26, 28, 30]:
418            assert mesh.boundary[(k,2)] == 'right'
419
420        for k in [7, 15, 23, 31]:
421            assert mesh.boundary[(k,1)] == 'top'
422        for k in [0, 8, 16, 24]:
423            assert mesh.boundary[(k,1)] == 'bottom'
424
425
426
427    def test_rectangular_mesh(self):
428        M=4
429        N=16
430        len1 = 100.0
431        len2 = 17.0
432
433        points, vertices, boundary = rectangular(M, N, len1, len2)
434        mesh = Mesh(points, vertices, boundary)
435
436        assert len(mesh) == 2*M*N
437
438        for i in range(len(mesh)):
439            assert mesh.areas[i] == len1*len2/(2*M*N)
440
441            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
442            assert mesh.edgelengths[i, 0] == hypo
443            assert mesh.edgelengths[i, 1] == len1/M #x direction
444            assert mesh.edgelengths[i, 2] == len2/N #y direction
445
446        #Test that points are arranged in a counter clock wise order
447        mesh.check_integrity()
448
449
450    def test_rectangular_mesh2(self):
451        #Check that integers don't cause trouble
452        N = 16
453
454        points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
455        mesh = Mesh(points, vertices, boundary)
456
457
458
459    def test_surrogate_neighbours(self):
460        a = [0.0, 0.0]
461        b = [0.0, 2.0]
462        c = [2.0,0.0]
463        d = [0.0, 4.0]
464        e = [2.0, 2.0]
465        f = [4.0,0.0]
466
467        points = [a, b, c, d, e, f]
468
469        #bac, bce, ecf, dbe
470        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
471        mesh = Mesh(points, vertices)
472        mesh.check_integrity()
473
474
475        T = mesh
476        tid = 0
477        assert T.number_of_boundaries[tid] == 2
478        assert T.surrogate_neighbours[tid, 0] == tid
479        assert T.surrogate_neighbours[tid, 1] == 1
480        assert T.surrogate_neighbours[tid, 2] == tid
481
482        tid = 1
483        assert T.number_of_boundaries[tid] == 0
484        assert T.surrogate_neighbours[tid, 0] == 2
485        assert T.surrogate_neighbours[tid, 1] == 3
486        assert T.surrogate_neighbours[tid, 2] == 0
487
488        tid = 2
489        assert T.number_of_boundaries[tid] == 2
490        assert T.surrogate_neighbours[tid, 0] == tid
491        assert T.surrogate_neighbours[tid, 1] == tid
492        assert T.surrogate_neighbours[tid, 2] == 1
493
494        tid = 3
495        assert T.number_of_boundaries[tid] == 2
496        assert T.surrogate_neighbours[tid, 0] == 1
497        assert T.surrogate_neighbours[tid, 1] == tid
498        assert T.surrogate_neighbours[tid, 2] == tid
499
500
501    def test_boundary_inputs(self):
502        a = [0.0, 0.0]
503        b = [0.0, 2.0]
504        c = [2.0,0.0]
505        d = [0.0, 4.0]
506        e = [2.0, 2.0]
507        f = [4.0,0.0]
508
509        points = [a, b, c, d, e, f]
510
511        #bac, bce, ecf, dbe
512        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
513
514        boundary = { (0, 0): 'First',
515                     (0, 2): 'Second',
516                     (2, 0): 'Third',
517                     (2, 1): 'Fourth',
518                     (3, 1): 'Fifth',
519                     (3, 2): 'Sixth'}
520
521
522        mesh = Mesh(points, vertices, boundary)
523        mesh.check_integrity()
524
525
526        #Check enumeration
527        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
528        #    b = -k-1
529        #    assert mesh.neighbours[vol_id, edge_id] == b
530
531
532
533    def test_boundary_inputs_using_one_default(self):
534        a = [0.0, 0.0]
535        b = [0.0, 2.0]
536        c = [2.0,0.0]
537        d = [0.0, 4.0]
538        e = [2.0, 2.0]
539        f = [4.0,0.0]
540
541        points = [a, b, c, d, e, f]
542
543        #bac, bce, ecf, dbe
544        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
545
546        boundary = { (0, 0): 'First',
547                     (0, 2): 'Second',
548                     (2, 0): 'Third',
549                     (2, 1): 'Fourth',
550                     #(3, 1): 'Fifth',  #Skip this
551                     (3, 2): 'Sixth'}
552
553
554        mesh = Mesh(points, vertices, boundary)
555        mesh.check_integrity()
556
557        from config import default_boundary_tag
558        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
559
560
561        #Check enumeration
562        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
563        #    b = -k-1
564        #    assert mesh.neighbours[vol_id, edge_id] == b
565
566    def test_boundary_inputs_using_all_defaults(self):
567        a = [0.0, 0.0]
568        b = [0.0, 2.0]
569        c = [2.0,0.0]
570        d = [0.0, 4.0]
571        e = [2.0, 2.0]
572        f = [4.0,0.0]
573
574        points = [a, b, c, d, e, f]
575
576        #bac, bce, ecf, dbe
577        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
578
579        boundary = { (0, 0): 'First',
580                     (0, 2): 'Second',
581                     (2, 0): 'Third',
582                     (2, 1): 'Fourth',
583                     #(3, 1): 'Fifth',  #Skip this
584                     (3, 2): 'Sixth'}
585
586
587        mesh = Mesh(points, vertices) #, boundary)
588        mesh.check_integrity()
589
590        from config import default_boundary_tag
591        assert mesh.boundary[ (0, 0) ] == default_boundary_tag
592        assert mesh.boundary[ (0, 2) ] == default_boundary_tag
593        assert mesh.boundary[ (2, 0) ] == default_boundary_tag
594        assert mesh.boundary[ (2, 1) ] == default_boundary_tag
595        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
596        assert mesh.boundary[ (3, 2) ] == default_boundary_tag
597
598
599        #Check enumeration
600        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
601        #    b = -k-1
602        #    assert mesh.neighbours[vol_id, edge_id] == b
603
604
605
606
607
608
609    def test_inputs(self):
610        a = [0.0, 0.0]
611        b = [0.0, 2.0]
612        c = [2.0,0.0]
613        d = [0.0, 4.0]
614        e = [2.0, 2.0]
615        f = [4.0,0.0]
616
617        points = [a, b, c, d, e, f]
618
619        #bac, bce, ecf, dbe
620        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
621
622        #Too few points
623        try:
624            mesh = Mesh([points[0]], vertices)
625        except AssertionError:
626            pass
627        else:
628            raise 'Should have raised an exception'
629
630        #Too few points - 1 element
631        try:
632            mesh = Mesh([points[0]], [vertices[0]])
633        except AssertionError:
634            pass
635        else:
636            raise 'Should have raised an exception'
637
638        #Wrong dimension of vertices
639        try:
640            mesh = Mesh(points, vertices[0])
641        except AssertionError:
642            pass
643        else:
644            raise 'Should have raised an exception'
645
646        #Unsubscriptable coordinates object raises exception
647        try:
648            mesh = Mesh(points[0], [vertices[0]])
649        except AssertionError:
650            pass
651        else:
652            raise 'Should have raised an exception'
653
654        #FIXME: This has been commented out pending a decision
655        #whether to allow partial boundary tags or not
656        #
657        #Not specifying all boundary tags
658        #try:
659        #    mesh = Mesh(points, vertices, {(3,0): 'x'})
660        #except AssertionError:
661        #    pass
662        #else:
663        #    raise 'Should have raised an exception'
664
665        #Specifying wrong non existing segment
666        try:
667            mesh = Mesh(points, vertices, {(5,0): 'x'})
668        except AssertionError:
669            pass
670        else:
671            raise 'Should have raised an exception'
672
673
674
675
676    def test_internal_boundaries(self):
677        """
678        get values based on triangle lists.
679        """
680        from mesh_factory import rectangular
681        from shallow_water import Domain
682        from Numeric import zeros, Float
683
684        #Create basic mesh
685        points, vertices, boundary = rectangular(1, 3)
686
687        # Add an internal boundary
688        boundary[(2,0)] = 'internal'
689        boundary[(1,0)] = 'internal'
690
691        #Create shallow water domain
692        domain = Domain(points, vertices, boundary)
693        domain.build_tagged_elements_dictionary({'bottom':[0,1],
694                                                 'top':[4,5],
695                                                 'all':[0,1,2,3,4,5]})
696
697
698    def test_boundary_polygon(self):
699        from mesh_factory import rectangular
700        from mesh import Mesh
701        from Numeric import zeros, Float
702        from util import inside_polygon
703
704        #Create basic mesh
705        points, vertices, boundary = rectangular(2, 2)
706        mesh = Mesh(points, vertices, boundary)
707
708
709        P = mesh.get_boundary_polygon()
710
711        assert len(P) == 8
712        assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
713                            [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
714                            [0.0, 1.0], [0.0, 0.5]])
715        for p in points:
716            #print p, P
717            assert inside_polygon(p, P)
718
719
720    def test_boundary_polygon_II(self):
721        from mesh import Mesh
722        from Numeric import zeros, Float
723        from util import inside_polygon
724
725        #Points
726        a = [0.0, 0.0] #0
727        b = [0.0, 0.5] #1
728        c = [0.0, 1.0] #2
729        d = [0.5, 0.0] #3
730        e = [0.5, 0.5] #4
731        f = [1.0, 0.0] #5
732        g = [1.0, 0.5] #6
733        h = [1.0, 1.0] #7
734        i = [1.5, 0.5] #8
735
736        points = [a, b, c, d, e, f, g, h, i]
737
738        #dea, bae, bec, fgd,
739        #edg, ghe, gfi, gih
740        vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
741                     [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
742
743        mesh = Mesh(points, vertices)
744
745        mesh.check_integrity()
746
747        P = mesh.get_boundary_polygon()
748
749        assert len(P) == 8
750        assert allclose(P, [a, d, f, i, h, e, c, b])
751
752        for p in points:
753            #print p, P
754            assert inside_polygon(p, P)
755
756
757    def test_boundary_polygon_III(self):
758        """Same as II but vertices ordered differently
759        """
760
761        from mesh import Mesh
762        from Numeric import zeros, Float
763        from util import inside_polygon
764
765        #Points
766        a = [0.0, 0.0] #0
767        b = [0.0, 0.5] #1
768        c = [0.0, 1.0] #2
769        d = [0.5, 0.0] #3
770        e = [0.5, 0.5] #4
771        f = [1.0, 0.0] #5
772        g = [1.0, 0.5] #6
773        h = [1.0, 1.0] #7
774        i = [1.5, 0.5] #8
775
776        points = [a, b, c, d, e, f, g, h, i]
777
778        #edg, ghe, gfi, gih
779        #dea, bae, bec, fgd,
780        vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
781                    [3,4,0], [1,0,4], [1,4,2], [5,6,3]]
782
783
784        mesh = Mesh(points, vertices)
785        mesh.check_integrity()
786
787
788        P = mesh.get_boundary_polygon()
789
790        assert len(P) == 8
791        assert allclose(P, [a, d, f, i, h, e, c, b])
792
793        for p in points:
794            assert inside_polygon(p, P)
795
796
797
798#-------------------------------------------------------------
799if __name__ == "__main__":
800    suite = unittest.makeSuite(Test_Mesh,'test')
801    runner = unittest.TextTestRunner()
802    runner.run(suite)
803
804
805
806
Note: See TracBrowser for help on using the repository browser.