source: inundation/pyvolution/test_mesh.py @ 2709

Last change on this file since 2709 was 2709, checked in by ole, 19 years ago

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

File size: 25.2 KB
Line 
1#!/usr/bin/env python
2
3
4
5#FIXME: Seperate the tests for mesh and general_mesh
6
7#FIXME (Ole): Maxe this test independent of anything that inherits from General_mesh (namely shallow_water)
8
9import unittest
10from math import sqrt
11
12from mesh import *
13from mesh_factory import rectangular
14from config import epsilon
15from Numeric import allclose, array
16
17from utilities.polygon import inside_polygon
18
19def distance(x, y):
20    return sqrt( sum( (array(x)-array(y))**2 ))
21
22class Test_Mesh(unittest.TestCase):
23    def setUp(self):
24        pass
25
26    def tearDown(self):
27        pass
28
29    def test_triangle_inputs(self):
30        points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
31        vertices = [0,1,2] #Wrong
32
33        try:
34            mesh = Mesh(points, vertices)
35        except:
36            pass
37        else:
38            msg = 'Should have raised exception'
39            raise msg
40
41
42    def test_basic_triangle(self):
43
44        a = [0.0, 0.0]
45        b = [4.0, 0.0]
46        c = [0.0, 3.0]
47
48        points = [a, b, c]
49        vertices = [[0,1,2]]
50        mesh = Mesh(points, vertices)
51
52        #Centroid
53        centroid = mesh.centroid_coordinates[0]
54        assert centroid[0] == 4.0/3
55        assert centroid[1] == 1.0
56
57        #Area
58        assert mesh.areas[0] == 6.0,\
59               'Area was %f, should have been 6.0' %mesh.areas[0]
60
61        #Normals
62        normals = mesh.get_normals()
63        assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
64        assert allclose(normals[0, 2:4], [-1.0, 0.0])
65        assert allclose(normals[0, 4:6], [0.0, -1.0])
66
67        assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
68        assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
69        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])
70
71        #Edge lengths
72        assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
73
74
75        #Vertex coordinates
76        V = mesh.get_vertex_coordinates()
77        assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
78
79        V = mesh.get_vertex_coordinates(obj=True)
80        assert allclose(V, [ [0.0, 0.0],
81                             [4.0, 0.0],
82                             [0.0, 3.0] ])
83
84        V0 = mesh.get_vertex_coordinate(0, 0)
85        assert allclose(V0, [0.0, 0.0])
86
87        V1 = mesh.get_vertex_coordinate(0, 1)
88        assert allclose(V1, [4.0, 0.0])
89
90        V2 = mesh.get_vertex_coordinate(0, 2)
91        assert allclose(V2, [0.0, 3.0])
92
93
94        #General tests:
95
96        #Test that points are arranged in a counter clock wise order etc
97        mesh.check_integrity()
98
99
100        #Test that the centroid is located 2/3 of the way
101        #from each vertex to the midpoint of the opposite side
102
103        V = mesh.get_vertex_coordinates()
104
105        x0 = V[0,0]
106        y0 = V[0,1]
107        x1 = V[0,2]
108        y1 = V[0,3]
109        x2 = V[0,4]
110        y2 = V[0,5]
111
112        m0 = [(x1 + x2)/2, (y1 + y2)/2]
113        m1 = [(x0 + x2)/2, (y0 + y2)/2]
114        m2 = [(x1 + x0)/2, (y1 + y0)/2]
115
116        d0 = distance(centroid, [x0, y0])
117        d1 = distance(m0, [x0, y0])
118        assert d0 == 2*d1/3
119        #
120        d0 = distance(centroid, [x1, y1])
121        d1 = distance(m1, [x1, y1])
122        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
123
124        d0 = distance(centroid, [x2, y2])
125        d1 = distance(m2, [x2, y2])
126        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
127
128        #Radius
129        d0 = distance(centroid, m0)
130        assert d0 == 5.0/6
131
132        d1 = distance(centroid, m1)
133        assert d1 == sqrt(73.0/36)
134
135        d2 = distance(centroid, m2)
136        assert d2 == sqrt(13.0/9)
137
138        assert mesh.radii[0] == min(d0, d1, d2)
139        assert mesh.radii[0] == 5.0/6
140
141
142        #Let x be the centroid of triangle abc.
143        #Test that areas of the three triangles axc, cxb, and bxa are equal.
144        points = [a, b, c, centroid]
145        vertices = [[0,3,2], [2,3,1], [1,3,0]]
146        new_mesh = Mesh(points, vertices)
147
148        assert new_mesh.areas[0] == new_mesh.areas[1]
149        assert new_mesh.areas[1] == new_mesh.areas[2]
150        assert new_mesh.areas[1] == new_mesh.areas[2]
151
152        assert new_mesh.areas[1] == mesh.areas[0]/3
153
154
155
156    def test_general_triangle(self):
157        a = [2.0, 1.0]
158        b = [6.0, 2.0]
159        c = [1.0, 3.0]
160
161        points = [a, b, c]
162        vertices = [[0,1,2]]
163
164        mesh = Mesh(points, vertices)
165        centroid = mesh.centroid_coordinates[0]
166
167
168        #Test that the centroid is located 2/3 of the way
169        #from each vertex to the midpoint of the opposite side
170
171        V = mesh.get_vertex_coordinates()
172
173        x0 = V[0,0]
174        y0 = V[0,1]
175        x1 = V[0,2]
176        y1 = V[0,3]
177        x2 = V[0,4]
178        y2 = V[0,5]
179
180        m0 = [(x1 + x2)/2, (y1 + y2)/2]
181        m1 = [(x0 + x2)/2, (y0 + y2)/2]
182        m2 = [(x1 + x0)/2, (y1 + y0)/2]
183
184        d0 = distance(centroid, [x0, y0])
185        d1 = distance(m0, [x0, y0])
186        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
187        #
188        d0 = distance(centroid, [x1, y1])
189        d1 = distance(m1, [x1, y1])
190        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
191
192        d0 = distance(centroid, [x2, y2])
193        d1 = distance(m2, [x2, y2])
194        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
195
196        #Radius
197        d0 = distance(centroid, m0)
198        d1 = distance(centroid, m1)
199        d2 = distance(centroid, m2)
200        assert mesh.radii[0] == min(d0, d1, d2)
201
202
203
204        #Let x be the centroid of triangle abc.
205        #Test that areas of the three triangles axc, cxb, and bxa are equal.
206
207        points = [a, b, c, centroid]
208        vertices = [[0,3,2], [2,3,1], [1,3,0]]
209        new_mesh = Mesh(points, vertices)
210
211        assert new_mesh.areas[0] == new_mesh.areas[1]
212        assert new_mesh.areas[1] == new_mesh.areas[2]
213        assert new_mesh.areas[1] == new_mesh.areas[2]
214
215        assert new_mesh.areas[1] == mesh.areas[0]/3
216
217
218        #Test that points are arranged in a counter clock wise order
219        mesh.check_integrity()
220
221    def test_inscribed_circle_equilateral(self):
222        """test that the radius is calculated correctly by mesh in the case of an equilateral triangle"""
223        a = [0.0, 0.0]
224        b = [2.0, 0.0]
225        c = [1.0, sqrt(3.0)]
226
227        points = [a, b, c]
228        vertices = [[0,1,2]]
229
230        mesh = Mesh(points, vertices,use_inscribed_circle=False)
231        assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
232
233        mesh = Mesh(points, vertices,use_inscribed_circle=True)
234        assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
235
236    def test_inscribed_circle_rightangle_triangle(self):
237        """test that the radius is calculated correctly by mesh in the case of a right-angled triangle"""
238        a = [0.0, 0.0]
239        b = [4.0, 0.0]
240        c = [0.0, 3.0]
241
242        points = [a, b, c]
243        vertices = [[0,1,2]]
244
245        mesh = Mesh(points, vertices,use_inscribed_circle=False)
246        assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
247
248        mesh = Mesh(points, vertices,use_inscribed_circle=True)
249        assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
250
251
252    def test_two_triangles(self):
253        a = [0.0, 0.0]
254        b = [0.0, 2.0]
255        c = [2.0,0.0]
256        e = [2.0, 2.0]
257        points = [a, b, c, e]
258        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
259        mesh = Mesh(points, vertices)
260
261        assert mesh.areas[0] == 2.0
262
263        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
264
265
266        #Test that points are arranged in a counter clock wise order
267        mesh.check_integrity()
268
269
270
271    def test_more_triangles(self):
272
273        a = [0.0, 0.0]
274        b = [0.0, 2.0]
275        c = [2.0, 0.0]
276        d = [0.0, 4.0]
277        e = [2.0, 2.0]
278        f = [4.0, 0.0]
279
280        points = [a, b, c, d, e, f]
281        #bac, bce, ecf, dbe, daf, dae
282        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
283        mesh = Mesh(points, vertices)
284
285        #Test that points are arranged in a counter clock wise order
286        mesh.check_integrity()
287
288        assert mesh.areas[0] == 2.0
289        assert mesh.areas[1] == 2.0
290        assert mesh.areas[2] == 2.0
291        assert mesh.areas[3] == 2.0
292
293        assert mesh.edgelengths[1,0] == 2.0
294        assert mesh.edgelengths[1,1] == 2.0
295        assert mesh.edgelengths[1,2] == sqrt(8.0)
296
297        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
298        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
299        assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
300        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
301
302    def test_mesh_and_neighbours(self):
303        a = [0.0, 0.0]
304        b = [0.0, 2.0]
305        c = [2.0,0.0]
306        d = [0.0, 4.0]
307        e = [2.0, 2.0]
308        f = [4.0,0.0]
309
310
311        points = [a, b, c, d, e, f]
312
313        #bac, bce, ecf, dbe
314        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
315        mesh = Mesh(points, vertices)
316
317        mesh.check_integrity()
318
319
320        T = mesh
321        tid = 0
322        assert T.number_of_boundaries[tid] == 2
323        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
324        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
325        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)
326
327        tid = 1
328        assert T.number_of_boundaries[tid] == 0
329        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
330        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
331        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
332
333        tid = 2
334        assert T.number_of_boundaries[tid] == 2
335        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
336        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
337        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)
338
339        tid = 3
340        assert T.number_of_boundaries[tid] == 2
341        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
342        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
343        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
344
345        #Neighbouring edges
346        tid = 0
347        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
348        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
349        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)
350
351        tid = 1
352        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
353        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
354        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
355
356        tid = 2
357        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
358        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
359        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
360
361        tid = 3
362        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
363        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
364        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
365
366
367    def test_build_neighbour_structure_duplicates(self):
368        p0 = [-66.0, 14.0]
369        p1 = [14.0, -66.0]
370        p2 = [14.0, 14.0]
371        p3 = [60.0, 20.0]
372        p4 = [10.0, 60.0]
373        p5 = [60.0, 60.0]
374
375        points = [p0, p1, p2, p3, p4, p5]
376        triangles = [ [0, 1, 2],
377                      [3, 2, 1],
378                      [0, 2, 4],
379                      [0, 2, 4],
380                      [4, 2, 5],
381                      [5, 2, 3]]
382        try:
383            mesh = Mesh(points, triangles)
384        except:
385            pass
386        else:
387            raise "triangle edge duplicates not caught"
388
389    def test_rectangular_mesh_basic(self):
390        M=1
391        N=1
392
393        points, vertices, boundary = rectangular(M, N)
394        mesh = Mesh(points, vertices, boundary)
395
396        #Test that points are arranged in a counter clock wise order
397        mesh.check_integrity()
398
399        M=2
400        N=2
401        points, vertices, boundary = rectangular(M, N)
402        mesh = Mesh(points, vertices, boundary)
403
404        #Test that points are arranged in a counter clock wise order
405        mesh.check_integrity()
406
407        #assert mesh.boundary[(7,1)] == 2 # top
408        assert mesh.boundary[(7,1)] == 'top' # top
409        assert mesh.boundary[(3,1)] == 'top' # top
410
411
412    def test_boundary_tags(self):
413
414
415        points, vertices, boundary = rectangular(4, 4)
416        mesh = Mesh(points, vertices, boundary)
417
418
419        #Test that points are arranged in a counter clock wise order
420        mesh.check_integrity()
421
422        #print mesh.get_boundary_tags()
423        #print mesh.boundary
424
425        for k in [1, 3, 5, 7]:
426            assert mesh.boundary[(k,2)] == 'left'
427
428        for k in [24, 26, 28, 30]:
429            assert mesh.boundary[(k,2)] == 'right'
430
431        for k in [7, 15, 23, 31]:
432            assert mesh.boundary[(k,1)] == 'top'
433        for k in [0, 8, 16, 24]:
434            assert mesh.boundary[(k,1)] == 'bottom'
435
436
437
438    def test_rectangular_mesh(self):
439        M=4
440        N=16
441        len1 = 100.0
442        len2 = 17.0
443
444        points, vertices, boundary = rectangular(M, N, len1, len2)
445        mesh = Mesh(points, vertices, boundary)
446
447        assert len(mesh) == 2*M*N
448
449        for i in range(len(mesh)):
450            assert mesh.areas[i] == len1*len2/(2*M*N)
451
452            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
453            assert mesh.edgelengths[i, 0] == hypo
454            assert mesh.edgelengths[i, 1] == len1/M #x direction
455            assert mesh.edgelengths[i, 2] == len2/N #y direction
456
457        #Test that points are arranged in a counter clock wise order
458        mesh.check_integrity()
459
460
461    def test_rectangular_mesh2(self):
462        #Check that integers don't cause trouble
463        N = 16
464
465        points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
466        mesh = Mesh(points, vertices, boundary)
467
468
469
470    def test_surrogate_neighbours(self):
471        a = [0.0, 0.0]
472        b = [0.0, 2.0]
473        c = [2.0,0.0]
474        d = [0.0, 4.0]
475        e = [2.0, 2.0]
476        f = [4.0,0.0]
477
478        points = [a, b, c, d, e, f]
479
480        #bac, bce, ecf, dbe
481        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
482        mesh = Mesh(points, vertices)
483        mesh.check_integrity()
484
485
486        T = mesh
487        tid = 0
488        assert T.number_of_boundaries[tid] == 2
489        assert T.surrogate_neighbours[tid, 0] == tid
490        assert T.surrogate_neighbours[tid, 1] == 1
491        assert T.surrogate_neighbours[tid, 2] == tid
492
493        tid = 1
494        assert T.number_of_boundaries[tid] == 0
495        assert T.surrogate_neighbours[tid, 0] == 2
496        assert T.surrogate_neighbours[tid, 1] == 3
497        assert T.surrogate_neighbours[tid, 2] == 0
498
499        tid = 2
500        assert T.number_of_boundaries[tid] == 2
501        assert T.surrogate_neighbours[tid, 0] == tid
502        assert T.surrogate_neighbours[tid, 1] == tid
503        assert T.surrogate_neighbours[tid, 2] == 1
504
505        tid = 3
506        assert T.number_of_boundaries[tid] == 2
507        assert T.surrogate_neighbours[tid, 0] == 1
508        assert T.surrogate_neighbours[tid, 1] == tid
509        assert T.surrogate_neighbours[tid, 2] == tid
510
511
512    def test_boundary_inputs(self):
513        a = [0.0, 0.0]
514        b = [0.0, 2.0]
515        c = [2.0,0.0]
516        d = [0.0, 4.0]
517        e = [2.0, 2.0]
518        f = [4.0,0.0]
519
520        points = [a, b, c, d, e, f]
521
522        #bac, bce, ecf, dbe
523        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
524
525        boundary = { (0, 0): 'First',
526                     (0, 2): 'Second',
527                     (2, 0): 'Third',
528                     (2, 1): 'Fourth',
529                     (3, 1): 'Fifth',
530                     (3, 2): 'Sixth'}
531
532
533        mesh = Mesh(points, vertices, boundary)
534        mesh.check_integrity()
535
536
537        #Check enumeration
538        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
539        #    b = -k-1
540        #    assert mesh.neighbours[vol_id, edge_id] == b
541
542
543
544    def test_boundary_inputs_using_one_default(self):
545        a = [0.0, 0.0]
546        b = [0.0, 2.0]
547        c = [2.0,0.0]
548        d = [0.0, 4.0]
549        e = [2.0, 2.0]
550        f = [4.0,0.0]
551
552        points = [a, b, c, d, e, f]
553
554        #bac, bce, ecf, dbe
555        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
556
557        boundary = { (0, 0): 'First',
558                     (0, 2): 'Second',
559                     (2, 0): 'Third',
560                     (2, 1): 'Fourth',
561                     #(3, 1): 'Fifth',  #Skip this
562                     (3, 2): 'Sixth'}
563
564
565        mesh = Mesh(points, vertices, boundary)
566        mesh.check_integrity()
567
568        from config import default_boundary_tag
569        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
570
571
572        #Check enumeration
573        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
574        #    b = -k-1
575        #    assert mesh.neighbours[vol_id, edge_id] == b
576
577    def test_boundary_inputs_using_all_defaults(self):
578        a = [0.0, 0.0]
579        b = [0.0, 2.0]
580        c = [2.0,0.0]
581        d = [0.0, 4.0]
582        e = [2.0, 2.0]
583        f = [4.0,0.0]
584
585        points = [a, b, c, d, e, f]
586
587        #bac, bce, ecf, dbe
588        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
589
590        boundary = { (0, 0): 'First',
591                     (0, 2): 'Second',
592                     (2, 0): 'Third',
593                     (2, 1): 'Fourth',
594                     #(3, 1): 'Fifth',  #Skip this
595                     (3, 2): 'Sixth'}
596
597
598        mesh = Mesh(points, vertices) #, boundary)
599        mesh.check_integrity()
600
601        from config import default_boundary_tag
602        assert mesh.boundary[ (0, 0) ] == default_boundary_tag
603        assert mesh.boundary[ (0, 2) ] == default_boundary_tag
604        assert mesh.boundary[ (2, 0) ] == default_boundary_tag
605        assert mesh.boundary[ (2, 1) ] == default_boundary_tag
606        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
607        assert mesh.boundary[ (3, 2) ] == default_boundary_tag
608
609
610        #Check enumeration
611        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
612        #    b = -k-1
613        #    assert mesh.neighbours[vol_id, edge_id] == b
614
615
616
617
618
619
620    def test_inputs(self):
621        a = [0.0, 0.0]
622        b = [0.0, 2.0]
623        c = [2.0,0.0]
624        d = [0.0, 4.0]
625        e = [2.0, 2.0]
626        f = [4.0,0.0]
627
628        points = [a, b, c, d, e, f]
629
630        #bac, bce, ecf, dbe
631        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
632
633        #Too few points
634        try:
635            mesh = Mesh([points[0]], vertices)
636        except AssertionError:
637            pass
638        else:
639            raise 'Should have raised an exception'
640
641        #Too few points - 1 element
642        try:
643            mesh = Mesh([points[0]], [vertices[0]])
644        except AssertionError:
645            pass
646        else:
647            raise 'Should have raised an exception'
648
649        #Wrong dimension of vertices
650        try:
651            mesh = Mesh(points, vertices[0])
652        except AssertionError:
653            pass
654        else:
655            raise 'Should have raised an exception'
656
657        #Unsubscriptable coordinates object raises exception
658        try:
659            mesh = Mesh(points[0], [vertices[0]])
660        except AssertionError:
661            pass
662        else:
663            raise 'Should have raised an exception'
664
665        #FIXME: This has been commented out pending a decision
666        #whether to allow partial boundary tags or not
667        #
668        #Not specifying all boundary tags
669        #try:
670        #    mesh = Mesh(points, vertices, {(3,0): 'x'})
671        #except AssertionError:
672        #    pass
673        #else:
674        #    raise 'Should have raised an exception'
675
676        #Specifying wrong non existing segment
677        try:
678            mesh = Mesh(points, vertices, {(5,0): 'x'})
679        except AssertionError:
680            pass
681        else:
682            raise 'Should have raised an exception'
683
684
685
686
687    def test_internal_boundaries(self):
688        """
689        get values based on triangle lists.
690        """
691        from mesh_factory import rectangular
692        from shallow_water import Domain
693        from Numeric import zeros, Float
694
695        #Create basic mesh
696        points, vertices, boundary = rectangular(1, 3)
697
698        # Add an internal boundary
699        boundary[(2,0)] = 'internal'
700        boundary[(1,0)] = 'internal'
701
702        #Create shallow water domain
703        domain = Domain(points, vertices, boundary)
704        domain.build_tagged_elements_dictionary({'bottom':[0,1],
705                                                 'top':[4,5],
706                                                 'all':[0,1,2,3,4,5]})
707
708
709    def test_boundary_polygon(self):
710        from mesh_factory import rectangular
711        from mesh import Mesh
712        from Numeric import zeros, Float
713
714        #Create basic mesh
715        points, vertices, boundary = rectangular(2, 2)
716        mesh = Mesh(points, vertices, boundary)
717
718
719        P = mesh.get_boundary_polygon()
720
721        assert len(P) == 8
722        assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
723                            [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
724                            [0.0, 1.0], [0.0, 0.5]])
725        for p in points:
726            #print p, P
727            assert inside_polygon(p, P)
728
729
730    def test_boundary_polygon_II(self):
731        from mesh import Mesh
732        from Numeric import zeros, Float
733       
734
735        #Points
736        a = [0.0, 0.0] #0
737        b = [0.0, 0.5] #1
738        c = [0.0, 1.0] #2
739        d = [0.5, 0.0] #3
740        e = [0.5, 0.5] #4
741        f = [1.0, 0.0] #5
742        g = [1.0, 0.5] #6
743        h = [1.0, 1.0] #7
744        i = [1.5, 0.5] #8
745
746        points = [a, b, c, d, e, f, g, h, i]
747
748        #dea, bae, bec, fgd,
749        #edg, ghe, gfi, gih
750        vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
751                     [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
752
753        mesh = Mesh(points, vertices)
754
755        mesh.check_integrity()
756
757        P = mesh.get_boundary_polygon()
758
759        assert len(P) == 8
760        assert allclose(P, [a, d, f, i, h, e, c, b])
761
762        for p in points:
763            #print p, P
764            assert inside_polygon(p, P)
765
766
767    def test_boundary_polygon_III(self):
768        """Same as II but vertices ordered differently
769        """
770
771        from mesh import Mesh
772        from Numeric import zeros, Float
773
774
775        #Points
776        a = [0.0, 0.0] #0
777        b = [0.0, 0.5] #1
778        c = [0.0, 1.0] #2
779        d = [0.5, 0.0] #3
780        e = [0.5, 0.5] #4
781        f = [1.0, 0.0] #5
782        g = [1.0, 0.5] #6
783        h = [1.0, 1.0] #7
784        i = [1.5, 0.5] #8
785
786        points = [a, b, c, d, e, f, g, h, i]
787
788        #edg, ghe, gfi, gih
789        #dea, bae, bec, fgd,
790        vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
791                    [3,4,0], [1,0,4], [1,4,2], [5,6,3]]
792
793
794        mesh = Mesh(points, vertices)
795        mesh.check_integrity()
796
797
798        P = mesh.get_boundary_polygon()
799
800        assert len(P) == 8
801        assert allclose(P, [a, d, f, i, h, e, c, b])
802
803        for p in points:
804            assert inside_polygon(p, P)
805
806
807    def test_boundary_polygon_IV(self):
808        """Reproduce test test_spatio_temporal_file_function_time
809        from test_util.py that looked as if it produced the wrong boundary
810        """
811
812        from mesh import Mesh
813        from Numeric import zeros, Float
814        from mesh_factory import rectangular       
815
816        #Create a domain to hold test grid
817        #(0:15, -20:10)
818        points, vertices, boundary =\
819                rectangular(4, 4, 15, 30, origin = (0, -20))       
820
821        #####
822        mesh = Mesh(points, vertices)
823        mesh.check_integrity()
824
825        P = mesh.get_boundary_polygon()
826
827        #print P
828        assert len(P) == 16
829        for p in points:
830            assert inside_polygon(p, P)
831
832
833
834        #####
835        mesh = Mesh(points, vertices, boundary)
836        mesh.check_integrity()
837
838        P = mesh.get_boundary_polygon()
839
840       
841        #print P, len(P)
842        assert len(P) == 16
843
844        for p in points:
845            assert inside_polygon(p, P)
846
847        #print mesh.statistics()   
848
849
850
851    def test_boundary_polygon_V(self):
852        """Create a discontinuous mesh (duplicate vertices)
853        and check that boundary is as expected
854       
855        """
856        from mesh import Mesh
857        from Numeric import zeros, Float
858       
859
860        #Points
861        a = [0.0, 0.0] #0
862        b = [0.0, 0.5] #1
863        c = [0.0, 1.0] #2
864        d = [0.5, 0.0] #3
865        e = [0.5, 0.5] #4
866        f = [1.0, 0.0] #5
867        g = [1.0, 0.5] #6
868        h = [1.0, 1.0] #7
869        i = [1.5, 0.5] #8
870
871        #Duplicate points for triangles edg [4,3,6] (central) and
872        #gid [6,8,7] (top right boundary) to them disconnected
873        #from the others
874
875        e0 = [0.5, 0.5] #9
876        d0 = [0.5, 0.0] #10
877        g0 = [1.0, 0.5] #11
878        i0 = [1.5, 0.5] #12       
879       
880
881        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
882
883
884
885        #dea, bae, bec, fgd,
886        #edg, ghe, gfi, gih
887        #vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
888        #             [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
889
890
891        #dea, bae, bec, fgd,
892        #e0d0g0, ghe, gfi, g0i0h
893        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
894                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]       
895
896        mesh = Mesh(points, vertices)
897
898        mesh.check_integrity()
899
900        P = mesh.get_boundary_polygon()
901
902        #print P
903       
904        assert len(P) == 8
905        assert allclose(P, [a, d, f, i, h, e, c, b])
906        assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
907       
908
909        for p in points:
910            #print p, P
911            assert inside_polygon(p, P)
912
913
914#-------------------------------------------------------------
915if __name__ == "__main__":
916    suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
917    #suite = unittest.makeSuite(Test_Mesh,'test')
918    runner = unittest.TextTestRunner()
919    runner.run(suite)
920
921
922
923
Note: See TracBrowser for help on using the repository browser.