source: branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 15 years ago

numpy changes.

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