source: anuga_core/source/anuga/pyvolution/test_neighbour_mesh.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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 anuga.config import epsilon
15from Numeric import allclose, array
16
17from anuga.coordinate_transforms.geo_reference import Geo_reference
18from anuga.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 anuga.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 anuga.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.