source: inundation/pyvolution/test_mesh.py @ 1812

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

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