source: anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py @ 3688

Last change on this file since 3688 was 3688, checked in by ole, 17 years ago

Fixed get_boundary polygon based on information from new unit test.
It now works for all examples of meshes provided.

File size: 33.9 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 neighbour_mesh import *
13from mesh_factory import rectangular
14from anuga.config import epsilon
15from Numeric import allclose, array, Int
16
17from anuga.coordinate_transforms.geo_reference import Geo_reference
18from anuga.utilities.polygon import is_inside_polygon
19from anuga.utilities.numerical_tools import ensure_numeric
20
21def distance(x, y):
22    return sqrt( sum( (array(x)-array(y))**2 ))
23
24class Test_Mesh(unittest.TestCase):
25    def setUp(self):
26        pass
27
28    def tearDown(self):
29        pass
30
31    def test_triangle_inputs(self):
32        points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
33        vertices = [0,1,2] #Wrong
34
35        try:
36            mesh = Mesh(points, vertices)
37        except:
38            pass
39        else:
40            msg = 'Should have raised exception'
41            raise msg
42
43
44    def test_basic_triangle(self):
45
46        a = [0.0, 0.0]
47        b = [4.0, 0.0]
48        c = [0.0, 3.0]
49
50        points = [a, b, c]
51        vertices = [[0,1,2]]
52        mesh = Mesh(points, vertices)
53
54        #Centroid
55        centroid = mesh.centroid_coordinates[0]
56        assert centroid[0] == 4.0/3
57        assert centroid[1] == 1.0
58
59        #Area
60        assert mesh.areas[0] == 6.0,\
61               'Area was %f, should have been 6.0' %mesh.areas[0]
62
63        #Normals
64        normals = mesh.get_normals()
65        assert allclose(normals[0, 0:2], [3.0/5, 4.0/5])
66        assert allclose(normals[0, 2:4], [-1.0, 0.0])
67        assert allclose(normals[0, 4:6], [0.0, -1.0])
68
69        assert allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
70        assert allclose(mesh.get_normal(0,1), [-1.0, 0.0])
71        assert allclose(mesh.get_normal(0,2), [0.0, -1.0])
72
73        #Edge lengths
74        assert allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
75
76
77        #Vertex coordinates
78        V = mesh.get_vertex_coordinates()
79        assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
80
81        V = mesh.get_vertex_coordinates(obj=True)
82        assert allclose(V, [ [0.0, 0.0],
83                             [4.0, 0.0],
84                             [0.0, 3.0] ])
85
86        V0 = mesh.get_vertex_coordinate(0, 0)
87        assert allclose(V0, [0.0, 0.0])
88
89        V1 = mesh.get_vertex_coordinate(0, 1)
90        assert allclose(V1, [4.0, 0.0])
91
92        V2 = mesh.get_vertex_coordinate(0, 2)
93        assert allclose(V2, [0.0, 3.0])
94
95
96        #General tests:
97
98        #Test that points are arranged in a counter clock wise order etc
99        mesh.check_integrity()
100
101
102        #Test that the centroid is located 2/3 of the way
103        #from each vertex to the midpoint of the opposite side
104
105        V = mesh.get_vertex_coordinates()
106
107        x0 = V[0,0]
108        y0 = V[0,1]
109        x1 = V[0,2]
110        y1 = V[0,3]
111        x2 = V[0,4]
112        y2 = V[0,5]
113
114        m0 = [(x1 + x2)/2, (y1 + y2)/2]
115        m1 = [(x0 + x2)/2, (y0 + y2)/2]
116        m2 = [(x1 + x0)/2, (y1 + y0)/2]
117
118        d0 = distance(centroid, [x0, y0])
119        d1 = distance(m0, [x0, y0])
120        assert d0 == 2*d1/3
121        #
122        d0 = distance(centroid, [x1, y1])
123        d1 = distance(m1, [x1, y1])
124        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
125
126        d0 = distance(centroid, [x2, y2])
127        d1 = distance(m2, [x2, y2])
128        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
129
130        #Radius
131        d0 = distance(centroid, m0)
132        assert d0 == 5.0/6
133
134        d1 = distance(centroid, m1)
135        assert d1 == sqrt(73.0/36)
136
137        d2 = distance(centroid, m2)
138        assert d2 == sqrt(13.0/9)
139
140        assert mesh.radii[0] == min(d0, d1, d2)
141        assert mesh.radii[0] == 5.0/6
142
143
144        #Let x be the centroid of triangle abc.
145        #Test that areas of the three triangles axc, cxb, and bxa are equal.
146        points = [a, b, c, centroid]
147        vertices = [[0,3,2], [2,3,1], [1,3,0]]
148        new_mesh = Mesh(points, vertices)
149
150        assert new_mesh.areas[0] == new_mesh.areas[1]
151        assert new_mesh.areas[1] == new_mesh.areas[2]
152        assert new_mesh.areas[1] == new_mesh.areas[2]
153
154        assert new_mesh.areas[1] == mesh.areas[0]/3
155
156
157
158    def test_general_triangle(self):
159        a = [2.0, 1.0]
160        b = [6.0, 2.0]
161        c = [1.0, 3.0]
162
163        points = [a, b, c]
164        vertices = [[0,1,2]]
165
166        mesh = Mesh(points, vertices)
167        centroid = mesh.centroid_coordinates[0]
168
169
170        #Test that the centroid is located 2/3 of the way
171        #from each vertex to the midpoint of the opposite side
172
173        V = mesh.get_vertex_coordinates()
174
175        x0 = V[0,0]
176        y0 = V[0,1]
177        x1 = V[0,2]
178        y1 = V[0,3]
179        x2 = V[0,4]
180        y2 = V[0,5]
181
182        m0 = [(x1 + x2)/2, (y1 + y2)/2]
183        m1 = [(x0 + x2)/2, (y0 + y2)/2]
184        m2 = [(x1 + x0)/2, (y1 + y0)/2]
185
186        d0 = distance(centroid, [x0, y0])
187        d1 = distance(m0, [x0, y0])
188        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
189        #
190        d0 = distance(centroid, [x1, y1])
191        d1 = distance(m1, [x1, y1])
192        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
193
194        d0 = distance(centroid, [x2, y2])
195        d1 = distance(m2, [x2, y2])
196        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
197
198        #Radius
199        d0 = distance(centroid, m0)
200        d1 = distance(centroid, m1)
201        d2 = distance(centroid, m2)
202        assert mesh.radii[0] == min(d0, d1, d2)
203
204
205
206        #Let x be the centroid of triangle abc.
207        #Test that areas of the three triangles axc, cxb, and bxa are equal.
208
209        points = [a, b, c, centroid]
210        vertices = [[0,3,2], [2,3,1], [1,3,0]]
211        new_mesh = Mesh(points, vertices)
212
213        assert new_mesh.areas[0] == new_mesh.areas[1]
214        assert new_mesh.areas[1] == new_mesh.areas[2]
215        assert new_mesh.areas[1] == new_mesh.areas[2]
216
217        assert new_mesh.areas[1] == mesh.areas[0]/3
218
219
220        #Test that points are arranged in a counter clock wise order
221        mesh.check_integrity()
222
223    def test_inscribed_circle_equilateral(self):
224        """test that the radius is calculated correctly by mesh in the case of an equilateral triangle"""
225        a = [0.0, 0.0]
226        b = [2.0, 0.0]
227        c = [1.0, sqrt(3.0)]
228
229        points = [a, b, c]
230        vertices = [[0,1,2]]
231
232        mesh = Mesh(points, vertices,use_inscribed_circle=False)
233        assert allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
234
235        mesh = Mesh(points, vertices,use_inscribed_circle=True)
236        assert allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
237
238    def test_inscribed_circle_rightangle_triangle(self):
239        """test that the radius is calculated correctly by mesh in the case of a right-angled triangle"""
240        a = [0.0, 0.0]
241        b = [4.0, 0.0]
242        c = [0.0, 3.0]
243
244        points = [a, b, c]
245        vertices = [[0,1,2]]
246
247        mesh = Mesh(points, vertices,use_inscribed_circle=False)
248        assert allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
249
250        mesh = Mesh(points, vertices,use_inscribed_circle=True)
251        assert allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
252
253
254    def test_two_triangles(self):
255        a = [0.0, 0.0]
256        b = [0.0, 2.0]
257        c = [2.0,0.0]
258        e = [2.0, 2.0]
259        points = [a, b, c, e]
260        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
261        mesh = Mesh(points, vertices)
262
263        assert mesh.areas[0] == 2.0
264
265        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
266
267
268        #Test that points are arranged in a counter clock wise order
269        mesh.check_integrity()
270
271
272
273    def test_more_triangles(self):
274
275        a = [0.0, 0.0]
276        b = [0.0, 2.0]
277        c = [2.0, 0.0]
278        d = [0.0, 4.0]
279        e = [2.0, 2.0]
280        f = [4.0, 0.0]
281
282        points = [a, b, c, d, e, f]
283        #bac, bce, ecf, dbe, daf, dae
284        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
285        mesh = Mesh(points, vertices)
286
287        #Test that points are arranged in a counter clock wise order
288        mesh.check_integrity()
289
290        assert mesh.areas[0] == 2.0
291        assert mesh.areas[1] == 2.0
292        assert mesh.areas[2] == 2.0
293        assert mesh.areas[3] == 2.0
294
295        assert mesh.edgelengths[1,0] == 2.0
296        assert mesh.edgelengths[1,1] == 2.0
297        assert mesh.edgelengths[1,2] == sqrt(8.0)
298
299        assert allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
300        assert allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
301        assert allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
302        assert allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
303
304    def test_mesh_and_neighbours(self):
305        a = [0.0, 0.0]
306        b = [0.0, 2.0]
307        c = [2.0,0.0]
308        d = [0.0, 4.0]
309        e = [2.0, 2.0]
310        f = [4.0,0.0]
311
312
313        points = [a, b, c, d, e, f]
314
315        #bac, bce, ecf, dbe
316        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
317        mesh = Mesh(points, vertices)
318
319        mesh.check_integrity()
320
321
322        T = mesh
323        tid = 0
324        assert T.number_of_boundaries[tid] == 2
325        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
326        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
327        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)
328
329        tid = 1
330        assert T.number_of_boundaries[tid] == 0
331        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
332        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
333        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
334
335        tid = 2
336        assert T.number_of_boundaries[tid] == 2
337        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
338        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
339        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)
340
341        tid = 3
342        assert T.number_of_boundaries[tid] == 2
343        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
344        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
345        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
346
347        #Neighbouring edges
348        tid = 0
349        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
350        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
351        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)
352
353        tid = 1
354        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
355        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
356        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
357
358        tid = 2
359        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
360        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
361        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
362
363        tid = 3
364        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
365        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
366        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
367
368
369    def test_build_neighbour_structure_duplicates(self):
370        p0 = [-66.0, 14.0]
371        p1 = [14.0, -66.0]
372        p2 = [14.0, 14.0]
373        p3 = [60.0, 20.0]
374        p4 = [10.0, 60.0]
375        p5 = [60.0, 60.0]
376
377        points = [p0, p1, p2, p3, p4, p5]
378        triangles = [ [0, 1, 2],
379                      [3, 2, 1],
380                      [0, 2, 4],
381                      [0, 2, 4],
382                      [4, 2, 5],
383                      [5, 2, 3]]
384        try:
385            mesh = Mesh(points, triangles)
386        except:
387            pass
388        else:
389            raise "triangle edge duplicates not caught"
390
391    def test_rectangular_mesh_basic(self):
392        M=1
393        N=1
394
395        points, vertices, boundary = rectangular(M, N)
396        mesh = Mesh(points, vertices, boundary)
397
398        #Test that points are arranged in a counter clock wise order
399        mesh.check_integrity()
400
401        M=2
402        N=2
403        points, vertices, boundary = rectangular(M, N)
404        mesh = Mesh(points, vertices, boundary)
405
406        #Test that points are arranged in a counter clock wise order
407        mesh.check_integrity()
408
409        #assert mesh.boundary[(7,1)] == 2 # top
410        assert mesh.boundary[(7,1)] == 'top' # top
411        assert mesh.boundary[(3,1)] == 'top' # top
412
413
414    def test_boundary_tags(self):
415
416
417        points, vertices, boundary = rectangular(4, 4)
418        mesh = Mesh(points, vertices, boundary)
419
420
421        #Test that points are arranged in a counter clock wise order
422        mesh.check_integrity()
423
424        #print mesh.get_boundary_tags()
425        #print mesh.boundary
426
427        for k in [1, 3, 5, 7]:
428            assert mesh.boundary[(k,2)] == 'left'
429
430        for k in [24, 26, 28, 30]:
431            assert mesh.boundary[(k,2)] == 'right'
432
433        for k in [7, 15, 23, 31]:
434            assert mesh.boundary[(k,1)] == 'top'
435        for k in [0, 8, 16, 24]:
436            assert mesh.boundary[(k,1)] == 'bottom'
437
438
439
440    def test_rectangular_mesh(self):
441        M=4
442        N=16
443        len1 = 100.0
444        len2 = 17.0
445
446        points, vertices, boundary = rectangular(M, N, len1, len2)
447        mesh = Mesh(points, vertices, boundary)
448
449        assert len(mesh) == 2*M*N
450
451        for i in range(len(mesh)):
452            assert mesh.areas[i] == len1*len2/(2*M*N)
453
454            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
455            assert mesh.edgelengths[i, 0] == hypo
456            assert mesh.edgelengths[i, 1] == len1/M #x direction
457            assert mesh.edgelengths[i, 2] == len2/N #y direction
458
459        #Test that points are arranged in a counter clock wise order
460        mesh.check_integrity()
461
462
463    def test_rectangular_mesh2(self):
464        #Check that integers don't cause trouble
465        N = 16
466
467        points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
468        mesh = Mesh(points, vertices, boundary)
469
470
471
472    def test_surrogate_neighbours(self):
473        a = [0.0, 0.0]
474        b = [0.0, 2.0]
475        c = [2.0,0.0]
476        d = [0.0, 4.0]
477        e = [2.0, 2.0]
478        f = [4.0,0.0]
479
480        points = [a, b, c, d, e, f]
481
482        #bac, bce, ecf, dbe
483        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
484        mesh = Mesh(points, vertices)
485        mesh.check_integrity()
486
487
488        T = mesh
489        tid = 0
490        assert T.number_of_boundaries[tid] == 2
491        assert T.surrogate_neighbours[tid, 0] == tid
492        assert T.surrogate_neighbours[tid, 1] == 1
493        assert T.surrogate_neighbours[tid, 2] == tid
494
495        tid = 1
496        assert T.number_of_boundaries[tid] == 0
497        assert T.surrogate_neighbours[tid, 0] == 2
498        assert T.surrogate_neighbours[tid, 1] == 3
499        assert T.surrogate_neighbours[tid, 2] == 0
500
501        tid = 2
502        assert T.number_of_boundaries[tid] == 2
503        assert T.surrogate_neighbours[tid, 0] == tid
504        assert T.surrogate_neighbours[tid, 1] == tid
505        assert T.surrogate_neighbours[tid, 2] == 1
506
507        tid = 3
508        assert T.number_of_boundaries[tid] == 2
509        assert T.surrogate_neighbours[tid, 0] == 1
510        assert T.surrogate_neighbours[tid, 1] == tid
511        assert T.surrogate_neighbours[tid, 2] == tid
512
513
514    def test_boundary_inputs(self):
515        a = [0.0, 0.0]
516        b = [0.0, 2.0]
517        c = [2.0,0.0]
518        d = [0.0, 4.0]
519        e = [2.0, 2.0]
520        f = [4.0,0.0]
521
522        points = [a, b, c, d, e, f]
523
524        #bac, bce, ecf, dbe
525        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
526
527        boundary = { (0, 0): 'First',
528                     (0, 2): 'Second',
529                     (2, 0): 'Third',
530                     (2, 1): 'Fourth',
531                     (3, 1): 'Fifth',
532                     (3, 2): 'Sixth'}
533
534
535        mesh = Mesh(points, vertices, boundary)
536        mesh.check_integrity()
537
538
539        #Check enumeration
540        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
541        #    b = -k-1
542        #    assert mesh.neighbours[vol_id, edge_id] == b
543
544
545
546    def test_boundary_inputs_using_one_default(self):
547        a = [0.0, 0.0]
548        b = [0.0, 2.0]
549        c = [2.0,0.0]
550        d = [0.0, 4.0]
551        e = [2.0, 2.0]
552        f = [4.0,0.0]
553
554        points = [a, b, c, d, e, f]
555
556        #bac, bce, ecf, dbe
557        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
558
559        boundary = { (0, 0): 'First',
560                     (0, 2): 'Second',
561                     (2, 0): 'Third',
562                     (2, 1): 'Fourth',
563                     #(3, 1): 'Fifth',  #Skip this
564                     (3, 2): 'Sixth'}
565
566
567        mesh = Mesh(points, vertices, boundary)
568        mesh.check_integrity()
569
570        from anuga.config import default_boundary_tag
571        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
572
573
574        #Check enumeration
575        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
576        #    b = -k-1
577        #    assert mesh.neighbours[vol_id, edge_id] == b
578
579    def test_boundary_inputs_using_all_defaults(self):
580        a = [0.0, 0.0]
581        b = [0.0, 2.0]
582        c = [2.0,0.0]
583        d = [0.0, 4.0]
584        e = [2.0, 2.0]
585        f = [4.0,0.0]
586
587        points = [a, b, c, d, e, f]
588
589        #bac, bce, ecf, dbe
590        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
591
592        boundary = { (0, 0): 'First',
593                     (0, 2): 'Second',
594                     (2, 0): 'Third',
595                     (2, 1): 'Fourth',
596                     #(3, 1): 'Fifth',  #Skip this
597                     (3, 2): 'Sixth'}
598
599
600        mesh = Mesh(points, vertices) #, boundary)
601        mesh.check_integrity()
602
603        from anuga.config import default_boundary_tag
604        assert mesh.boundary[ (0, 0) ] == default_boundary_tag
605        assert mesh.boundary[ (0, 2) ] == default_boundary_tag
606        assert mesh.boundary[ (2, 0) ] == default_boundary_tag
607        assert mesh.boundary[ (2, 1) ] == default_boundary_tag
608        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
609        assert mesh.boundary[ (3, 2) ] == default_boundary_tag
610
611
612        #Check enumeration
613        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
614        #    b = -k-1
615        #    assert mesh.neighbours[vol_id, edge_id] == b
616
617
618
619
620
621
622    def test_inputs(self):
623        a = [0.0, 0.0]
624        b = [0.0, 2.0]
625        c = [2.0,0.0]
626        d = [0.0, 4.0]
627        e = [2.0, 2.0]
628        f = [4.0,0.0]
629
630        points = [a, b, c, d, e, f]
631
632        #bac, bce, ecf, dbe
633        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
634
635        #Too few points
636        try:
637            mesh = Mesh([points[0]], vertices)
638        except AssertionError:
639            pass
640        else:
641            raise 'Should have raised an exception'
642
643        #Too few points - 1 element
644        try:
645            mesh = Mesh([points[0]], [vertices[0]])
646        except AssertionError:
647            pass
648        else:
649            raise 'Should have raised an exception'
650
651        #Wrong dimension of vertices
652        try:
653            mesh = Mesh(points, vertices[0])
654        except AssertionError:
655            pass
656        else:
657            raise 'Should have raised an exception'
658
659        #Unsubscriptable coordinates object raises exception
660        try:
661            mesh = Mesh(points[0], [vertices[0]])
662        except AssertionError:
663            pass
664        else:
665            raise 'Should have raised an exception'
666
667        #FIXME: This has been commented out pending a decision
668        #whether to allow partial boundary tags or not
669        #
670        #Not specifying all boundary tags
671        #try:
672        #    mesh = Mesh(points, vertices, {(3,0): 'x'})
673        #except AssertionError:
674        #    pass
675        #else:
676        #    raise 'Should have raised an exception'
677
678        #Specifying wrong non existing segment
679        try:
680            mesh = Mesh(points, vertices, {(5,0): 'x'})
681        except AssertionError:
682            pass
683        else:
684            raise 'Should have raised an exception'
685
686
687
688
689    def test_internal_boundaries(self):
690        """
691        get values based on triangle lists.
692        """
693        from mesh_factory import rectangular
694        from Numeric import zeros, Float
695
696        #Create basic mesh
697        points, vertices, boundary = rectangular(1, 3)
698
699        # Add an internal boundary
700        boundary[(2,0)] = 'internal'
701        boundary[(1,0)] = 'internal'
702
703        #Create shallow water domain
704        domain = Mesh(points, vertices, boundary)
705        domain.build_tagged_elements_dictionary({'bottom':[0,1],
706                                                 'top':[4,5],
707                                                 'all':[0,1,2,3,4,5]})
708
709
710    def test_boundary_polygon(self):
711        from mesh_factory import rectangular
712        #from mesh import Mesh
713        from Numeric import zeros, Float
714
715        #Create basic mesh
716        points, vertices, boundary = rectangular(2, 2)
717        mesh = Mesh(points, vertices, boundary)
718
719
720        P = mesh.get_boundary_polygon()
721
722        assert len(P) == 8
723        assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
724                            [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
725                            [0.0, 1.0], [0.0, 0.5]])
726        for p in points:
727            #print p, P
728            assert is_inside_polygon(p, P)
729
730
731    def test_boundary_polygon_II(self):
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 is_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 Numeric import zeros, Float
772
773
774        #Points
775        a = [0.0, 0.0] #0
776        b = [0.0, 0.5] #1
777        c = [0.0, 1.0] #2
778        d = [0.5, 0.0] #3
779        e = [0.5, 0.5] #4
780        f = [1.0, 0.0] #5
781        g = [1.0, 0.5] #6
782        h = [1.0, 1.0] #7
783        i = [1.5, 0.5] #8
784
785        points = [a, b, c, d, e, f, g, h, i]
786
787        #edg, ghe, gfi, gih
788        #dea, bae, bec, fgd,
789        vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
790                    [3,4,0], [1,0,4], [1,4,2], [5,6,3]]
791
792
793        mesh = Mesh(points, vertices)
794        mesh.check_integrity()
795
796
797        P = mesh.get_boundary_polygon()
798
799        assert len(P) == 8
800        assert allclose(P, [a, d, f, i, h, e, c, b])
801
802        for p in points:
803            assert is_inside_polygon(p, P)
804
805
806    def test_boundary_polygon_IV(self):
807        """Reproduce test test_spatio_temporal_file_function_time
808        from test_util.py that looked as if it produced the wrong boundary
809        """
810
811        from Numeric import zeros, Float
812        from mesh_factory import rectangular       
813
814        #Create a domain to hold test grid
815        #(0:15, -20:10)
816        points, vertices, boundary =\
817                rectangular(4, 4, 15, 30, origin = (0, -20))       
818
819        #####
820        mesh = Mesh(points, vertices)
821        mesh.check_integrity()
822
823        P = mesh.get_boundary_polygon()
824
825        #print P
826        assert len(P) == 16
827        for p in points:
828            assert is_inside_polygon(p, P)
829
830
831
832        #####
833        mesh = Mesh(points, vertices, boundary)
834        mesh.check_integrity()
835
836        P = mesh.get_boundary_polygon()
837
838       
839        #print P, len(P)
840        assert len(P) == 16
841
842        for p in points:
843            assert is_inside_polygon(p, P)
844
845        #print mesh.statistics()   
846
847
848
849    def test_boundary_polygon_V(self):
850        """Create a discontinuous mesh (duplicate vertices)
851        and check that boundary is as expected
852       
853        """
854        from Numeric import zeros, Float
855       
856
857        #Points
858        a = [0.0, 0.0] #0
859        b = [0.0, 0.5] #1
860        c = [0.0, 1.0] #2
861        d = [0.5, 0.0] #3
862        e = [0.5, 0.5] #4
863        f = [1.0, 0.0] #5
864        g = [1.0, 0.5] #6
865        h = [1.0, 1.0] #7
866        i = [1.5, 0.5] #8
867
868        #Duplicate points for triangles edg [4,3,6] (central) and
869        #gid [6,8,7] (top right boundary) to them disconnected
870        #from the others
871
872        e0 = [0.5, 0.5] #9
873        d0 = [0.5, 0.0] #10
874        g0 = [1.0, 0.5] #11
875        i0 = [1.5, 0.5] #12       
876       
877
878        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
879
880
881
882        #dea, bae, bec, fgd,
883        #edg, ghe, gfi, gih
884        #vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
885        #             [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
886
887
888        #dea, bae, bec, fgd,
889        #e0d0g0, ghe, gfi, g0i0h
890        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
891                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]       
892
893        mesh = Mesh(points, vertices)
894
895        mesh.check_integrity()
896
897        P = mesh.get_boundary_polygon()
898
899        #print P
900       
901        assert len(P) == 8
902        assert allclose(P, [a, d, f, i, h, e, c, b])
903        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)])
904       
905
906        for p in points:
907            #print p, P
908            assert is_inside_polygon(p, P)
909
910
911
912    def test_boundary_polygon_VI(self):
913        """test_boundary_polygon_VI(self)
914
915        Create a discontinuous mesh (duplicate vertices) from a real situation that failed
916        and check that boundary is as expected
917        """
918
919
920        from anuga.utilities.polygon import plot_polygons
921
922        # First do the continuous version of mesh
923
924        points = [[   6626.85400391,      0.        ],
925                  [      0.        ,  38246.4140625 ],
926                  [   9656.2734375 ,  68351.265625  ],
927                  [  20827.25585938,  77818.203125  ],
928                  [  32755.59375   ,  58126.9765625 ],
929                  [  35406.3359375 ,  79332.9140625 ],
930                  [  31998.23828125,  88799.84375   ],
931                  [  23288.65820313, 104704.296875  ],
932                  [  32187.57617188, 109816.4375    ],
933                  [  50364.08984375, 110763.1328125 ],
934                  [  80468.9453125 ,  96184.0546875 ],
935                  [  86149.1015625 , 129886.34375   ],
936                  [ 118715.359375  , 129886.34375   ],
937                  [ 117768.6640625 ,  85770.4296875 ],
938                  [ 101485.5390625 ,  45251.9453125 ],
939                  [  49985.4140625 ,   2272.06396484],
940                  [  51737.94140625,  90559.2109375 ],
941                  [  56659.0703125 ,  65907.6796875 ],
942                  [  75735.4765625 ,  23762.00585938],
943                  [  52341.70703125,  38563.39453125]]
944
945        ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
946
947        triangles = [[19, 0,15],
948                     [ 2, 4, 3],
949                     [ 4, 2, 1],
950                     [ 1,19, 4],
951                     [15,18,19],
952                     [18,14,17],
953                     [19, 1, 0],
954                     [ 6, 8, 7],
955                     [ 8, 6,16],
956                     [10, 9,16],
957                     [17, 5, 4],
958                     [16,17,10],
959                     [17,19,18],
960                     [ 5,17,16],
961                     [10,14,13],
962                     [10,17,14],
963                     [ 8,16, 9],
964                     [12,11,10],
965                     [10,13,12],
966                     [19,17, 4],
967                     [16, 6, 5]]
968
969        mesh = Mesh(points, triangles)
970        mesh.check_integrity()
971        Pref = mesh.get_boundary_polygon()
972
973        #plot_polygons([ensure_numeric(Pref)], 'goodP')
974
975        for p in points:
976            assert is_inside_polygon(p, Pref)
977       
978       
979        # Then do the discontinuous version
980        import warnings
981        warnings.filterwarnings('ignore')
982
983       
984        points = [[  52341.70703125,  38563.39453125],
985                  [   6626.85400391,      0.        ],
986                  [  49985.4140625 ,   2272.06396484],
987                  [   9656.2734375 ,  68351.265625  ],
988                  [  32755.59375   ,  58126.9765625 ],
989                  [  20827.25585938,  77818.203125  ],
990                  [  32755.59375   ,  58126.9765625 ],
991                  [   9656.2734375 ,  68351.265625  ],
992                  [      0.        ,  38246.4140625 ],
993                  [      0.        ,  38246.4140625 ],
994                  [  52341.70703125,  38563.39453125],
995                  [  32755.59375   ,  58126.9765625 ],
996                  [  49985.4140625 ,   2272.06396484],
997                  [  75735.4765625 ,  23762.00585938],
998                  [  52341.70703125,  38563.39453125],
999                  [  75735.4765625 ,  23762.00585938],
1000                  [ 101485.5390625 ,  45251.9453125 ],
1001                  [  56659.0703125 ,  65907.6796875 ],
1002                  [  52341.70703125,  38563.39453125],
1003                  [      0.        ,  38246.4140625 ],
1004                  [   6626.85400391,      0.        ],
1005                  [  31998.23828125,  88799.84375   ],
1006                  [  32187.57617188, 109816.4375    ],
1007                  [  23288.65820313, 104704.296875  ],
1008                  [  32187.57617188, 109816.4375    ],
1009                  [  31998.23828125,  88799.84375   ],
1010                  [  51737.94140625,  90559.2109375 ],
1011                  [  80468.9453125 ,  96184.0546875 ],
1012                  [  50364.08984375, 110763.1328125 ],
1013                  [  51737.94140625,  90559.2109375 ],
1014                  [  56659.0703125 ,  65907.6796875 ],
1015                  [  35406.3359375 ,  79332.9140625 ],
1016                  [  32755.59375   ,  58126.9765625 ],
1017                  [  51737.94140625,  90559.2109375 ],
1018                  [  56659.0703125 ,  65907.6796875 ],
1019                  [  80468.9453125 ,  96184.0546875 ],
1020                  [  56659.0703125 ,  65907.6796875 ],
1021                  [  52341.70703125,  38563.39453125],
1022                  [  75735.4765625 ,  23762.00585938],
1023                  [  35406.3359375 ,  79332.9140625 ],
1024                  [  56659.0703125 ,  65907.6796875 ],
1025                  [  51737.94140625,  90559.2109375 ],
1026                  [  80468.9453125 ,  96184.0546875 ],
1027                  [ 101485.5390625 ,  45251.9453125 ],
1028                  [ 117768.6640625 ,  85770.4296875 ],
1029                  [  80468.9453125 ,  96184.0546875 ],
1030                  [  56659.0703125 ,  65907.6796875 ],
1031                  [ 101485.5390625 ,  45251.9453125 ],
1032                  [  32187.57617188, 109816.4375    ],
1033                  [  51737.94140625,  90559.2109375 ],
1034                  [  50364.08984375, 110763.1328125 ],
1035                  [ 118715.359375  , 129886.34375   ],
1036                  [  86149.1015625 , 129886.34375   ],
1037                  [  80468.9453125 ,  96184.0546875 ],
1038                  [  80468.9453125 ,  96184.0546875 ],
1039                  [ 117768.6640625 ,  85770.4296875 ],
1040                  [ 118715.359375  , 129886.34375   ],
1041                  [  52341.70703125,  38563.39453125],
1042                  [  56659.0703125 ,  65907.6796875 ],
1043                  [  32755.59375   ,  58126.9765625 ],
1044                  [  51737.94140625,  90559.2109375 ],
1045                  [  31998.23828125,  88799.84375   ],
1046                  [  35406.3359375 ,  79332.9140625 ]]
1047
1048        scaled_points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
1049
1050        triangles = [[ 0, 1, 2],
1051                     [ 3, 4, 5],
1052                     [ 6, 7, 8],
1053                     [ 9,10,11],
1054                     [12,13,14],
1055                     [15,16,17],
1056                     [18,19,20],
1057                     [21,22,23],
1058                     [24,25,26],
1059                     [27,28,29],
1060                     [30,31,32],
1061                     [33,34,35],
1062                     [36,37,38],
1063                     [39,40,41],
1064                     [42,43,44],
1065                     [45,46,47],
1066                     [48,49,50],
1067                     [51,52,53],
1068                     [54,55,56],
1069                     [57,58,59],
1070                     [60,61,62]]
1071
1072
1073        # First use scaled points for ease of debugging
1074        mesh = Mesh(scaled_points, triangles)
1075        mesh.check_integrity()
1076        P = mesh.get_boundary_polygon()
1077
1078        for p in scaled_points:
1079            assert is_inside_polygon(p, P)           
1080
1081        # Then use original points and test       
1082        mesh = Mesh(points, triangles)
1083        mesh.check_integrity()
1084        P = mesh.get_boundary_polygon()
1085
1086        for p in points:
1087            assert is_inside_polygon(p, P)           
1088
1089        assert allclose(P, Pref)   
1090
1091    def test_lone_vertices(self):
1092        a = [2.0, 1.0]
1093        b = [6.0, 2.0]
1094        c = [1.0, 3.0]
1095        d = [2.0, 4.0]
1096
1097        points = [a, b, c, d]
1098        vertices = [[0,1,2]]
1099
1100        mesh = Mesh(points, vertices)
1101        mesh.check_integrity()
1102        loners = mesh.get_lone_vertices()
1103        self.failUnless(loners==[3],
1104                        'FAILED!')
1105
1106       
1107        a = [2.0, 1.0]
1108        b = [6.0, 2.0]
1109        c = [1.0, 3.0]
1110        d = [2.0, 4.0]
1111
1112        points = [d, a, b, c]
1113        vertices = [[3,1,2]]
1114
1115        mesh = Mesh(points, vertices)
1116        mesh.check_integrity()
1117        loners = mesh.get_lone_vertices()
1118        self.failUnless(loners==[0],
1119                        'FAILED!') 
1120
1121    def test_mesh_get_boundary_polygon_with_georeferencing(self):
1122     
1123        # test
1124        a = [0.0, 0.0]
1125        b = [4.0, 0.0]
1126        c = [0.0, 4.0]
1127
1128        absolute_points = [a, b, c]
1129        vertices = [[0,1,2]]
1130       
1131        geo = Geo_reference(56,67,-56)
1132
1133        relative_points = geo.change_points_geo_ref(absolute_points)
1134
1135        #print 'Relative', relative_points
1136        #print 'Absolute', absolute_points       
1137
1138        mesh = Mesh(relative_points, vertices, geo_reference=geo)
1139        boundary_polygon = mesh.get_boundary_polygon()
1140
1141        assert allclose(absolute_points, boundary_polygon)
1142       
1143#-------------------------------------------------------------
1144if __name__ == "__main__":
1145    #suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
1146    suite = unittest.makeSuite(Test_Mesh,'test')
1147    runner = unittest.TextTestRunner()
1148    runner.run(suite)
1149
1150
1151
1152
Note: See TracBrowser for help on using the repository browser.