source: branches/inundation-numpy-branch/pyvolution/test_mesh.py @ 7248

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