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

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

Implemeted and tested ticket:175 - flow_through_cross_section.
Also updated manual

File size: 54.1 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    def test_get_triangle_neighbours(self):
1234        a = [0.0, 0.0]
1235        b = [0.0, 2.0]
1236        c = [2.0,0.0]
1237        e = [2.0, 2.0]
1238        points = [a, b, c, e]
1239        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
1240        mesh = Mesh(points, vertices)
1241
1242        neighbours = mesh.get_triangle_neighbours(0)
1243        assert allclose(neighbours, [-1,1,-1])
1244        neighbours = mesh.get_triangle_neighbours(-10)
1245        assert neighbours == []
1246        neighbours = mesh.get_triangle_neighbours(2)
1247        assert neighbours == []
1248
1249
1250    def test_get_intersecting_segments1(self):
1251        """test_get_intersecting_segments(self):
1252
1253        Very simple test (horizontal lines)
1254        """
1255
1256        # Build test mesh
1257       
1258        # Create basic mesh
1259        # 9 points at (0,0), (0, 1), ..., (2,2)
1260        # 8 triangles enumerated from left bottom to right top.
1261        points, vertices, boundary = rectangular(2, 2, 2, 2)
1262        mesh = Mesh(points, vertices, boundary)
1263
1264        # Very simple horizontal line intersecting
1265        #
1266
1267
1268        for y_line in [0.1, 0.2, 0.314159, 0.41, 0.6, 0.99, 1.01, 1.5, 1.77, 1.9]:
1269            if y_line < 1:
1270                ceiling = 1
1271                floor = 0
1272                intersected_triangles = [0,1,4,5]
1273            elif y_line > 1:
1274                ceiling = 2
1275                floor = 1
1276                intersected_triangles = [2,3,6,7]
1277            else:
1278                raise Exception, 'this test is not for parallel lines'
1279
1280
1281            line = [[-1,y_line], [3,y_line]]
1282
1283            L = mesh.get_intersecting_segments(line)
1284            assert len(L) == 4
1285
1286           
1287
1288            # Check all normals point straight down etc
1289            total_length = 0
1290            for x in L:
1291                if x.triangle_id % 2 == 0:
1292                    assert allclose(x.length, ceiling-y_line)
1293                else:
1294                    assert allclose(x.length, y_line-floor)               
1295
1296               
1297                assert allclose(x.normal, [0,-1])
1298
1299                assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
1300                assert allclose(x.segment[0][1], y_line)
1301                assert allclose(x.segment[1][1], y_line)               
1302
1303                assert x.triangle_id in intersected_triangles
1304
1305                total_length += x.length
1306
1307            msg = 'Segments do not add up'
1308            assert allclose(total_length, 2), msg
1309           
1310
1311    def test_get_intersecting_segments_coinciding(self):
1312        """test_get_intersecting_segments_coinciding(self):
1313
1314        Test that lines coinciding with triangle edges work.
1315        """
1316
1317        # Build test mesh
1318       
1319        # Create basic mesh
1320        # 9 points at (0,0), (0, 1), ..., (2,2)
1321        # 8 triangles enumerated from left bottom to right top.
1322        points, vertices, boundary = rectangular(2, 2, 2, 2)
1323        mesh = Mesh(points, vertices, boundary)
1324        intersected_triangles = [1,5]
1325
1326        # Very simple horizontal line intersecting
1327        #
1328
1329        y_line = 1.0
1330       
1331        line = [[-1,y_line], [3,y_line]]
1332
1333        L = mesh.get_intersecting_segments(line)
1334
1335
1336        msg = 'Only two triangles should be returned'   
1337        assert len(L) == 2, msg   
1338           
1339
1340        # Check all
1341        total_length = 0
1342        for x in L:
1343            assert allclose(x.length, 1.0)
1344            assert allclose(x.normal, [0,-1])
1345
1346            assert allclose(x.segment[1][0], x.segment[0][0] + x.length)
1347            assert allclose(x.segment[0][1], y_line)
1348            assert allclose(x.segment[1][1], y_line)                           
1349
1350
1351
1352            assert x.triangle_id in intersected_triangles
1353
1354            total_length += x.length
1355
1356        msg = 'Segments do not add up'
1357        assert allclose(total_length, 2), msg
1358       
1359
1360    def test_get_intersecting_segments_partially_coinciding(self):
1361        """test_get_intersecting_segments_partially_coinciding(self):
1362
1363        Test that line coinciding with triangle edges work.
1364        But this ones only coincide with parts of the edge.
1365        """
1366
1367        # Build test mesh
1368       
1369        # Create basic mesh
1370        # 9 points at (0,0), (0, 1), ..., (2,2)
1371        # 8 triangles enumerated from left bottom to right top.
1372        points, vertices, boundary = rectangular(2, 2, 2, 2)
1373        mesh = Mesh(points, vertices, boundary)
1374        intersected_triangles = [1,5]
1375
1376        # Horizontal line intersecting along center but stopping
1377        # parway through second triangle's edge
1378        #
1379
1380        y_line = 1.0
1381       
1382        #line = [[0, y_line], [2, y_line]]
1383        line = [[0, y_line], [1.5, y_line]]
1384
1385        L = mesh.get_intersecting_segments(line)
1386        #for x in L:
1387        #    print x
1388
1389        msg = 'Two triangles should be returned'   
1390        assert len(L) == 2, msg   
1391           
1392
1393        # Check all
1394        total_length = 0
1395        for x in L:
1396            if x.triangle_id == 1:
1397                assert allclose(x.length, 1)       
1398                assert allclose(x.normal, [0, -1])
1399               
1400            if x.triangle_id == 5:
1401                assert allclose(x.length, 0.5)
1402                assert allclose(x.normal, [0, -1])
1403
1404
1405            assert x.triangle_id in intersected_triangles
1406
1407            total_length += x.length
1408
1409        msg = 'Segments do not add up'
1410        assert allclose(total_length, 1.5), msg           
1411
1412
1413
1414    def test_get_intersecting_segments2(self):
1415        """test_get_intersecting_segments(self):
1416
1417        Lines with a slope
1418        """
1419
1420        s2 = sqrt(2.0)/2
1421       
1422
1423        # Build test mesh
1424       
1425        # Create basic mesh
1426        # 9 points at (0,0), (0, 1), ..., (2,2)
1427        # 8 triangles enumerated from left bottom to right top.
1428        points, vertices, boundary = rectangular(2, 2, 2, 2)
1429        mesh = Mesh(points, vertices, boundary)
1430       
1431
1432        # Diagonal cutting through a vertex and hypothenuses
1433        line = [[0, 2], [2, 0]]
1434        intersected_triangles = [3,2,5,4]               
1435
1436        L = mesh.get_intersecting_segments(line)
1437        assert len(L) == 4
1438
1439        #print L
1440       
1441        # Check all segments
1442        total_length = 0
1443        for i, x in enumerate(L): 
1444            assert allclose(x.length, s2)
1445            assert allclose(x.normal, [-s2, -s2])
1446            assert allclose(sum(x.normal**2), 1)
1447           
1448            assert x.triangle_id in intersected_triangles
1449
1450            total_length += x.length
1451
1452        msg = 'Segments do not add up'
1453        assert allclose(total_length, 4*s2), msg
1454
1455
1456        # Diagonal cutting through a vertex and hypothenuses (reversed)
1457        line = [[2, 0], [0, 2]]
1458        intersected_triangles = [3,2,5,4]               
1459
1460        L = mesh.get_intersecting_segments(line)
1461        assert len(L) == 4
1462
1463        #print L
1464       
1465        # Check all segments
1466        total_length = 0
1467        for i, x in enumerate(L): 
1468            assert allclose(x.length, s2)
1469            assert allclose(x.normal, [s2, s2])
1470            assert allclose(sum(x.normal**2), 1)
1471           
1472            assert x.triangle_id in intersected_triangles
1473
1474            total_length += x.length
1475
1476        msg = 'Segments do not add up'
1477        assert allclose(total_length, 4*s2), msg                   
1478
1479
1480
1481        # Diagonal coinciding with hypothenuses
1482        line = [[2, 2], [0, 0]]
1483        intersected_triangles = [6,0]               
1484
1485        L = mesh.get_intersecting_segments(line)
1486        assert len(L) == 2
1487
1488        #print L
1489       
1490        # Check all segments
1491        total_length = 0
1492        for i, x in enumerate(L): 
1493            assert allclose(x.length, 2*s2)
1494            assert allclose(x.normal, [-s2, s2])
1495            assert allclose(sum(x.normal**2), 1)
1496           
1497            assert x.triangle_id in intersected_triangles
1498
1499            total_length += x.length
1500
1501        msg = 'Segments do not add up'
1502        assert allclose(total_length, 4*s2), msg                       
1503
1504
1505        # Diagonal coinciding with hypothenuses (reversed)
1506        line = [[0, 0], [2, 2]]
1507        intersected_triangles = [6,0]               
1508
1509        L = mesh.get_intersecting_segments(line)
1510        assert len(L) == 2
1511
1512        #print L
1513       
1514        # Check all segments
1515        total_length = 0
1516        for i, x in enumerate(L): 
1517            assert allclose(x.length, 2*s2)
1518            assert allclose(x.normal, [s2, -s2])
1519            assert allclose(sum(x.normal**2), 1)
1520           
1521            assert x.triangle_id in intersected_triangles
1522
1523            total_length += x.length
1524
1525        msg = 'Segments do not add up'
1526        assert allclose(total_length, 4*s2), msg                       
1527
1528
1529
1530        # line with slope [1, -1] cutting through vertices of tri 7 and 6
1531        line = [[1, 2], [2, 1]]
1532        intersected_triangles = [7,6]               
1533
1534        L = mesh.get_intersecting_segments(line)
1535        assert len(L) == 2
1536
1537        #print L
1538       
1539        # Check all segments
1540        total_length = 0
1541        for i, x in enumerate(L): 
1542            assert allclose(x.length, s2)
1543            assert allclose(x.normal, [-s2, -s2])
1544            assert allclose(sum(x.normal**2), 1)
1545           
1546            assert x.triangle_id in intersected_triangles
1547
1548            total_length += x.length
1549
1550        msg = 'Segments do not add up'
1551        assert allclose(total_length, 2*s2), msg
1552
1553
1554        # Arbitrary line with slope [1, -1] cutting through tri 7 and 6
1555        line = [[1.1, 2], [2.1, 1]]
1556        intersected_triangles = [7,6]               
1557
1558        L = mesh.get_intersecting_segments(line)
1559        assert len(L) == 2
1560       
1561        # Check all segments
1562        total_length = 0
1563        for i, x in enumerate(L): 
1564            assert allclose(x.normal, [-s2, -s2])
1565            assert allclose(sum(x.normal**2), 1)
1566
1567            msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles)
1568            assert x.triangle_id in intersected_triangles, msg
1569           
1570
1571
1572    def test_get_intersecting_segments3(self):
1573        """test_get_intersecting_segments(self):
1574
1575        Check that line can stop inside a triangle
1576       
1577        """
1578
1579
1580
1581        s2 = sqrt(2.0)/2
1582       
1583
1584        # Build test mesh
1585       
1586        # Create basic mesh
1587        # 9 points at (0,0), (0, 1), ..., (2,2)
1588        # 8 triangles enumerated from left bottom to right top.
1589        points, vertices, boundary = rectangular(2, 2, 2, 2)
1590        mesh = Mesh(points, vertices, boundary)
1591       
1592
1593        # Line cutting through one triangle and ending on its edge
1594        line = [[0.5, 3], [0.5, 1.5]]
1595        intersected_triangles = [3]               
1596
1597        L = mesh.get_intersecting_segments(line)
1598        assert len(L) == 1
1599        assert L[0].triangle_id == 3
1600        assert allclose(L[0].length, 0.5)       
1601        assert allclose(L[0].normal, [-1,0])               
1602
1603
1604
1605        # Now try to shorten it so that its endpoint falls short of the far edge
1606        line = [[0.5, 3], [0.5, 1.6]]
1607        intersected_triangles = [3]               
1608
1609        L = mesh.get_intersecting_segments(line)
1610        assert len(L) == 1
1611        assert L[0].triangle_id == 3
1612        assert allclose(L[0].length, 0.4)
1613        assert allclose(L[0].normal, [-1,0])
1614
1615        intersected_triangles = [3]
1616
1617        # Now the same, but with direction changed
1618        line = [[0.5, 3], [0.5, 1.6]]
1619        line = [[0.5, 1.6], [0.5, 3]]       
1620        intersected_triangles = [3]               
1621
1622        L = mesh.get_intersecting_segments(line)
1623        assert len(L) == 1
1624        assert L[0].triangle_id == 3
1625        assert allclose(L[0].length, 0.4)
1626        assert allclose(L[0].normal, [1,0])               
1627       
1628
1629           
1630
1631    def test_get_intersecting_segments4(self):
1632        """test_get_intersecting_segments(self):
1633
1634        Simple poly line
1635       
1636        """
1637
1638
1639
1640        s2 = sqrt(2.0)/2
1641       
1642
1643        # Build test mesh
1644       
1645        # Create basic mesh
1646        # 9 points at (0,0), (0, 1), ..., (2,2)
1647        # 8 triangles enumerated from left bottom to right top.
1648        points, vertices, boundary = rectangular(2, 2, 2, 2)
1649        mesh = Mesh(points, vertices, boundary)
1650       
1651
1652        # Polyline with three segments cutting through mesh
1653        line = [[0.5, 3], [0.5, 1.5], [1,1]]
1654
1655        L = mesh.get_intersecting_segments(line)
1656        assert len(L) == 2
1657
1658        for x in L:
1659            if x.triangle_id == 3:
1660                assert allclose(x.length, 0.5)       
1661                assert allclose(x.normal, [-1,0])
1662               
1663            if x.triangle_id == 2:
1664                assert allclose(x.length, s2)
1665                assert allclose(x.normal, [-s2,-s2])
1666
1667
1668
1669    def test_get_intersecting_segments5(self):
1670        """test_get_intersecting_segments(self):
1671
1672        More complex poly line
1673       
1674        """
1675
1676
1677
1678        s2 = sqrt(2.0)/2
1679       
1680
1681        # Build test mesh
1682       
1683        # Create basic mesh
1684        # 9 points at (0,0), (0, 1), ..., (2,2)
1685        # 8 triangles enumerated from left bottom to right top.
1686        points, vertices, boundary = rectangular(2, 2, 2, 2)
1687        mesh = Mesh(points, vertices, boundary)
1688       
1689
1690        # Polyline with three segments cutting through mesh
1691        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75]] 
1692
1693        L = mesh.get_intersecting_segments(line)
1694        assert len(L) == 3
1695
1696        for x in L:
1697            if x.triangle_id == 3:
1698                assert allclose(x.length, 0.5)       
1699                assert allclose(x.normal, [-1,0])
1700               
1701            if x.triangle_id == 2:
1702                msg = str(x.length)
1703                assert allclose(x.length, s2), msg
1704                assert allclose(x.normal, [-s2,-s2])
1705
1706            if x.triangle_id == 5:
1707                segvec = array([line[2][0]-1,
1708                                line[2][1]-1])
1709                msg = str(x.length)
1710                assert allclose(x.length, sqrt(sum(segvec**2))), msg
1711                assert allclose(x.normal, [-s2,-s2])                                               
1712
1713
1714    def test_get_intersecting_segments6(self):
1715        """test_get_intersecting_segments(self):
1716
1717        Even more complex poly line, where line breaks within triangle 5
1718
1719        5 segments are returned even though only four triangles [3,2,5,6] are touched.
1720        Triangle 5 therefor has two segments in it.
1721       
1722        """
1723
1724
1725
1726        s2 = sqrt(2.0)/2
1727       
1728
1729        # Build test mesh
1730       
1731        # Create basic mesh
1732        # 9 points at (0,0), (0, 1), ..., (2,2)
1733        # 8 triangles enumerated from left bottom to right top.
1734        points, vertices, boundary = rectangular(2, 2, 2, 2)
1735        mesh = Mesh(points, vertices, boundary)
1736       
1737
1738        # Polyline with three segments cutting through mesh
1739        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75], [2.25, 1.75]]
1740
1741        L = mesh.get_intersecting_segments(line)
1742        #for x in L:
1743        #    print x
1744
1745        assert len(L) == 5
1746
1747        for x in L:
1748            if x.triangle_id == 3:
1749                assert allclose(x.length, 0.5)       
1750                assert allclose(x.normal, [-1,0])
1751               
1752            if x.triangle_id == 2:
1753                msg = str(x.length)
1754                assert allclose(x.length, s2), msg
1755                assert allclose(x.normal, [-s2,-s2])
1756
1757            if x.triangle_id == 5:
1758                if x.segment == ((1.0, 1.0), (1.25, 0.75)):                   
1759                    segvec = array([line[2][0]-1,
1760                                    line[2][1]-1])
1761                    msg = str(x.length)
1762                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
1763                    assert allclose(x.normal, [-s2,-s2])
1764                elif x.segment == ((1.25, 0.75), (1.5, 1.0)):
1765                    segvec = array([1.5-line[2][0],
1766                                    1.0-line[2][1]])
1767                   
1768                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
1769                    assert allclose(x.normal, [s2,-s2])
1770                else:
1771                    msg = 'Unknow segment: %s' %x.segment
1772                    raise Exception, msg
1773               
1774
1775                   
1776            if x.triangle_id == 6:
1777                assert allclose(x.normal, [s2,-s2])
1778                assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
1779
1780
1781      # Internal test that sum of line segments add up
1782      # to length of input line
1783      #
1784      # Could be useful perhaps
1785      #
1786      #xi1 = line[1][0]
1787      #eta1 = line[1][1]
1788      #linevector = array([xi1-xi0, eta1-eta0])
1789      #linelength = sqrt(sum(linevector**2))
1790      #
1791      #segmentlength = 0
1792      #for segment in triangle_intersections:
1793      #    vector = array([segment[1][0] - segment[0][0],
1794      #                    segment[1][1] - segment[0][1]])
1795      #    length = sqrt(sum(vector**2))     
1796      #    segmentlength += length
1797      #
1798      #msg = 'Sum of intersecting segments do not add up'   
1799      #assert allclose(segmentlength, linelength), msg   
1800
1801
1802
1803
1804    def test_get_intersecting_segments7(self):
1805        """test_get_intersecting_segments(self):
1806
1807        Check that line can stop inside a triangle - this is from
1808        flow throug a cross sections example in test_datamanager.
1809       
1810        """
1811
1812        # Build test mesh
1813        width = 5
1814        length = 100
1815        t_end = 1
1816        points, vertices, boundary = rectangular(length, width,
1817                                                 length, width)
1818
1819        mesh = Mesh(points, vertices, boundary)
1820       
1821
1822        # A range of partial lines
1823        x = length/2.
1824        for i in range(10):
1825            start_point = [length/2., i*width/10.]
1826            #print
1827            #print start_point
1828                           
1829            line = [start_point, [length/2., width]]
1830 
1831            L = mesh.get_intersecting_segments(line)
1832
1833            if start_point[1] < 1:
1834                assert len(L) == 5
1835               
1836           
1837            total_length = 0   
1838            for x in L:
1839                total_length += x.length
1840               
1841
1842            ref_length = line[1][1] - line[0][1]
1843            #print ref_length, total_length
1844            assert allclose(total_length, ref_length)
1845
1846
1847
1848       
1849#-------------------------------------------------------------
1850if __name__ == "__main__":
1851    #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
1852    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
1853    #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
1854    suite = unittest.makeSuite(Test_Mesh,'test')
1855    runner = unittest.TextTestRunner()
1856    runner.run(suite)
1857
1858
1859
1860
Note: See TracBrowser for help on using the repository browser.