source: inundation/pyvolution/test_mesh.py @ 3072

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

Made boundary polygon return absolute UTM coordinates as per ticket:81
File size: 26.7 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 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 mesh import Mesh
733        from Numeric import zeros, Float
734       
735
736        #Points
737        a = [0.0, 0.0] #0
738        b = [0.0, 0.5] #1
739        c = [0.0, 1.0] #2
740        d = [0.5, 0.0] #3
741        e = [0.5, 0.5] #4
742        f = [1.0, 0.0] #5
743        g = [1.0, 0.5] #6
744        h = [1.0, 1.0] #7
745        i = [1.5, 0.5] #8
746
747        points = [a, b, c, d, e, f, g, h, i]
748
749        #dea, bae, bec, fgd,
750        #edg, ghe, gfi, gih
751        vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
752                     [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
753
754        mesh = Mesh(points, vertices)
755
756        mesh.check_integrity()
757
758        P = mesh.get_boundary_polygon()
759
760        assert len(P) == 8
761        assert allclose(P, [a, d, f, i, h, e, c, b])
762
763        for p in points:
764            #print p, P
765            assert is_inside_polygon(p, P)
766
767
768    def test_boundary_polygon_III(self):
769        """Same as II but vertices ordered differently
770        """
771
772        from mesh import Mesh
773        from Numeric import zeros, Float
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 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_IV(self):
809        """Reproduce test test_spatio_temporal_file_function_time
810        from test_util.py that looked as if it produced the wrong boundary
811        """
812
813        from mesh import Mesh
814        from Numeric import zeros, Float
815        from mesh_factory import rectangular       
816
817        #Create a domain to hold test grid
818        #(0:15, -20:10)
819        points, vertices, boundary =\
820                rectangular(4, 4, 15, 30, origin = (0, -20))       
821
822        #####
823        mesh = Mesh(points, vertices)
824        mesh.check_integrity()
825
826        P = mesh.get_boundary_polygon()
827
828        #print P
829        assert len(P) == 16
830        for p in points:
831            assert is_inside_polygon(p, P)
832
833
834
835        #####
836        mesh = Mesh(points, vertices, boundary)
837        mesh.check_integrity()
838
839        P = mesh.get_boundary_polygon()
840
841       
842        #print P, len(P)
843        assert len(P) == 16
844
845        for p in points:
846            assert is_inside_polygon(p, P)
847
848        #print mesh.statistics()   
849
850
851
852    def test_boundary_polygon_V(self):
853        """Create a discontinuous mesh (duplicate vertices)
854        and check that boundary is as expected
855       
856        """
857        from mesh import Mesh
858        from Numeric import zeros, Float
859       
860
861        #Points
862        a = [0.0, 0.0] #0
863        b = [0.0, 0.5] #1
864        c = [0.0, 1.0] #2
865        d = [0.5, 0.0] #3
866        e = [0.5, 0.5] #4
867        f = [1.0, 0.0] #5
868        g = [1.0, 0.5] #6
869        h = [1.0, 1.0] #7
870        i = [1.5, 0.5] #8
871
872        #Duplicate points for triangles edg [4,3,6] (central) and
873        #gid [6,8,7] (top right boundary) to them disconnected
874        #from the others
875
876        e0 = [0.5, 0.5] #9
877        d0 = [0.5, 0.0] #10
878        g0 = [1.0, 0.5] #11
879        i0 = [1.5, 0.5] #12       
880       
881
882        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
883
884
885
886        #dea, bae, bec, fgd,
887        #edg, ghe, gfi, gih
888        #vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
889        #             [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
890
891
892        #dea, bae, bec, fgd,
893        #e0d0g0, ghe, gfi, g0i0h
894        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
895                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]       
896
897        mesh = Mesh(points, vertices)
898
899        mesh.check_integrity()
900
901        P = mesh.get_boundary_polygon()
902
903        #print P
904       
905        assert len(P) == 8
906        assert allclose(P, [a, d, f, i, h, e, c, b])
907        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)])
908       
909
910        for p in points:
911            #print p, P
912            assert is_inside_polygon(p, P)
913
914
915
916    def test_lone_vertices(self):
917        a = [2.0, 1.0]
918        b = [6.0, 2.0]
919        c = [1.0, 3.0]
920        d = [2.0, 4.0]
921
922        points = [a, b, c, d]
923        vertices = [[0,1,2]]
924
925        mesh = Mesh(points, vertices)
926        mesh.check_integrity()
927        loners = mesh.get_lone_vertices()
928        self.failUnless(loners==[3],
929                        'FAILED!')
930
931       
932        a = [2.0, 1.0]
933        b = [6.0, 2.0]
934        c = [1.0, 3.0]
935        d = [2.0, 4.0]
936
937        points = [d, a, b, c]
938        vertices = [[3,1,2]]
939
940        mesh = Mesh(points, vertices)
941        mesh.check_integrity()
942        loners = mesh.get_lone_vertices()
943        self.failUnless(loners==[0],
944                        'FAILED!') 
945
946    def test_mesh_get_boundary_polygon_with_georeferencing(self):
947     
948        # test
949        a = [0.0, 0.0]
950        b = [4.0, 0.0]
951        c = [0.0, 4.0]
952
953        absolute_points = [a, b, c]
954        vertices = [[0,1,2]]
955       
956        geo = Geo_reference(56,67,-56)
957
958        relative_points = geo.change_points_geo_ref(absolute_points)
959
960        #print 'Relative', relative_points
961        #print 'Absolute', absolute_points       
962
963        mesh = Mesh(relative_points, vertices, geo_reference=geo)
964        boundary_polygon = mesh.get_boundary_polygon()
965
966        assert allclose(absolute_points, boundary_polygon)
967       
968#-------------------------------------------------------------
969if __name__ == "__main__":
970    #suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
971    suite = unittest.makeSuite(Test_Mesh,'test')
972    runner = unittest.TextTestRunner()
973    runner.run(suite)
974
975
976
977
Note: See TracBrowser for help on using the repository browser.