source: inundation/pyvolution/test_neighbour_mesh.py @ 3354

Last change on this file since 3354 was 3086, checked in by ole, 19 years ago

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