source: inundation/pyvolution/test_mesh.py @ 1911

Last change on this file since 1911 was 1911, checked in by ole, 18 years ago

Refactored pyvolution to use polygon functionality from new utilities package

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