source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py @ 8124

Last change on this file since 8124 was 8124, checked in by wilsonr, 13 years ago

Made changes described in ticket 359.

File size: 54.0 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 mesh_factory import rectangular_periodic
15from anuga.config import epsilon
16
17from anuga.coordinate_transforms.geo_reference import Geo_reference
18from anuga.geometry.polygon import is_inside_polygon
19from anuga.utilities.numerical_tools import ensure_numeric
20
21import numpy as num
22
23
24def distance(x, y):
25    return sqrt(num.sum((num.array(x)-num.array(y))**2))
26
27class Test_Mesh(unittest.TestCase):
28    def setUp(self):
29        pass
30
31    def tearDown(self):
32        pass
33
34    def test_triangle_inputs(self):
35        points = [[0.0, 0.0], [4.0, 0.0], [0.0, 3.0]]
36        vertices = [0,1,2] #Wrong
37
38        try:
39            mesh = Mesh(points, vertices)
40        except:
41            pass
42        else:
43            msg = 'Should have raised exception'
44            raise Exception(msg)
45
46
47    def test_basic_triangle(self):
48
49        a = [0.0, 0.0]
50        b = [4.0, 0.0]
51        c = [0.0, 3.0]
52
53        points = [a, b, c]
54        vertices = [[0,1,2]]
55        mesh = Mesh(points, vertices)
56
57        #Centroid
58        centroid = mesh.centroid_coordinates[0]
59        assert centroid[0] == 4.0/3
60        assert centroid[1] == 1.0
61
62        #Area
63        assert mesh.areas[0] == 6.0,\
64               'Area was %f, should have been 6.0' %mesh.areas[0]
65
66        #Normals
67        normals = mesh.get_normals()
68        assert num.allclose(normals[0, 0:2], [3.0/5, 4.0/5])
69        assert num.allclose(normals[0, 2:4], [-1.0, 0.0])
70        assert num.allclose(normals[0, 4:6], [0.0, -1.0])
71
72        assert num.allclose(mesh.get_normal(0,0), [3.0/5, 4.0/5])
73        assert num.allclose(mesh.get_normal(0,1), [-1.0, 0.0])
74        assert num.allclose(mesh.get_normal(0,2), [0.0, -1.0])
75
76        #Edge lengths
77        assert num.allclose(mesh.edgelengths[0], [5.0, 3.0, 4.0])
78
79
80        #Vertex coordinates
81        #V = mesh.get_vertex_coordinates()
82        #assert allclose(V[0], [0.0, 0.0, 4.0, 0.0, 0.0, 3.0])
83       
84
85        V = mesh.get_vertex_coordinates()
86        assert num.allclose(V, [ [0.0, 0.0],
87                             [4.0, 0.0],
88                             [0.0, 3.0] ])
89
90        V0 = mesh.get_vertex_coordinate(0, 0)
91        assert num.allclose(V0, [0.0, 0.0])
92
93        V1 = mesh.get_vertex_coordinate(0, 1)
94        assert num.allclose(V1, [4.0, 0.0])
95
96        V2 = mesh.get_vertex_coordinate(0, 2)
97        assert num.allclose(V2, [0.0, 3.0])
98
99
100        #General tests:
101
102        #Test that points are arranged in a counter clock wise order etc
103        mesh.check_integrity()
104
105
106        #Test that the centroid is located 2/3 of the way
107        #from each vertex to the midpoint of the opposite side
108
109        V = mesh.get_vertex_coordinates()
110        x0 = V[0, 0]; y0 = V[0, 1]
111        x1 = V[1, 0]; y1 = V[1, 1]
112        x2 = V[2, 0]; y2 = V[2, 1]
113        #x0 = V[0,0]
114        #y0 = V[0,1]
115        #x1 = V[0,2]
116        #y1 = V[0,3]
117        #x2 = V[0,4]
118        #y2 = V[0,5]
119
120        m0 = [(x1 + x2)/2, (y1 + y2)/2]
121        m1 = [(x0 + x2)/2, (y0 + y2)/2]
122        m2 = [(x1 + x0)/2, (y1 + y0)/2]
123
124        d0 = distance(centroid, [x0, y0])
125        d1 = distance(m0, [x0, y0])
126        assert d0 == 2*d1/3
127        #
128        d0 = distance(centroid, [x1, y1])
129        d1 = distance(m1, [x1, y1])
130        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
131
132        d0 = distance(centroid, [x2, y2])
133        d1 = distance(m2, [x2, y2])
134        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
135
136        #Radius
137        d0 = distance(centroid, m0)
138        assert d0 == 5.0/6
139
140        d1 = distance(centroid, m1)
141        assert d1 == sqrt(73.0/36)
142
143        d2 = distance(centroid, m2)
144        assert d2 == sqrt(13.0/9)
145
146        assert mesh.radii[0] == min(d0, d1, d2)
147        assert mesh.radii[0] == 5.0/6
148
149
150        #Let x be the centroid of triangle abc.
151        #Test that areas of the three triangles axc, cxb, and bxa are equal.
152        points = [a, b, c, centroid]
153        vertices = [[0,3,2], [2,3,1], [1,3,0]]
154        new_mesh = Mesh(points, vertices)
155
156        assert new_mesh.areas[0] == new_mesh.areas[1]
157        assert new_mesh.areas[1] == new_mesh.areas[2]
158        assert new_mesh.areas[1] == new_mesh.areas[2]
159
160        assert new_mesh.areas[1] == mesh.areas[0]/3
161
162
163
164    def test_general_triangle(self):
165        a = [2.0, 1.0]
166        b = [6.0, 2.0]
167        c = [1.0, 3.0]
168
169        points = [a, b, c]
170        vertices = [[0,1,2]]
171
172        mesh = Mesh(points, vertices)
173        centroid = mesh.centroid_coordinates[0]
174
175
176        #Test that the centroid is located 2/3 of the way
177        #from each vertex to the midpoint of the opposite side
178
179        V = mesh.get_vertex_coordinates()
180        x0 = V[0, 0]; y0 = V[0, 1]
181        x1 = V[1, 0]; y1 = V[1, 1]
182        x2 = V[2, 0]; y2 = V[2, 1]       
183
184        #x0 = V[0,0]
185        #y0 = V[0,1]
186        #x1 = V[0,2]
187        #y1 = V[0,3]
188        #x2 = V[0,4]
189        #y2 = V[0,5]
190
191        m0 = [(x1 + x2)/2, (y1 + y2)/2]
192        m1 = [(x0 + x2)/2, (y0 + y2)/2]
193        m2 = [(x1 + x0)/2, (y1 + y0)/2]
194
195        d0 = distance(centroid, [x0, y0])
196        d1 = distance(m0, [x0, y0])
197        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
198        #
199        d0 = distance(centroid, [x1, y1])
200        d1 = distance(m1, [x1, y1])
201        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
202
203        d0 = distance(centroid, [x2, y2])
204        d1 = distance(m2, [x2, y2])
205        assert abs(d0 - 2*d1/3) < epsilon, '%e, %e' %(d0, 2*d1/3)
206
207        #Radius
208        d0 = distance(centroid, m0)
209        d1 = distance(centroid, m1)
210        d2 = distance(centroid, m2)
211        assert mesh.radii[0] == min(d0, d1, d2)
212
213
214
215        #Let x be the centroid of triangle abc.
216        #Test that areas of the three triangles axc, cxb, and bxa are equal.
217
218        points = [a, b, c, centroid]
219        vertices = [[0,3,2], [2,3,1], [1,3,0]]
220        new_mesh = Mesh(points, vertices)
221
222        assert new_mesh.areas[0] == new_mesh.areas[1]
223        assert new_mesh.areas[1] == new_mesh.areas[2]
224        assert new_mesh.areas[1] == new_mesh.areas[2]
225
226        assert new_mesh.areas[1] == mesh.areas[0]/3
227
228
229        #Test that points are arranged in a counter clock wise order
230        mesh.check_integrity()
231
232    def test_inscribed_circle_equilateral(self):
233        """test that the radius is calculated correctly by mesh in the case of an equilateral triangle"""
234        a = [0.0, 0.0]
235        b = [2.0, 0.0]
236        c = [1.0, sqrt(3.0)]
237
238        points = [a, b, c]
239        vertices = [[0,1,2]]
240
241        mesh = Mesh(points, vertices,use_inscribed_circle=False)
242        assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'Steve''s doesn''t work'
243
244        mesh = Mesh(points, vertices,use_inscribed_circle=True)
245        assert num.allclose(mesh.radii[0],sqrt(3.0)/3),'inscribed circle doesn''t work'
246
247    def test_inscribed_circle_rightangle_triangle(self):
248        """test that the radius is calculated correctly by mesh in the case of a right-angled triangle"""
249        a = [0.0, 0.0]
250        b = [4.0, 0.0]
251        c = [0.0, 3.0]
252
253        points = [a, b, c]
254        vertices = [[0,1,2]]
255
256        mesh = Mesh(points, vertices,use_inscribed_circle=False)
257        assert num.allclose(mesh.radii[0],5.0/6),'Steve''s doesn''t work'
258
259        mesh = Mesh(points, vertices,use_inscribed_circle=True)
260        assert num.allclose(mesh.radii[0],1.0),'inscribed circle doesn''t work'
261
262
263    def test_two_triangles(self):
264        a = [0.0, 0.0]
265        b = [0.0, 2.0]
266        c = [2.0,0.0]
267        e = [2.0, 2.0]
268        points = [a, b, c, e]
269        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
270        mesh = Mesh(points, vertices)
271
272        assert mesh.areas[0] == 2.0
273
274        assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
275
276
277        #Test that points are arranged in a counter clock wise order
278        mesh.check_integrity()
279
280
281
282    def test_more_triangles(self):
283
284        a = [0.0, 0.0]
285        b = [0.0, 2.0]
286        c = [2.0, 0.0]
287        d = [0.0, 4.0]
288        e = [2.0, 2.0]
289        f = [4.0, 0.0]
290
291        points = [a, b, c, d, e, f]
292        #bac, bce, ecf, dbe, daf, dae
293        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
294        mesh = Mesh(points, vertices)
295
296        #Test that points are arranged in a counter clock wise order
297        mesh.check_integrity()
298
299        assert mesh.areas[0] == 2.0
300        assert mesh.areas[1] == 2.0
301        assert mesh.areas[2] == 2.0
302        assert mesh.areas[3] == 2.0
303
304        assert mesh.edgelengths[1,0] == 2.0
305        assert mesh.edgelengths[1,1] == 2.0
306        assert mesh.edgelengths[1,2] == sqrt(8.0)
307
308        assert num.allclose(mesh.centroid_coordinates[0], [2.0/3, 2.0/3])
309        assert num.allclose(mesh.centroid_coordinates[1], [4.0/3, 4.0/3])
310        assert num.allclose(mesh.centroid_coordinates[2], [8.0/3, 2.0/3])
311        assert num.allclose(mesh.centroid_coordinates[3], [2.0/3, 8.0/3])
312
313    def test_mesh_and_neighbours(self):
314        a = [0.0, 0.0]
315        b = [0.0, 2.0]
316        c = [2.0,0.0]
317        d = [0.0, 4.0]
318        e = [2.0, 2.0]
319        f = [4.0,0.0]
320
321
322        points = [a, b, c, d, e, f]
323
324        #bac, bce, ecf, dbe
325        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
326        mesh = Mesh(points, vertices)
327
328        mesh.check_integrity()
329
330
331        T = mesh
332        tid = 0
333        assert T.number_of_boundaries[tid] == 2
334        assert T.neighbours[tid, 0] < 0  #Opposite point b (0,2)
335        assert T.neighbours[tid, 1] == 1 #Opposite point a (0,0)
336        assert T.neighbours[tid, 2] < 0  #Opposite point c (2,0)
337
338        tid = 1
339        assert T.number_of_boundaries[tid] == 0
340        assert T.neighbours[tid, 0] == 2  #Opposite point b (0,2)
341        assert T.neighbours[tid, 1] == 3  #Opposite point c (2,0)
342        assert T.neighbours[tid, 2] == 0  #Opposite point e (2,2)
343
344        tid = 2
345        assert T.number_of_boundaries[tid] == 2
346        assert T.neighbours[tid, 0] < 0   #Opposite point e (2,2)
347        assert T.neighbours[tid, 1] < 0   #Opposite point c (2,0)
348        assert T.neighbours[tid, 2] == 1  #Opposite point f (4,0)
349
350        tid = 3
351        assert T.number_of_boundaries[tid] == 2
352        assert T.neighbours[tid, 0] == 1  #Opposite point d (0,4)
353        assert T.neighbours[tid, 1] < 0   #Opposite point b (0,3)
354        assert T.neighbours[tid, 2] < 0   #Opposite point e (2,2)
355
356        #Neighbouring edges
357        tid = 0
358        assert T.neighbour_edges[tid, 0] < 0  #Opposite point b (0,2)
359        assert T.neighbour_edges[tid, 1] == 2 #Opposite point a (0,0)
360        assert T.neighbour_edges[tid, 2] < 0  #Opposite point c (2,0)
361
362        tid = 1
363        assert T.neighbour_edges[tid, 0] == 2 #Opposite point b (0,2)
364        assert T.neighbour_edges[tid, 1] == 0 #Opposite point c (2,0)
365        assert T.neighbour_edges[tid, 2] == 1 #Opposite point e (2,2)
366
367        tid = 2
368        assert T.neighbour_edges[tid, 0] < 0  #Opposite point e (2,2)
369        assert T.neighbour_edges[tid, 1] < 0  #Opposite point c (2,0)
370        assert T.neighbour_edges[tid, 2] == 0 #Opposite point f (4,0)
371
372        tid = 3
373        assert T.neighbour_edges[tid, 0] == 1 #Opposite point d (0,4)
374        assert T.neighbour_edges[tid, 1] < 0  #Opposite point b (0,3)
375        assert T.neighbour_edges[tid, 2] < 0  #Opposite point e (2,2)
376
377
378    def test_build_neighbour_structure_duplicates(self):
379        p0 = [-66.0, 14.0]
380        p1 = [14.0, -66.0]
381        p2 = [14.0, 14.0]
382        p3 = [60.0, 20.0]
383        p4 = [10.0, 60.0]
384        p5 = [60.0, 60.0]
385
386        points = [p0, p1, p2, p3, p4, p5]
387        triangles = [ [0, 1, 2],
388                      [3, 2, 1],
389                      [0, 2, 4],
390                      [0, 2, 4],
391                      [4, 2, 5],
392                      [5, 2, 3]]
393        try:
394            mesh = Mesh(points, triangles)
395        except:
396            pass
397        else:
398            raise Exception("triangle edge duplicates not caught")
399
400
401    def test_rectangular_mesh_basic(self):
402        M=1
403        N=1
404
405        points, vertices, boundary = rectangular(M, N)
406        mesh = Mesh(points, vertices, boundary)
407
408        #Test that points are arranged in a counter clock wise order
409        mesh.check_integrity()
410
411        M=2
412        N=2
413        points, vertices, boundary = rectangular(M, N)
414        mesh = Mesh(points, vertices, boundary)
415
416        #Test that points are arranged in a counter clock wise order
417        mesh.check_integrity()
418
419        #assert mesh.boundary[(7,1)] == 2 # top
420        assert mesh.boundary[(7,1)] == 'top' # top
421        assert mesh.boundary[(3,1)] == 'top' # top
422
423
424
425
426
427
428    def test_boundary_tags(self):
429
430
431        points, vertices, boundary = rectangular(4, 4)
432        mesh = Mesh(points, vertices, boundary)
433
434
435        #Test that points are arranged in a counter clock wise order
436        mesh.check_integrity()
437
438        #print mesh.get_boundary_tags()
439        #print mesh.boundary
440
441        for k in [1, 3, 5, 7]:
442            assert mesh.boundary[(k,2)] == 'left'
443
444        for k in [24, 26, 28, 30]:
445            assert mesh.boundary[(k,2)] == 'right'
446
447        for k in [7, 15, 23, 31]:
448            assert mesh.boundary[(k,1)] == 'top'
449        for k in [0, 8, 16, 24]:
450            assert mesh.boundary[(k,1)] == 'bottom'
451
452
453
454    def test_rectangular_mesh(self):
455        M=4
456        N=16
457        len1 = 100.0
458        len2 = 17.0
459
460        points, vertices, boundary = rectangular(M, N, len1, len2)
461        mesh = Mesh(points, vertices, boundary)
462
463        assert len(mesh) == 2*M*N
464
465        for i in range(len(mesh)):
466            assert mesh.areas[i] == len1*len2/(2*M*N)
467
468            hypo = sqrt((len1/M)**2 + (len2/N)**2) #hypothenuse
469            assert mesh.edgelengths[i, 0] == hypo
470            assert mesh.edgelengths[i, 1] == len1/M #x direction
471            assert mesh.edgelengths[i, 2] == len2/N #y direction
472
473        #Test that points are arranged in a counter clock wise order
474        mesh.check_integrity()
475
476
477    def test_rectangular_mesh2(self):
478        #Check that integers don't cause trouble
479        N = 16
480
481        points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
482        mesh = Mesh(points, vertices, boundary)
483
484
485    def test_surrogate_neighbours(self):
486        a = [0.0, 0.0]
487        b = [0.0, 2.0]
488        c = [2.0,0.0]
489        d = [0.0, 4.0]
490        e = [2.0, 2.0]
491        f = [4.0,0.0]
492
493        points = [a, b, c, d, e, f]
494
495        #bac, bce, ecf, dbe
496        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
497        mesh = Mesh(points, vertices)
498        mesh.check_integrity()
499
500
501        T = mesh
502        tid = 0
503        assert T.number_of_boundaries[tid] == 2
504        assert T.surrogate_neighbours[tid, 0] == tid
505        assert T.surrogate_neighbours[tid, 1] == 1
506        assert T.surrogate_neighbours[tid, 2] == tid
507
508        tid = 1
509        assert T.number_of_boundaries[tid] == 0
510        assert T.surrogate_neighbours[tid, 0] == 2
511        assert T.surrogate_neighbours[tid, 1] == 3
512        assert T.surrogate_neighbours[tid, 2] == 0
513
514        tid = 2
515        assert T.number_of_boundaries[tid] == 2
516        assert T.surrogate_neighbours[tid, 0] == tid
517        assert T.surrogate_neighbours[tid, 1] == tid
518        assert T.surrogate_neighbours[tid, 2] == 1
519
520        tid = 3
521        assert T.number_of_boundaries[tid] == 2
522        assert T.surrogate_neighbours[tid, 0] == 1
523        assert T.surrogate_neighbours[tid, 1] == tid
524        assert T.surrogate_neighbours[tid, 2] == tid
525
526
527    def test_boundary_inputs(self):
528        a = [0.0, 0.0]
529        b = [0.0, 2.0]
530        c = [2.0,0.0]
531        d = [0.0, 4.0]
532        e = [2.0, 2.0]
533        f = [4.0,0.0]
534
535        points = [a, b, c, d, e, f]
536
537        #bac, bce, ecf, dbe
538        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
539
540        boundary = { (0, 0): 'First',
541                     (0, 2): 'Second',
542                     (2, 0): 'Third',
543                     (2, 1): 'Fourth',
544                     (3, 1): 'Fifth',
545                     (3, 2): 'Sixth'}
546
547
548        mesh = Mesh(points, vertices, boundary)
549        mesh.check_integrity()
550
551
552        #Check enumeration
553        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
554        #    b = -k-1
555        #    assert mesh.neighbours[vol_id, edge_id] == b
556
557
558
559    def test_boundary_inputs_using_one_default(self):
560        a = [0.0, 0.0]
561        b = [0.0, 2.0]
562        c = [2.0,0.0]
563        d = [0.0, 4.0]
564        e = [2.0, 2.0]
565        f = [4.0,0.0]
566
567        points = [a, b, c, d, e, f]
568
569        #bac, bce, ecf, dbe
570        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
571
572        boundary = { (0, 0): 'First',
573                     (0, 2): 'Second',
574                     (2, 0): 'Third',
575                     (2, 1): 'Fourth',
576                     #(3, 1): 'Fifth',  #Skip this
577                     (3, 2): 'Sixth'}
578
579
580        mesh = Mesh(points, vertices, boundary)
581        mesh.check_integrity()
582
583        from anuga.config import default_boundary_tag
584        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
585
586
587        #Check enumeration
588        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
589        #    b = -k-1
590        #    assert mesh.neighbours[vol_id, edge_id] == b
591
592    def test_boundary_inputs_using_all_defaults(self):
593        a = [0.0, 0.0]
594        b = [0.0, 2.0]
595        c = [2.0,0.0]
596        d = [0.0, 4.0]
597        e = [2.0, 2.0]
598        f = [4.0,0.0]
599
600        points = [a, b, c, d, e, f]
601
602        #bac, bce, ecf, dbe
603        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
604
605        boundary = { (0, 0): 'First',
606                     (0, 2): 'Second',
607                     (2, 0): 'Third',
608                     (2, 1): 'Fourth',
609                     #(3, 1): 'Fifth',  #Skip this
610                     (3, 2): 'Sixth'}
611
612
613        mesh = Mesh(points, vertices) #, boundary)
614        mesh.check_integrity()
615
616        from anuga.config import default_boundary_tag
617        assert mesh.boundary[ (0, 0) ] == default_boundary_tag
618        assert mesh.boundary[ (0, 2) ] == default_boundary_tag
619        assert mesh.boundary[ (2, 0) ] == default_boundary_tag
620        assert mesh.boundary[ (2, 1) ] == default_boundary_tag
621        assert mesh.boundary[ (3, 1) ] == default_boundary_tag
622        assert mesh.boundary[ (3, 2) ] == default_boundary_tag
623
624
625        #Check enumeration
626        #for k, (vol_id, edge_id) in enumerate(mesh.boundary_segments):
627        #    b = -k-1
628        #    assert mesh.neighbours[vol_id, edge_id] == b
629
630
631
632
633
634
635    def test_inputs(self):
636        a = [0.0, 0.0]
637        b = [0.0, 2.0]
638        c = [2.0,0.0]
639        d = [0.0, 4.0]
640        e = [2.0, 2.0]
641        f = [4.0,0.0]
642
643        points = [a, b, c, d, e, f]
644
645        #bac, bce, ecf, dbe
646        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
647
648        #Too few points
649        try:
650            mesh = Mesh([points[0]], vertices)
651        except AssertionError:
652            pass
653        else:
654            raise Exception('Should have raised an exception')
655
656        #Too few points - 1 element
657        try:
658            mesh = Mesh([points[0]], [vertices[0]])
659        except AssertionError:
660            pass
661        else:
662            raise Exception('Should have raised an exception')
663
664        #Wrong dimension of vertices
665        try:
666            mesh = Mesh(points, vertices[0])
667        except AssertionError:
668            pass
669        else:
670            raise Exception('Should have raised an exception')
671
672        #Unsubscriptable coordinates object raises exception
673        try:
674            mesh = Mesh(points[0], [vertices[0]])
675        except AssertionError:
676            pass
677        else:
678            raise Exception('Should have raised an exception')
679
680        #FIXME: This has been commented out pending a decision
681        #whether to allow partial boundary tags or not
682        #
683        #Not specifying all boundary tags
684        #try:
685        #    mesh = Mesh(points, vertices, {(3,0): 'x'})
686        #except AssertionError:
687        #    pass
688        #else:
689        #    raise Exception('Should have raised an exception')
690
691        #Specifying wrong non existing segment
692        try:
693            mesh = Mesh(points, vertices, {(5,0): 'x'})
694        except AssertionError:
695            pass
696        else:
697            raise Exception('Should have raised an exception')
698
699
700
701
702    def test_internal_boundaries(self):
703        """
704        get values based on triangle lists.
705        """
706        from mesh_factory import rectangular
707
708        #Create basic mesh
709        points, vertices, boundary = rectangular(1, 3)
710
711        # Add an internal boundary
712        boundary[(2,0)] = 'internal'
713        boundary[(1,0)] = 'internal'
714
715        #Create shallow water domain
716        domain = Mesh(points, vertices, boundary)
717        domain.build_tagged_elements_dictionary({'bottom':[0,1],
718                                                 'top':[4,5],
719                                                 'all':[0,1,2,3,4,5]})
720
721
722    def test_boundary_polygon(self):
723        from mesh_factory import rectangular
724        #from mesh import Mesh
725
726        #Create basic mesh
727        points, vertices, boundary = rectangular(2, 2)
728        mesh = Mesh(points, vertices, boundary)
729
730
731        P = mesh.get_boundary_polygon()
732
733        assert len(P) == 8
734        assert num.allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
735                                [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
736                                [0.0, 1.0], [0.0, 0.5]])
737        for p in points:
738            #print p, P
739            assert is_inside_polygon(p, P)
740
741
742    def test_boundary_polygon_II(self):
743
744        #Points
745        a = [0.0, 0.0] #0
746        b = [0.0, 0.5] #1
747        c = [0.0, 1.0] #2
748        d = [0.5, 0.0] #3
749        e = [0.5, 0.5] #4
750        f = [1.0, 0.0] #5
751        g = [1.0, 0.5] #6
752        h = [1.0, 1.0] #7
753        i = [1.5, 0.5] #8
754
755        points = [a, b, c, d, e, f, g, h, i]
756
757        #dea, bae, bec, fgd,
758        #edg, ghe, gfi, gih
759        vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
760                     [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
761
762        mesh = Mesh(points, vertices)
763
764        mesh.check_integrity()
765
766        P = mesh.get_boundary_polygon()
767
768        assert len(P) == 8
769        assert num.allclose(P, [a, d, f, i, h, e, c, b])
770
771        for p in points:
772            #print p, P
773            assert is_inside_polygon(p, P)
774
775
776    def test_boundary_polygon_III(self):
777        """Same as II but vertices ordered differently
778        """
779
780
781        #Points
782        a = [0.0, 0.0] #0
783        b = [0.0, 0.5] #1
784        c = [0.0, 1.0] #2
785        d = [0.5, 0.0] #3
786        e = [0.5, 0.5] #4
787        f = [1.0, 0.0] #5
788        g = [1.0, 0.5] #6
789        h = [1.0, 1.0] #7
790        i = [1.5, 0.5] #8
791
792        points = [a, b, c, d, e, f, g, h, i]
793
794        #edg, ghe, gfi, gih
795        #dea, bae, bec, fgd,
796        vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
797                    [3,4,0], [1,0,4], [1,4,2], [5,6,3]]
798
799
800        mesh = Mesh(points, vertices)
801        mesh.check_integrity()
802
803
804        P = mesh.get_boundary_polygon()
805
806        assert len(P) == 8
807        assert num.allclose(P, [a, d, f, i, h, e, c, b])
808
809        for p in points:
810            assert is_inside_polygon(p, P)
811
812           
813    def test_boundary_polygon_IIIa(self):
814        """test_boundary_polygon_IIIa - Check pathological situation where
815        one triangle has no neighbours. This may be the case if a mesh
816        is partitioned using pymetis.
817        """
818
819
820        #Points
821        a = [0.0, 0.0] #0
822        b = [0.0, 0.5] #1
823        c = [0.0, 1.0] #2
824        d = [0.5, 0.0] #3
825        e = [0.5, 0.5] #4
826        f = [1.0, 0.0] #5
827        g = [1.0, 0.5] #6
828        h = [1.0, 1.0] #7
829       
830        # Add pathological triangle with no neighbours to an otherwise
831        # trivial mesh
832       
833        points = [a, b, c, d, e, f, g, h]
834
835        #cbe, aeb, dea, fed, ghe (pathological triangle)
836        vertices = [[2,1,4], [0,4,1], [3,4,0], [5,4,3],
837                    [6,7,4]]   
838                   
839        mesh = Mesh(points, vertices)
840        mesh.check_integrity()
841
842        P = mesh.get_boundary_polygon(verbose=False)
843
844       
845        assert len(P) == 9
846       
847        # Note that point e appears twice!
848        assert num.allclose(P, [a, d, f, e, g, h, e, c, b])
849
850        for p in points:
851            msg = 'Point %s is not inside polygon %s'\
852            %(p, P)     
853            assert is_inside_polygon(p, P), msg             
854
855
856                                       
857           
858           
859
860    def test_boundary_polygon_IV(self):
861        """Reproduce test test_spatio_temporal_file_function_time
862        from test_util.py that looked as if it produced the wrong boundary
863        """
864
865        from mesh_factory import rectangular       
866
867        #Create a domain to hold test grid
868        #(0:15, -20:10)
869        points, vertices, boundary =\
870                rectangular(4, 4, 15, 30, origin = (0, -20))       
871
872        #####
873        mesh = Mesh(points, vertices)
874        mesh.check_integrity()
875
876        P = mesh.get_boundary_polygon()
877
878        #print P
879        assert len(P) == 16
880        for p in points:
881            assert is_inside_polygon(p, P)
882
883
884
885        #####
886        mesh = Mesh(points, vertices, boundary)
887        mesh.check_integrity()
888
889        P = mesh.get_boundary_polygon()
890
891       
892        #print P, len(P)
893        assert len(P) == 16
894
895        for p in points:
896            assert is_inside_polygon(p, P)
897
898        #print mesh.statistics()   
899
900
901
902    def test_boundary_polygon_V(self):
903        """Create a discontinuous mesh (duplicate vertices)
904        and check that boundary is as expected
905       
906        """
907
908        #Points
909        a = [0.0, 0.0] #0
910        b = [0.0, 0.5] #1
911        c = [0.0, 1.0] #2
912        d = [0.5, 0.0] #3
913        e = [0.5, 0.5] #4
914        f = [1.0, 0.0] #5
915        g = [1.0, 0.5] #6
916        h = [1.0, 1.0] #7
917        i = [1.5, 0.5] #8
918
919        #Duplicate points for triangles edg [4,3,6] (central) and
920        #gid [6,8,7] (top right boundary) to them disconnected
921        #from the others
922
923        e0 = [0.5, 0.5] #9
924        d0 = [0.5, 0.0] #10
925        g0 = [1.0, 0.5] #11
926        i0 = [1.5, 0.5] #12       
927       
928
929        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
930
931
932
933        #dea, bae, bec, fgd,
934        #edg, ghe, gfi, gih
935        #vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
936        #             [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
937
938
939        #dea, bae, bec, fgd,
940        #e0d0g0, ghe, gfi, g0i0h
941        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
942                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]       
943
944        mesh = Mesh(points, vertices)
945
946        mesh.check_integrity()
947
948        P = mesh.get_boundary_polygon()
949
950        #print P
951       
952        assert len(P) == 8
953        assert num.allclose(P, [a, d, f, i, h, e, c, b])
954        assert num.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)])
955       
956
957        for p in points:
958            #print p, P
959            assert is_inside_polygon(p, P)
960
961
962
963    def test_boundary_polygon_VI(self):
964        """test_boundary_polygon_VI(self)
965
966        Create a discontinuous mesh (duplicate vertices) from a real situation that failed
967        and check that boundary is as expected
968        """
969
970        # First do the continuous version of mesh
971
972        points = [[   6626.85400391,      0.        ],
973                  [      0.        ,  38246.4140625 ],
974                  [   9656.2734375 ,  68351.265625  ],
975                  [  20827.25585938,  77818.203125  ],
976                  [  32755.59375   ,  58126.9765625 ],
977                  [  35406.3359375 ,  79332.9140625 ],
978                  [  31998.23828125,  88799.84375   ],
979                  [  23288.65820313, 104704.296875  ],
980                  [  32187.57617188, 109816.4375    ],
981                  [  50364.08984375, 110763.1328125 ],
982                  [  80468.9453125 ,  96184.0546875 ],
983                  [  86149.1015625 , 129886.34375   ],
984                  [ 118715.359375  , 129886.34375   ],
985                  [ 117768.6640625 ,  85770.4296875 ],
986                  [ 101485.5390625 ,  45251.9453125 ],
987                  [  49985.4140625 ,   2272.06396484],
988                  [  51737.94140625,  90559.2109375 ],
989                  [  56659.0703125 ,  65907.6796875 ],
990                  [  75735.4765625 ,  23762.00585938],
991                  [  52341.70703125,  38563.39453125]]
992
993        triangles = [[19, 0,15],
994                     [ 2, 4, 3],
995                     [ 4, 2, 1],
996                     [ 1,19, 4],
997                     [15,18,19],
998                     [18,14,17],
999                     [19, 1, 0],
1000                     [ 6, 8, 7],
1001                     [ 8, 6,16],
1002                     [10, 9,16],
1003                     [17, 5, 4],
1004                     [16,17,10],
1005                     [17,19,18],
1006                     [ 5,17,16],
1007                     [10,14,13],
1008                     [10,17,14],
1009                     [ 8,16, 9],
1010                     [12,11,10],
1011                     [10,13,12],
1012                     [19,17, 4],
1013                     [16, 6, 5]]
1014
1015        mesh = Mesh(points, triangles)
1016        mesh.check_integrity()
1017        Pref = mesh.get_boundary_polygon()
1018
1019        #plot_polygons([ensure_numeric(Pref)], 'goodP')
1020
1021        for p in points:
1022            assert is_inside_polygon(p, Pref)
1023       
1024       
1025        # Then do the discontinuous version
1026        import warnings
1027        warnings.filterwarnings('ignore')
1028
1029       
1030        points = [[  52341.70703125,  38563.39453125],
1031                  [   6626.85400391,      0.        ],
1032                  [  49985.4140625 ,   2272.06396484],
1033                  [   9656.2734375 ,  68351.265625  ],
1034                  [  32755.59375   ,  58126.9765625 ],
1035                  [  20827.25585938,  77818.203125  ],
1036                  [  32755.59375   ,  58126.9765625 ],
1037                  [   9656.2734375 ,  68351.265625  ],
1038                  [      0.        ,  38246.4140625 ],
1039                  [      0.        ,  38246.4140625 ],
1040                  [  52341.70703125,  38563.39453125],
1041                  [  32755.59375   ,  58126.9765625 ],
1042                  [  49985.4140625 ,   2272.06396484],
1043                  [  75735.4765625 ,  23762.00585938],
1044                  [  52341.70703125,  38563.39453125],
1045                  [  75735.4765625 ,  23762.00585938],
1046                  [ 101485.5390625 ,  45251.9453125 ],
1047                  [  56659.0703125 ,  65907.6796875 ],
1048                  [  52341.70703125,  38563.39453125],
1049                  [      0.        ,  38246.4140625 ],
1050                  [   6626.85400391,      0.        ],
1051                  [  31998.23828125,  88799.84375   ],
1052                  [  32187.57617188, 109816.4375    ],
1053                  [  23288.65820313, 104704.296875  ],
1054                  [  32187.57617188, 109816.4375    ],
1055                  [  31998.23828125,  88799.84375   ],
1056                  [  51737.94140625,  90559.2109375 ],
1057                  [  80468.9453125 ,  96184.0546875 ],
1058                  [  50364.08984375, 110763.1328125 ],
1059                  [  51737.94140625,  90559.2109375 ],
1060                  [  56659.0703125 ,  65907.6796875 ],
1061                  [  35406.3359375 ,  79332.9140625 ],
1062                  [  32755.59375   ,  58126.9765625 ],
1063                  [  51737.94140625,  90559.2109375 ],
1064                  [  56659.0703125 ,  65907.6796875 ],
1065                  [  80468.9453125 ,  96184.0546875 ],
1066                  [  56659.0703125 ,  65907.6796875 ],
1067                  [  52341.70703125,  38563.39453125],
1068                  [  75735.4765625 ,  23762.00585938],
1069                  [  35406.3359375 ,  79332.9140625 ],
1070                  [  56659.0703125 ,  65907.6796875 ],
1071                  [  51737.94140625,  90559.2109375 ],
1072                  [  80468.9453125 ,  96184.0546875 ],
1073                  [ 101485.5390625 ,  45251.9453125 ],
1074                  [ 117768.6640625 ,  85770.4296875 ],
1075                  [  80468.9453125 ,  96184.0546875 ],
1076                  [  56659.0703125 ,  65907.6796875 ],
1077                  [ 101485.5390625 ,  45251.9453125 ],
1078                  [  32187.57617188, 109816.4375    ],
1079                  [  51737.94140625,  90559.2109375 ],
1080                  [  50364.08984375, 110763.1328125 ],
1081                  [ 118715.359375  , 129886.34375   ],
1082                  [  86149.1015625 , 129886.34375   ],
1083                  [  80468.9453125 ,  96184.0546875 ],
1084                  [  80468.9453125 ,  96184.0546875 ],
1085                  [ 117768.6640625 ,  85770.4296875 ],
1086                  [ 118715.359375  , 129886.34375   ],
1087                  [  52341.70703125,  38563.39453125],
1088                  [  56659.0703125 ,  65907.6796875 ],
1089                  [  32755.59375   ,  58126.9765625 ],
1090                  [  51737.94140625,  90559.2109375 ],
1091                  [  31998.23828125,  88799.84375   ],
1092                  [  35406.3359375 ,  79332.9140625 ]]
1093
1094        scaled_points = ensure_numeric(points, num.int)/1000  # Simplify for ease of interpretation
1095
1096        triangles = [[ 0, 1, 2],
1097                     [ 3, 4, 5],
1098                     [ 6, 7, 8],
1099                     [ 9,10,11],
1100                     [12,13,14],
1101                     [15,16,17],
1102                     [18,19,20],
1103                     [21,22,23],
1104                     [24,25,26],
1105                     [27,28,29],
1106                     [30,31,32],
1107                     [33,34,35],
1108                     [36,37,38],
1109                     [39,40,41],
1110                     [42,43,44],
1111                     [45,46,47],
1112                     [48,49,50],
1113                     [51,52,53],
1114                     [54,55,56],
1115                     [57,58,59],
1116                     [60,61,62]]
1117
1118
1119        # First use scaled points for ease of debugging
1120        mesh = Mesh(scaled_points, triangles)
1121        mesh.check_integrity()
1122        P = mesh.get_boundary_polygon()
1123
1124        for p in scaled_points:
1125            assert is_inside_polygon(p, P)           
1126
1127        # Then use original points and test       
1128        mesh = Mesh(points, triangles)
1129        mesh.check_integrity()
1130        P = mesh.get_boundary_polygon()
1131
1132        for p in points:
1133            assert is_inside_polygon(p, P)           
1134
1135        assert num.allclose(P, Pref)   
1136
1137    def test_lone_vertices(self):
1138        a = [2.0, 1.0]
1139        b = [6.0, 2.0]
1140        c = [1.0, 3.0]
1141        d = [2.0, 4.0]
1142
1143        points = [a, b, c, d]
1144        vertices = [[0,1,2]]
1145
1146        mesh = Mesh(points, vertices)
1147        mesh.check_integrity()
1148        loners = mesh.get_lone_vertices()
1149        self.failUnless(loners==[3],
1150                        'FAILED!')
1151
1152       
1153        a = [2.0, 1.0]
1154        b = [6.0, 2.0]
1155        c = [1.0, 3.0]
1156        d = [2.0, 4.0]
1157
1158        points = [d, a, b, c]
1159        vertices = [[3,1,2]]
1160
1161        mesh = Mesh(points, vertices)
1162        mesh.check_integrity()
1163        loners = mesh.get_lone_vertices()
1164        self.failUnless(loners==[0],
1165                        'FAILED!') 
1166
1167    def test_mesh_get_boundary_polygon_with_georeferencing(self):
1168        """test_mesh_get_boundary_polygon_with_georeferencing
1169       
1170        Test that get_boundary_polygon returns absolute coordinates
1171        """
1172       
1173        # test
1174        a = [0.0, 0.0]
1175        b = [4.0, 0.0]
1176        c = [0.0, 4.0]
1177
1178        absolute_points = [a, b, c]
1179        vertices = [[0, 1, 2]]
1180       
1181        geo = Geo_reference(56, 67, -56)
1182
1183        relative_points = geo.change_points_geo_ref(absolute_points)
1184
1185        mesh = Mesh(relative_points, vertices, geo_reference=geo)
1186        boundary_polygon = mesh.get_boundary_polygon()
1187
1188        assert num.allclose(absolute_points, boundary_polygon)
1189
1190    def test_get_triangle_containing_point(self):
1191
1192        a = [0.0, 0.0]
1193        b = [0.0, 2.0]
1194        c = [2.0, 0.0]
1195        d = [0.0, 4.0]
1196        e = [2.0, 2.0]
1197        f = [4.0, 0.0]
1198
1199        points = [a, b, c, d, e, f]
1200        #bac, bce, ecf, dbe
1201        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
1202        mesh = Mesh(points, vertices)
1203       
1204        mesh.check_integrity()
1205
1206
1207        try:
1208            id = mesh.get_triangle_containing_point([3.0, 5.0])
1209        except:
1210            pass
1211        else:
1212            msg = 'Should have caught point outside polygon (Non)'           
1213            raise Exception(msg)
1214           
1215        id = mesh.get_triangle_containing_point([0.5, 1.0])
1216        assert id == 0
1217
1218        id = mesh.get_triangle_containing_point([1.0, 3.0])
1219        assert id == 3       
1220
1221        for i, point in enumerate(mesh.get_centroid_coordinates()):
1222            id = mesh.get_triangle_containing_point(point)
1223            assert id == i       
1224
1225    def test_get_triangle_neighbours(self):
1226        a = [0.0, 0.0]
1227        b = [0.0, 2.0]
1228        c = [2.0,0.0]
1229        e = [2.0, 2.0]
1230        points = [a, b, c, e]
1231        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce
1232        mesh = Mesh(points, vertices)
1233
1234        neighbours = mesh.get_triangle_neighbours(0)
1235        assert num.allclose(neighbours, [-1,1,-1])
1236        neighbours = mesh.get_triangle_neighbours(-10)
1237        assert neighbours == []
1238        neighbours = mesh.get_triangle_neighbours(2)
1239        assert neighbours == []
1240
1241
1242    def test_get_intersecting_segments1(self):
1243        """test_get_intersecting_segments(self):
1244
1245        Very simple test (horizontal lines)
1246        """
1247
1248        # Build test mesh
1249       
1250        # Create basic mesh
1251        # 9 points at (0,0), (0, 1), ..., (2,2)
1252        # 8 triangles enumerated from left bottom to right top.
1253        points, vertices, boundary = rectangular(2, 2, 2, 2)
1254        mesh = Mesh(points, vertices, boundary)
1255
1256        # Very simple horizontal line intersecting
1257        #
1258
1259
1260        for y_line in [0.1, 0.2, 0.314159, 0.41, 0.6, 0.99, 1.01, 1.5, 1.77, 1.9]:
1261            if y_line < 1:
1262                ceiling = 1
1263                floor = 0
1264                intersected_triangles = [0,1,4,5]
1265            elif y_line > 1:
1266                ceiling = 2
1267                floor = 1
1268                intersected_triangles = [2,3,6,7]
1269            else:
1270                raise Exception('this test is not for parallel lines')
1271
1272
1273            line = [[-1,y_line], [3,y_line]]
1274
1275            L = mesh.get_intersecting_segments(line)
1276            assert len(L) == 4
1277
1278           
1279
1280            # Check all normals point straight down etc
1281            total_length = 0
1282            for x in L:
1283                if x.triangle_id % 2 == 0:
1284                    assert num.allclose(x.length, ceiling-y_line)
1285                else:
1286                    assert num.allclose(x.length, y_line-floor)               
1287
1288               
1289                assert num.allclose(x.normal, [0,-1])
1290
1291                assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length)
1292                assert num.allclose(x.segment[0][1], y_line)
1293                assert num.allclose(x.segment[1][1], y_line)               
1294
1295                assert x.triangle_id in intersected_triangles
1296
1297                total_length += x.length
1298
1299            msg = 'Segments do not add up'
1300            assert num.allclose(total_length, 2), msg
1301           
1302
1303    def test_get_intersecting_segments_coinciding(self):
1304        """test_get_intersecting_segments_coinciding(self):
1305
1306        Test that lines coinciding with triangle edges work.
1307        """
1308
1309        # Build test mesh
1310       
1311        # Create basic mesh
1312        # 9 points at (0,0), (0, 1), ..., (2,2)
1313        # 8 triangles enumerated from left bottom to right top.
1314        points, vertices, boundary = rectangular(2, 2, 2, 2)
1315        mesh = Mesh(points, vertices, boundary)
1316        intersected_triangles = [1,5]
1317
1318        # Very simple horizontal line intersecting
1319        #
1320
1321        y_line = 1.0
1322       
1323        line = [[-1,y_line], [3,y_line]]
1324
1325        L = mesh.get_intersecting_segments(line)
1326
1327
1328        msg = 'Only two triangles should be returned'   
1329        assert len(L) == 2, msg   
1330           
1331
1332        # Check all
1333        total_length = 0
1334        for x in L:
1335            assert num.allclose(x.length, 1.0)
1336            assert num.allclose(x.normal, [0,-1])
1337
1338            assert num.allclose(x.segment[1][0], x.segment[0][0] + x.length)
1339            assert num.allclose(x.segment[0][1], y_line)
1340            assert num.allclose(x.segment[1][1], y_line)                           
1341
1342
1343
1344            assert x.triangle_id in intersected_triangles
1345
1346            total_length += x.length
1347
1348        msg = 'Segments do not add up'
1349        assert num.allclose(total_length, 2), msg
1350       
1351
1352    def test_get_intersecting_segments_partially_coinciding(self):
1353        """test_get_intersecting_segments_partially_coinciding(self):
1354
1355        Test that line coinciding with triangle edges work.
1356        But this ones only coincide with parts of the edge.
1357        """
1358
1359        # Build test mesh
1360       
1361        # Create basic mesh
1362        # 9 points at (0,0), (0, 1), ..., (2,2)
1363        # 8 triangles enumerated from left bottom to right top.
1364        points, vertices, boundary = rectangular(2, 2, 2, 2)
1365        mesh = Mesh(points, vertices, boundary)
1366        intersected_triangles = [1,5]
1367
1368        # Horizontal line intersecting along center but stopping
1369        # parway through second triangle's edge
1370        #
1371
1372        y_line = 1.0
1373       
1374        #line = [[0, y_line], [2, y_line]]
1375        line = [[0, y_line], [1.5, y_line]]
1376
1377        L = mesh.get_intersecting_segments(line)
1378        #for x in L:
1379        #    print x
1380
1381        msg = 'Two triangles should be returned'   
1382        assert len(L) == 2, msg   
1383           
1384
1385        # Check all
1386        total_length = 0
1387        for x in L:
1388            if x.triangle_id == 1:
1389                assert num.allclose(x.length, 1)       
1390                assert num.allclose(x.normal, [0, -1])
1391               
1392            if x.triangle_id == 5:
1393                assert num.allclose(x.length, 0.5)
1394                assert num.allclose(x.normal, [0, -1])
1395
1396
1397            assert x.triangle_id in intersected_triangles
1398
1399            total_length += x.length
1400
1401        msg = 'Segments do not add up'
1402        assert num.allclose(total_length, 1.5), msg           
1403
1404
1405
1406    def test_get_intersecting_segments2(self):
1407        """test_get_intersecting_segments(self):
1408
1409        Lines with a slope
1410        """
1411
1412        s2 = sqrt(2.0)/2
1413       
1414
1415        # Build test mesh
1416       
1417        # Create basic mesh
1418        # 9 points at (0,0), (0, 1), ..., (2,2)
1419        # 8 triangles enumerated from left bottom to right top.
1420        points, vertices, boundary = rectangular(2, 2, 2, 2)
1421        mesh = Mesh(points, vertices, boundary)
1422       
1423
1424        # Diagonal cutting through a vertex and hypothenuses
1425        line = [[0, 2], [2, 0]]
1426        intersected_triangles = [3,2,5,4]               
1427
1428        L = mesh.get_intersecting_segments(line)
1429        assert len(L) == 4
1430
1431        #print L
1432       
1433        # Check all segments
1434        total_length = 0
1435        for i, x in enumerate(L): 
1436            assert num.allclose(x.length, s2)
1437            assert num.allclose(x.normal, [-s2, -s2])
1438            assert num.allclose(sum(x.normal**2), 1)
1439           
1440            assert x.triangle_id in intersected_triangles
1441
1442            total_length += x.length
1443
1444        msg = 'Segments do not add up'
1445        assert num.allclose(total_length, 4*s2), msg
1446
1447
1448        # Diagonal cutting through a vertex and hypothenuses (reversed)
1449        line = [[2, 0], [0, 2]]
1450        intersected_triangles = [3,2,5,4]               
1451
1452        L = mesh.get_intersecting_segments(line)
1453        assert len(L) == 4
1454
1455        #print L
1456       
1457        # Check all segments
1458        total_length = 0
1459        for i, x in enumerate(L): 
1460            assert num.allclose(x.length, s2)
1461            assert num.allclose(x.normal, [s2, s2])
1462            assert num.allclose(sum(x.normal**2), 1)
1463           
1464            assert x.triangle_id in intersected_triangles
1465
1466            total_length += x.length
1467
1468        msg = 'Segments do not add up'
1469        assert num.allclose(total_length, 4*s2), msg                   
1470
1471
1472
1473        # Diagonal coinciding with hypothenuses
1474        line = [[2, 2], [0, 0]]
1475        intersected_triangles = [6,0]               
1476
1477        L = mesh.get_intersecting_segments(line)
1478        assert len(L) == 2
1479
1480        #print L
1481       
1482        # Check all segments
1483        total_length = 0
1484        for i, x in enumerate(L): 
1485            assert num.allclose(x.length, 2*s2)
1486            assert num.allclose(x.normal, [-s2, s2])
1487            assert num.allclose(sum(x.normal**2), 1)
1488           
1489            assert x.triangle_id in intersected_triangles
1490
1491            total_length += x.length
1492
1493        msg = 'Segments do not add up'
1494        assert num.allclose(total_length, 4*s2), msg                       
1495
1496
1497        # Diagonal coinciding with hypothenuses (reversed)
1498        line = [[0, 0], [2, 2]]
1499        intersected_triangles = [6,0]               
1500
1501        L = mesh.get_intersecting_segments(line)
1502        assert len(L) == 2
1503
1504        #print L
1505       
1506        # Check all segments
1507        total_length = 0
1508        for i, x in enumerate(L): 
1509            assert num.allclose(x.length, 2*s2)
1510            assert num.allclose(x.normal, [s2, -s2])
1511            assert num.allclose(sum(x.normal**2), 1)
1512           
1513            assert x.triangle_id in intersected_triangles
1514
1515            total_length += x.length
1516
1517        msg = 'Segments do not add up'
1518        assert num.allclose(total_length, 4*s2), msg                       
1519
1520
1521
1522        # line with slope [1, -1] cutting through vertices of tri 7 and 6
1523        line = [[1, 2], [2, 1]]
1524        intersected_triangles = [7,6]               
1525
1526        L = mesh.get_intersecting_segments(line)
1527        assert len(L) == 2
1528
1529        #print L
1530       
1531        # Check all segments
1532        total_length = 0
1533        for i, x in enumerate(L): 
1534            assert num.allclose(x.length, s2)
1535            assert num.allclose(x.normal, [-s2, -s2])
1536            assert num.allclose(sum(x.normal**2), 1)
1537           
1538            assert x.triangle_id in intersected_triangles
1539
1540            total_length += x.length
1541
1542        msg = 'Segments do not add up'
1543        assert num.allclose(total_length, 2*s2), msg
1544
1545
1546        # Arbitrary line with slope [1, -1] cutting through tri 7 and 6
1547        line = [[1.1, 2], [2.1, 1]]
1548        intersected_triangles = [7,6]               
1549
1550        L = mesh.get_intersecting_segments(line)
1551        assert len(L) == 2
1552       
1553        # Check all segments
1554        total_length = 0
1555        for i, x in enumerate(L): 
1556            assert num.allclose(x.normal, [-s2, -s2])
1557            assert num.allclose(sum(x.normal**2), 1)
1558
1559            msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles)
1560            assert x.triangle_id in intersected_triangles, msg
1561           
1562
1563
1564    def test_get_intersecting_segments3(self):
1565        """test_get_intersecting_segments(self):
1566
1567        Check that line can stop inside a triangle
1568       
1569        """
1570
1571
1572
1573        s2 = sqrt(2.0)/2
1574       
1575
1576        # Build test mesh
1577       
1578        # Create basic mesh
1579        # 9 points at (0,0), (0, 1), ..., (2,2)
1580        # 8 triangles enumerated from left bottom to right top.
1581        points, vertices, boundary = rectangular(2, 2, 2, 2)
1582        mesh = Mesh(points, vertices, boundary)
1583       
1584
1585        # Line cutting through one triangle and ending on its edge
1586        line = [[0.5, 3], [0.5, 1.5]]
1587        intersected_triangles = [3]               
1588
1589        L = mesh.get_intersecting_segments(line)
1590        assert len(L) == 1
1591        assert L[0].triangle_id == 3
1592        assert num.allclose(L[0].length, 0.5)       
1593        assert num.allclose(L[0].normal, [-1,0])               
1594
1595
1596
1597        # Now try to shorten it so that its endpoint falls short of the far edge
1598        line = [[0.5, 3], [0.5, 1.6]]
1599        intersected_triangles = [3]               
1600
1601        L = mesh.get_intersecting_segments(line)
1602        assert len(L) == 1
1603        assert L[0].triangle_id == 3
1604        assert num.allclose(L[0].length, 0.4)
1605        assert num.allclose(L[0].normal, [-1,0])
1606
1607        intersected_triangles = [3]
1608
1609        # Now the same, but with direction changed
1610        line = [[0.5, 3], [0.5, 1.6]]
1611        line = [[0.5, 1.6], [0.5, 3]]       
1612        intersected_triangles = [3]               
1613
1614        L = mesh.get_intersecting_segments(line)
1615        assert len(L) == 1
1616        assert L[0].triangle_id == 3
1617        assert num.allclose(L[0].length, 0.4)
1618        assert num.allclose(L[0].normal, [1,0])               
1619       
1620
1621           
1622
1623    def test_get_intersecting_segments4(self):
1624        """test_get_intersecting_segments(self):
1625
1626        Simple poly line
1627       
1628        """
1629
1630
1631
1632        s2 = sqrt(2.0)/2
1633       
1634
1635        # Build test mesh
1636       
1637        # Create basic mesh
1638        # 9 points at (0,0), (0, 1), ..., (2,2)
1639        # 8 triangles enumerated from left bottom to right top.
1640        points, vertices, boundary = rectangular(2, 2, 2, 2)
1641        mesh = Mesh(points, vertices, boundary)
1642       
1643
1644        # Polyline with three segments cutting through mesh
1645        line = [[0.5, 3], [0.5, 1.5], [1,1]]
1646
1647        L = mesh.get_intersecting_segments(line)
1648        assert len(L) == 2
1649
1650        for x in L:
1651            if x.triangle_id == 3:
1652                assert num.allclose(x.length, 0.5)       
1653                assert num.allclose(x.normal, [-1,0])
1654               
1655            if x.triangle_id == 2:
1656                assert num.allclose(x.length, s2)
1657                assert num.allclose(x.normal, [-s2,-s2])
1658
1659
1660
1661    def test_get_intersecting_segments5(self):
1662        """test_get_intersecting_segments(self):
1663
1664        More complex poly line
1665       
1666        """
1667
1668
1669
1670        s2 = sqrt(2.0)/2
1671       
1672
1673        # Build test mesh
1674       
1675        # Create basic mesh
1676        # 9 points at (0,0), (0, 1), ..., (2,2)
1677        # 8 triangles enumerated from left bottom to right top.
1678        points, vertices, boundary = rectangular(2, 2, 2, 2)
1679        mesh = Mesh(points, vertices, boundary)
1680       
1681
1682        # Polyline with three segments cutting through mesh
1683        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75]] 
1684
1685        L = mesh.get_intersecting_segments(line)
1686        assert len(L) == 3
1687
1688        for x in L:
1689            if x.triangle_id == 3:
1690                assert num.allclose(x.length, 0.5)       
1691                assert num.allclose(x.normal, [-1,0])
1692               
1693            if x.triangle_id == 2:
1694                msg = str(x.length)
1695                assert num.allclose(x.length, s2), msg
1696                assert num.allclose(x.normal, [-s2,-s2])
1697
1698            if x.triangle_id == 5:
1699                segvec = num.array([line[2][0]-1,
1700                                    line[2][1]-1])
1701                msg = str(x.length)
1702                assert num.allclose(x.length, sqrt(sum(segvec**2))), msg
1703                assert num.allclose(x.normal, [-s2,-s2])                                               
1704
1705
1706    def test_get_intersecting_segments6(self):
1707        """test_get_intersecting_segments(self):
1708
1709        Even more complex poly line, where line breaks within triangle 5
1710
1711        5 segments are returned even though only four triangles [3,2,5,6] are touched.
1712        Triangle 5 therefor has two segments in it.
1713       
1714        """
1715
1716
1717
1718        s2 = sqrt(2.0)/2
1719       
1720
1721        # Build test mesh
1722       
1723        # Create basic mesh
1724        # 9 points at (0,0), (0, 1), ..., (2,2)
1725        # 8 triangles enumerated from left bottom to right top.
1726        points, vertices, boundary = rectangular(2, 2, 2, 2)
1727        mesh = Mesh(points, vertices, boundary)
1728       
1729
1730        # Polyline with three segments cutting through mesh
1731        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75], [2.25, 1.75]]
1732
1733        L = mesh.get_intersecting_segments(line)
1734        #for x in L:
1735        #    print x
1736
1737        assert len(L) == 5
1738
1739        for x in L:
1740            if x.triangle_id == 3:
1741                assert num.allclose(x.length, 0.5)       
1742                assert num.allclose(x.normal, [-1,0])
1743               
1744            if x.triangle_id == 2:
1745                msg = str(x.length)
1746                assert num.allclose(x.length, s2), msg
1747                assert num.allclose(x.normal, [-s2,-s2])
1748
1749            if x.triangle_id == 5:
1750                if x.segment == ((1.0, 1.0), (1.25, 0.75)):                   
1751                    segvec = num.array([line[2][0]-1,
1752                                        line[2][1]-1])
1753                    msg = str(x.length)
1754                    assert num.allclose(x.length, sqrt(sum(segvec**2))), msg
1755                    assert num.allclose(x.normal, [-s2,-s2])
1756                elif x.segment == ((1.25, 0.75), (1.5, 1.0)):
1757                    segvec = num.array([1.5-line[2][0],
1758                                        1.0-line[2][1]])
1759                   
1760                    assert num.allclose(x.length, sqrt(sum(segvec**2))), msg
1761                    assert num.allclose(x.normal, [s2,-s2])
1762                else:
1763                    msg = 'Unknown segment: %s' %x.segment
1764                    raise Exception(msg)
1765               
1766
1767                   
1768            if x.triangle_id == 6:
1769                assert num.allclose(x.normal, [s2,-s2])
1770                assert num.allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
1771
1772
1773      # Internal test that sum of line segments add up
1774      # to length of input line
1775      #
1776      # Could be useful perhaps
1777      #
1778      #xi1 = line[1][0]
1779      #eta1 = line[1][1]
1780      #linevector = num.array([xi1-xi0, eta1-eta0])
1781      #linelength = sqrt(sum(linevector**2))
1782      #
1783      #segmentlength = 0
1784      #for segment in triangle_intersections:
1785      #    vector = array([segment[1][0] - segment[0][0],
1786      #                    segment[1][1] - segment[0][1]])
1787      #    length = sqrt(sum(vector**2))     
1788      #    segmentlength += length
1789      #
1790      #msg = 'Sum of intersecting segments do not add up'   
1791      #assert allclose(segmentlength, linelength), msg   
1792
1793
1794
1795
1796    def test_get_intersecting_segments7(self):
1797        """test_get_intersecting_segments(self):
1798
1799        Check that line can stop inside a triangle - this is from
1800        flow throug a cross sections example in test_datamanager.
1801       
1802        """
1803
1804        # Build test mesh
1805        width = 5
1806        length = 100
1807        t_end = 1
1808        points, vertices, boundary = rectangular(length, width,
1809                                                 length, width)
1810
1811        mesh = Mesh(points, vertices, boundary)
1812       
1813
1814        # A range of partial lines
1815        x = length/2.
1816        for i in range(10):
1817            start_point = [length/2., i*width/10.]
1818            #print
1819            #print start_point
1820                           
1821            line = [start_point, [length/2., width]]
1822 
1823            L = mesh.get_intersecting_segments(line)
1824
1825            if start_point[1] < 1:
1826                assert len(L) == 5
1827               
1828           
1829            total_length = 0   
1830            for x in L:
1831                total_length += x.length
1832               
1833
1834            ref_length = line[1][1] - line[0][1]
1835            #print ref_length, total_length
1836            assert num.allclose(total_length, ref_length)
1837
1838
1839#-------------------------------------------------------------
1840
1841if __name__ == "__main__":
1842    suite = unittest.makeSuite(Test_Mesh, 'test')
1843    runner = unittest.TextTestRunner()
1844    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.