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

Last change on this file since 4836 was 4836, checked in by ole, 16 years ago

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