source: branches/numpy/anuga/pmesh/test_mesh_interface.py @ 7184

Last change on this file since 7184 was 7035, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk.
File size: 34.1 KB
Line 
1#!/usr/bin/env python
2
3import tempfile
4import unittest
5import os
6from anuga.pmesh.mesh import importMeshFromFile
7from anuga.pmesh.mesh_interface import create_mesh_from_regions
8
9from anuga.pmesh.mesh_interface import _create_mesh_from_regions
10
11from load_mesh.loadASCII import *
12from anuga.utilities.polygon import is_inside_polygon
13from anuga.coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
14
15
16class TestCase(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23    def test_create_mesh_from_regions(self):
24        x=-500
25        y=-1000
26        mesh_geo = geo_reference=Geo_reference(56, x, y)
27
28        # These are the absolute values
29        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
30
31        x_p = -10
32        y_p = -40
33        geo_ref_poly = Geo_reference(56, x_p, y_p)
34        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
35
36        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
37
38        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
39        inner1_polygon = geo_ref_poly.\
40                            change_points_geo_ref(inner1_polygon_absolute)
41
42        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
43        inner2_polygon = geo_ref_poly.\
44                            change_points_geo_ref(inner2_polygon_absolute)
45
46        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
47
48        m = create_mesh_from_regions(polygon,
49                                     boundary_tags,
50                                     10000000,
51                                     interior_regions=interior_regions,
52                                     poly_geo_reference=geo_ref_poly,
53                                     mesh_geo_reference=mesh_geo)
54
55        # Test the mesh instance
56        self.failUnless(len(m.regions)==3, 'FAILED!')
57        segs = m.getUserSegments()
58        self.failUnless(len(segs)==12, 'FAILED!')
59        self.failUnless(len(m.userVertices)==12, 'FAILED!')
60        self.failUnless(segs[0].tag=='walls', 'FAILED!')
61        self.failUnless(segs[1].tag=='walls', 'FAILED!')
62        self.failUnless(segs[2].tag=='bom', 'FAILED!')
63        self.failUnless(segs[3].tag=='bom', 'FAILED!')
64
65        # Assuming the order of the region points is known.
66        # (This isn't true, if you consider create_mesh_from_regions
67        # a black box)
68        poly_point = m.getRegions()[0]
69
70        # poly_point values are relative to the mesh geo-ref
71        # make them absolute
72        msg = ('Expected point (%s,%s) to be inside polygon %s'
73               % (str(poly_point.x+x), str(poly_point.y+y),
74                  str(polygon_absolute)))
75        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
76                                          polygon_absolute, closed=False),
77                        msg)
78#                        'FAILED!')
79
80        # Assuming the order of the region points is known.
81        # (This isn't true, if you consider create_mesh_from_regions
82        # a black box)
83        poly_point = m.getRegions()[1]
84
85        # poly_point values are relative to the mesh geo-ref
86        # make them absolute
87        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
88                                          inner1_polygon_absolute,
89                                          closed=False),
90                        'FAILED!')
91
92        # Assuming the order of the region points is known.
93        # (This isn't true, if you consider create_mesh_from_regions
94        # a black box)
95        poly_point = m.getRegions()[2]
96
97        # poly_point values are relative to the mesh geo-ref
98        # make them absolute
99        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
100                                          inner2_polygon_absolute,
101                                          closed=False),
102                        'FAILED!')
103
104    def test_create_mesh_from_regions_with_caching(self):
105        x=-500
106        y=-1000
107        mesh_geo = geo_reference=Geo_reference(56, x, y)
108
109        # These are the absolute values
110        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
111
112        x_p = -10
113        y_p = -40
114        geo_ref_poly = Geo_reference(56, x_p, y_p)
115        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
116
117        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
118
119        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
120        inner1_polygon = geo_ref_poly.\
121                            change_points_geo_ref(inner1_polygon_absolute)
122
123        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
124        inner2_polygon = geo_ref_poly.\
125                             change_points_geo_ref(inner2_polygon_absolute)
126
127        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
128
129        interior_holes = None
130
131        # Clear cache first
132        from anuga.caching import cache
133
134        cache(_create_mesh_from_regions,
135              (polygon, boundary_tags),
136              {'minimum_triangle_angle': 28.0,
137               'maximum_triangle_area': 10000000,
138               'interior_regions': interior_regions,
139               'interior_holes': interior_holes,
140               'poly_geo_reference': geo_ref_poly,
141               'mesh_geo_reference': mesh_geo,
142               'verbose': False},
143              verbose=False,
144              clear=1)
145
146        m = create_mesh_from_regions(polygon,
147                                     boundary_tags,
148                                     maximum_triangle_area=10000000,
149                                     interior_regions=interior_regions,
150                                     poly_geo_reference=geo_ref_poly,
151                                     mesh_geo_reference=mesh_geo,
152                                     verbose=False,
153                                     use_cache=True)
154
155
156        # Test the mesh instance
157        self.failUnless(len(m.regions)==3, 'FAILED!')
158        segs = m.getUserSegments()
159        self.failUnless(len(segs)==12, 'FAILED!')
160        self.failUnless(len(m.userVertices)==12, 'FAILED!')
161        self.failUnless(segs[0].tag=='walls', 'FAILED!')
162        self.failUnless(segs[1].tag=='walls', 'FAILED!')
163        self.failUnless(segs[2].tag=='bom', 'FAILED!')
164        self.failUnless(segs[3].tag=='bom', 'FAILED!')
165
166        # Assuming the order of the region points is known.
167        # (This isn't true, if you consider create_mesh_from_regions
168        # a black box)
169        poly_point = m.getRegions()[0]
170
171        # poly_point values are relative to the mesh geo-ref
172        # make them absolute
173        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
174                                          polygon_absolute,
175                                          closed=False),
176                        'FAILED!')
177
178        # Assuming the order of the region points is known.
179        # (This isn't true, if you consider create_mesh_from_regions
180        # a black box)
181        poly_point = m.getRegions()[1]
182
183        # poly_point values are relative to the mesh geo-ref
184        # make them absolute
185        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
186                                          inner1_polygon_absolute,
187                                          closed=False),
188                        'FAILED!')
189
190        # Assuming the order of the region points is known.
191        # (This isn't true, if you consider create_mesh_from_regions
192        # a black box)
193        poly_point = m.getRegions()[2]
194
195        # poly_point values are relative to the mesh geo-ref
196        # make them absolute
197        self.failUnless(is_inside_polygon([poly_point.x+x, poly_point.y+y],
198                                          inner2_polygon_absolute,
199                                          closed=False),
200                        'FAILED!')
201
202        # Now create m using cached values
203        m_cache = create_mesh_from_regions(polygon,
204                                           boundary_tags,
205                                           10000000,
206                                           interior_regions=interior_regions,
207                                           poly_geo_reference=geo_ref_poly,
208                                           mesh_geo_reference=mesh_geo,
209                                           verbose=False,
210                                           use_cache=True)
211
212    def test_create_mesh_from_regions2(self):
213        # These are the absolute values
214        min_x = -10
215        min_y = -88
216        polygon_absolute = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
217
218        x_p = -10
219        y_p = -40
220        zone = 808
221        geo_ref_poly = Geo_reference(zone, x_p, y_p)
222        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
223
224        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
225
226        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
227        inner1_polygon = geo_ref_poly.\
228                            change_points_geo_ref(inner1_polygon_absolute)
229
230        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
231        inner2_polygon = geo_ref_poly.\
232                            change_points_geo_ref(inner2_polygon_absolute)
233
234        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
235        m = create_mesh_from_regions(polygon,
236                                     boundary_tags,
237                                     10000000,
238                                     interior_regions=interior_regions,
239                                     poly_geo_reference=geo_ref_poly)
240
241        # Test the mesh instance
242        self.failUnless(len(m.regions)==3, 'FAILED!')
243        segs = m.getUserSegments()
244        self.failUnless(len(segs)==12, 'FAILED!')
245        self.failUnless(len(m.userVertices)==12, 'FAILED!')
246        self.failUnless(segs[0].tag=='walls', 'FAILED!')
247        self.failUnless(segs[1].tag=='walls', 'FAILED!')
248        self.failUnless(segs[2].tag=='bom', 'FAILED!')
249        self.failUnless(segs[3].tag=='bom', 'FAILED!')
250        self.failUnless(m.geo_reference.get_zone()==zone, 'FAILED!')
251        self.failUnless(m.geo_reference.get_xllcorner()==min_x, 'FAILED!')
252        self.failUnless(m.geo_reference.get_yllcorner()==min_y, 'FAILED!')
253
254    def test_create_mesh_from_regions3(self):
255        # These are the absolute values
256        min_x = -10
257        min_y = -88
258        polygon = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
259
260
261        x_p = -10
262        y_p = -40
263        geo_ref_poly = Geo_reference(56, x_p, y_p)
264
265        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
266
267        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
268        inner1_polygon = geo_ref_poly.\
269                            change_points_geo_ref(inner1_polygon_absolute)
270
271        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
272        inner2_polygon = geo_ref_poly.\
273                            change_points_geo_ref(inner2_polygon_absolute)
274
275        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
276        m = create_mesh_from_regions(polygon,
277                                     boundary_tags,
278                                     10000000,
279                                     interior_regions=interior_regions)
280
281        # Test the mesh instance
282        self.failUnless(len(m.regions) == 3, 'FAILED!')
283        segs = m.getUserSegments()
284        self.failUnless(len(segs) == 12, 'FAILED!')
285        self.failUnless(len(m.userVertices) == 12, 'FAILED!')
286        self.failUnless(segs[0].tag == 'walls', 'FAILED!')
287        self.failUnless(segs[1].tag == 'walls', 'FAILED!')
288        self.failUnless(segs[2].tag == 'bom', 'FAILED!')
289        self.failUnless(segs[3].tag == 'bom', 'FAILED!')
290        self.failUnless(m.geo_reference.get_zone() == DEFAULT_ZONE, 'FAILED!')
291        self.failUnless(m.geo_reference.get_xllcorner() == min_x, 'FAILED!')
292        self.failUnless(m.geo_reference.get_yllcorner() == min_y, 'FAILED!')
293
294    def test_create_mesh_from_regions4(self):
295        file_name = tempfile.mktemp('.tsh')
296
297        # These are the absolute values
298        density_outer = 1000
299        min_outer = 0
300        max_outer = 1000
301        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
302                         [max_outer,max_outer], [min_outer,max_outer]]
303
304        density_inner1 = 10000000
305        inner_buffer = 100
306        min_inner1 = min_outer + inner_buffer
307        max_inner1 = max_outer - inner_buffer
308        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
309                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
310
311        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
312
313        interior_regions = [(inner1_polygon, density_inner1)]
314        create_mesh_from_regions(polygon_outer,
315                                 boundary_tags,
316                                 density_outer,
317                                 interior_regions=interior_regions,
318                                 filename=file_name, verbose=False)
319
320        m = importMeshFromFile(file_name)
321
322        self.failUnless(len(m.getTriangulation()) <= 900,
323                        'Test mesh interface failed!')
324        self.failUnless(len(m.getTriangulation()) >= 200,
325                        'Test mesh interface failed!')
326
327        create_mesh_from_regions(polygon_outer,
328                                 boundary_tags,
329                                 interior_regions=interior_regions,
330                                 filename=file_name,
331                                 verbose=False)
332
333        m = importMeshFromFile(file_name)
334
335        self.failUnless(len(m.getTriangulation()) <= 100,
336                        'Test mesh interface failed!')
337
338        os.remove(file_name)
339
340    def test_create_mesh_from_regions5(self):
341        file_name = tempfile.mktemp('.tsh')
342
343        # These are the absolute values
344        density_outer = 10000000
345        min_outer = 0
346        max_outer = 1000
347        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
348                         [max_outer,max_outer], [min_outer,max_outer]]
349
350        density_inner1 = 1000
351        inner_buffer = 100
352        min_inner1 = min_outer + inner_buffer
353        max_inner1 = max_outer - inner_buffer
354        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
355                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
356
357        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
358
359        interior_regions = [(inner1_polygon, density_inner1)]
360        create_mesh_from_regions(polygon_outer,
361                                 boundary_tags,
362                                 density_outer,
363                                 interior_regions=interior_regions,
364                                 filename=file_name,
365                                 verbose=False)
366
367        m = importMeshFromFile(file_name)
368        self.failUnless(len(m.getTriangulation()) <= 2000,
369                        'Test mesh interface failed!')
370        self.failUnless(len(m.getTriangulation()) >= 900,
371                        'Test mesh interface failed!')
372
373        os.remove(file_name)
374
375    def test_create_mesh_from_regions6(self):
376        file_name = tempfile.mktemp('.tsh')
377
378        # These are the absolute values
379        density_outer = 1000
380        min_outer = 0
381        max_outer = 1000
382        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
383                         [max_outer,max_outer], [min_outer,max_outer]]
384
385        delta = 10
386        density_inner1 = 1000
387        min_inner1 = min_outer + delta
388        max_inner1 = max_outer - delta
389        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
390                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
391
392        density_inner2 = 10000000
393        min_inner2 = min_outer +  2*delta
394        max_inner2 = max_outer -  2*delta
395        inner2_polygon = [[min_inner2,min_inner2], [max_inner2,min_inner2],
396                          [max_inner2,max_inner2], [min_inner2,max_inner2]]
397
398        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
399
400        interior_regions = [(inner1_polygon, density_inner1),
401                            (inner2_polygon, density_inner2)]
402        create_mesh_from_regions(polygon_outer,
403                                 boundary_tags,
404                                 density_outer,
405                                 interior_regions=interior_regions,
406                                 filename=file_name,
407                                 verbose=False)
408
409        m = importMeshFromFile(file_name)
410        self.failUnless(len(m.getTriangulation()) <= 2000,
411                        'Test mesh interface failed!')
412        self.failUnless(len(m.getTriangulation()) >= 900,
413                        'Test mesh interface failed!')
414
415        os.remove(file_name)
416
417    def test_create_mesh_from_regions7(self):
418        file_name = tempfile.mktemp('.tsh')
419
420        # These are the absolute values
421        density_outer = 1001
422        min_outer = 0
423        max_outer = 1000
424        polygon_outer = [[min_outer,min_outer], [max_outer,min_outer],
425                         [max_outer,max_outer], [min_outer,max_outer]]
426
427        delta = 10
428        density_inner1 = 100000000
429        min_inner1 = min_outer + delta
430        max_inner1 = max_outer - delta
431        inner1_polygon = [[min_inner1,min_inner1], [max_inner1,min_inner1],
432                          [max_inner1,max_inner1], [min_inner1,max_inner1]]
433
434        density_inner2 = 1000
435        min_inner2 = min_outer +  2*delta
436        max_inner2 = max_outer -  2*delta
437        inner2_polygon = [[min_inner2,min_inner2], [max_inner2,min_inner2],
438                          [max_inner2,max_inner2], [min_inner2,max_inner2]]
439
440        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
441
442        # Note the list order is important
443        # The last region added will be the region triangle uses,
444        # if two regions points are in the same bounded area.
445        interior_regions = [(inner2_polygon, density_inner2),
446                            (inner1_polygon, density_inner1)]
447        create_mesh_from_regions(polygon_outer,
448                                 boundary_tags,
449                                 density_outer,
450                                 interior_regions=interior_regions,
451                                 filename=file_name,
452                                 verbose=False)
453
454        m = importMeshFromFile(file_name)
455        self.failUnless(len(m.getTriangulation()) <= 3000,
456                        'Test mesh interface failed!')
457        self.failUnless(len(m.getTriangulation()) >= 2000,
458                        'Test mesh interface failed!')
459
460        os.remove(file_name)
461
462    def test_create_mesh_from_regions_interior_regions(self):
463        '''Test that create_mesh_from_regions fails when an interior
464        region is outside bounding polygon.
465        '''
466
467        # These are the absolute values
468        min_x = 10
469        min_y = 88
470        polygon = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
471
472        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
473
474        # This one is inside bounding polygon - should pass
475        inner_polygon = [[800,400], [900,500], [800,600]]
476
477        interior_regions = [(inner_polygon, 5)]
478        m = create_mesh_from_regions(polygon,
479                                     boundary_tags,
480                                     10000000,
481                                     interior_regions=interior_regions)
482
483        # This one sticks outside bounding polygon - should fail
484        inner_polygon = [[800,400], [900,500], [800,600], [200, 995]]
485        inner_polygon1 = [[800,400], [1100,500], [800,600]]
486        interior_regions = [[inner_polygon, 50], [inner_polygon1, 50]]
487
488        try:
489            m = create_mesh_from_regions(polygon,
490                                         boundary_tags,
491                                         10000000,
492                                         interior_regions=interior_regions,
493                                         verbose=False)
494        except:
495            pass
496        else:
497            msg = 'Interior polygon sticking outside bounding polygon should '
498            msg += 'cause an Exception to be raised'
499            raise Exception, msg
500
501    def test_create_mesh_from_regions_interior_regions1(self):
502        '''Test that create_mesh_from_regions fails
503        when an interior region is outside bounding polygon.
504        '''
505
506        # These are the values
507        d0 = [310000, 7690000]
508        d1 = [280000, 7690000]
509        d2 = [270000, 7645000]
510        d3 = [240000, 7625000]
511        d4 = [270000, 7580000]
512        d5 = [300000, 7590000]
513        d6 = [340000, 7610000]
514        poly_all = [d0, d1, d2, d3, d4, d5, d6]
515
516        i0 = [304000, 7607000]
517        i1 = [302000, 7605000]
518        i2 = [304000, 7603000]
519        i3 = [307000, 7602000]
520        i4 = [309000, 7603000]
521#        i4 = [310000, 7580000]
522        i5 = [307000, 7606000]
523        poly_onslow = [i0, i1, i2, i3, i4, i5]
524
525        # Thevenard Island
526        j0 = [294000, 7629000]
527        j1 = [285000, 7625000]
528        j2 = [294000, 7621000]
529        j3 = [299000, 7625000]
530        poly_thevenard = [j0, j1, j2, j3]
531
532        # med res around onslow
533        l0 = [300000, 7610000]
534        l1 = [285000, 7600000]
535        l2 = [300000, 7597500]
536        l3 = [310000, 7770000] # this one is outside
537#        l3 = [310000, 7630000] # this one is NOT outside
538        l4 = [315000, 7610000]
539        poly_coast = [l0, l1, l2, l3, l4]
540
541        # general coast and local area to onslow region
542        m0 = [270000, 7581000]
543        m1 = [300000, 7591000]
544        m2 = [339000, 7610000]
545        m3 = [330000, 7630000]
546        m4 = [290000, 7640000]
547        m5 = [260000, 7600000]
548        poly_region = [m0, m1, m2, m3, m4, m5]
549
550        # This one sticks outside bounding polygon - should fail
551        interior_regions = [[poly_onslow, 50000], [poly_region, 50000],
552                            [poly_coast, 100000], [poly_thevenard, 100000]]
553        boundary_tags = {'walls': [0,1], 'bom': [2]}
554
555        try:
556            m = create_mesh_from_regions(poly_all,
557                                         boundary_tags,
558                                         10000000,
559                                         interior_regions=interior_regions,
560                                         verbose=False)
561        except:
562            pass
563        else:
564            msg = 'Interior polygon sticking outside bounding polygon should '
565            msg += 'cause an Exception to be raised'
566            raise Exception, msg
567
568    def FIXME_test_create_mesh_with_multiply_tagged_segments(self):
569        '''Test that create_mesh_from_regions fails when
570        segments are listed repeatedly in boundary_tags.
571        '''
572
573        # These are the absolute values
574        min_x = 10
575        min_y = 88
576        polygon = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
577        boundary_tags = {'walls': [0,1], 'bom': [1,2]}
578
579        # This one is inside bounding polygon - should pass
580        inner_polygon = [[800,400], [900,500], [800,600]]
581        interior_regions = [(inner_polygon, 5)]
582        m = create_mesh_from_regions(polygon,
583                                     boundary_tags,
584                                     10000000,
585                                     interior_regions=interior_regions,
586                                     verbose=False)
587
588        # This one sticks outside bounding polygon - should fail
589        inner_polygon = [[800,400], [900,500], [800,600]]
590        interior_regions = [(inner_polygon, 5)]
591
592        m = create_mesh_from_regions(polygon,
593                                     boundary_tags,
594                                     10000000,
595                                     interior_regions=interior_regions)
596        try:
597            m = create_mesh_from_regions(polygon,
598                                         boundary_tags,
599                                         10000000,
600                                         interior_regions=interior_regions)
601        except:
602            pass
603        else:
604            msg = 'Tags are listed repeatedly, but create mesh from regions '
605            msg += 'does not cause an Exception to be raised'
606            raise Exception, msg
607
608           
609    def test_create_mesh_with_segments_out_of_bounds(self):
610        """Test that create_mesh_from_regions fails when a segment is out of bounds.
611        """
612       
613        # These are the absolute values
614        min_x = 10 
615        min_y = 88
616        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
617
618       
619        boundary_tags = {'walls': [0,1], 'bom':[2,3], 'out': [5]}
620
621        # This one is inside bounding polygon - should pass
622        inner_polygon = [[800,400],[900,500],[800,600]]
623       
624        interior_regions = [(inner_polygon, 5)]
625
626
627        try:
628            m = create_mesh_from_regions(polygon,
629                                         boundary_tags,
630                                         10000000,
631                                         interior_regions=interior_regions)
632        except:
633            pass
634        else:
635            msg = 'Tags are listed repeatedly, but create mesh from regions '
636            msg += 'does not cause an Exception to be raised'
637            raise Exception, msg
638           
639                   
640
641
642    def test_create_mesh_from_regions_with_duplicate_verts(self):
643        # These are the absolute values
644        polygon_absolute = [[0.0, 0.0],
645                            [0, 4.0],
646                            [4.0, 4.0],
647                            [4.0, 0.0],
648                            [4.0, 0.0]]
649        x_p = -10
650        y_p = -40
651        zone = 808
652        geo_ref_poly = Geo_reference(zone, x_p, y_p)
653        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
654        boundary_tags = {'50': [0],
655                         '40': [1],
656                         '30': [2],
657                         'no where seg': [3],
658                         '20': [4]}
659        m = create_mesh_from_regions(polygon,
660                                     boundary_tags,
661                                     10000000,
662                                     poly_geo_reference=geo_ref_poly,
663                                     verbose=False)
664
665        fileName = 'badmesh.tsh'
666        #m.export_mesh_file(fileName)
667
668    def concept_create_mesh_from_regions_with_ungenerate(self):
669        x=0
670        y=0
671        mesh_geo = geo_reference=Geo_reference(56, x, y)
672
673        # These are the absolute values
674        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
675        x_p = -10
676        y_p = -40
677        geo_ref_poly = Geo_reference(56, x_p, y_p)
678        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
679
680        boundary_tags = {'walls': [0,1], 'bom': [2]}
681
682        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
683        inner1_polygon = geo_ref_poly.\
684                            change_points_geo_ref(inner1_polygon_absolute)
685
686        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
687        inner2_polygon = geo_ref_poly.\
688                            change_points_geo_ref(inner2_polygon_absolute)
689
690        max_area = 10000000
691        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
692        m = create_mesh_from_regions(polygon,
693                                     boundary_tags,
694                                     max_area,
695                                     interior_regions=interior_regions,
696                                     poly_geo_reference=geo_ref_poly,
697                                     mesh_geo_reference=mesh_geo)
698
699        m.export_mesh_file('a_test_mesh_iknterface.tsh')
700
701        fileName = tempfile.mktemp('.txt')
702        file = open(fileName, 'w')
703        file.write('         1       ??      ??\n\
704       90.0       90.0\n\
705       81.0       90.0\n\
706       81.0       81.0\n\
707       90.0       81.0\n\
708       90.0       90.0\n\
709END\n\
710         2      ?? ??\n\
711       10.0       80.0\n\
712       10.0       90.0\n\
713       20.0       90.0\n\
714       10.0       80.0\n\
715END\n\
716END\n')
717        file.close()
718
719        m.import_ungenerate_file(fileName, tag='wall')
720        os.remove(fileName)
721        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
722        m.export_mesh_file('b_test_mesh_iknterface.tsh')
723
724    def concept_ungenerateII(self):
725        from anuga.shallow_water import Domain, Reflective_boundary, \
726                            Dirichlet_boundary
727
728        x=0
729        y=0
730        mesh_geo = geo_reference=Geo_reference(56, x, y)
731
732        # These are the absolute values
733        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
734        x_p = -10
735        y_p = -40
736        geo_ref_poly = Geo_reference(56, x_p, y_p)
737        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
738
739        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
740
741        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
742        inner1_polygon = geo_ref_poly.\
743                            change_points_geo_ref(inner1_polygon_absolute)
744
745        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
746        inner2_polygon = geo_ref_poly.\
747                            change_points_geo_ref(inner2_polygon_absolute)
748
749        max_area = 1
750        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
751        m = create_mesh_from_regions(polygon,
752                                     boundary_tags,
753                                     max_area,
754                                     interior_regions=interior_regions,
755                                     poly_geo_reference=geo_ref_poly,
756                                     mesh_geo_reference=mesh_geo)
757
758        m.export_mesh_file('a_test_mesh_iknterface.tsh')
759
760        fileName = tempfile.mktemp('.txt')
761        file = open(fileName, 'w')
762        file.write('         1       ??      ??\n\
763       90.0       90.0\n\
764       81.0       90.0\n\
765       81.0       81.0\n\
766       90.0       81.0\n\
767       90.0       90.0\n\
768END\n\
769         2      ?? ??\n\
770       10.0       80.0\n\
771       10.0       90.0\n\
772       20.0       90.0\n\
773       10.0       80.0\n\
774END\n\
775END\n')
776        file.close()
777
778        m.import_ungenerate_file(fileName)      #, tag='wall')
779        os.remove(fileName)
780        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
781        mesh_filename = 'bento_b.tsh'
782        m.export_mesh_file(mesh_filename)
783
784        domain = Domain(mesh_filename, use_cache = False)
785
786        Br = Reflective_boundary(domain)
787        Bd = Dirichlet_boundary([3, 0, 0])
788        domain.set_boundary({'wall': Br, 'wave': Bd})
789        yieldstep = 0.1
790        finaltime = 10
791        for t in domain.evolve(yieldstep, finaltime):
792            domain.write_time()
793
794    def concept_ungenerateIII(self):
795        from anuga.shallow_water import Domain, Reflective_boundary, \
796                            Dirichlet_boundary
797        from anuga.pmesh.mesh_interface import create_mesh_from_regions
798
799        # These are the absolute values
800        polygon = [[0,0], [100,0], [100,100], [0,100]]
801
802        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
803        inner1_polygon = [[10,10], [20,10], [20,20], [10,20]]
804        inner2_polygon = [[30,30], [40,30], [40,40], [30,40]]
805
806        max_area = 1
807        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
808        m = create_mesh_from_regions(polygon,
809                                     boundary_tags,
810                                     max_area,
811                                     interior_regions=interior_regions)
812
813        fileName = tempfile.mktemp('.txt')
814        file = open(fileName, 'w')
815        file.write('         1       ??      ??\n\
816       90.0       90.0\n\
817       81.0       90.0\n\
818       81.0       81.0\n\
819       90.0       81.0\n\
820       90.0       90.0\n\
821END\n\
822         2      ?? ??\n\
823       10.0       80.0\n\
824       10.0       90.0\n\
825       20.0       90.0\n\
826       10.0       80.0\n\
827END\n\
828END\n')
829        file.close()
830
831        m.import_ungenerate_file(fileName)
832        os.remove(fileName)
833        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
834        mesh_filename = 'mesh.tsh'
835        m.export_mesh_file(mesh_filename)
836
837        domain = Domain(mesh_filename, use_cache=False)
838
839        Br = Reflective_boundary(domain)
840        Bd = Dirichlet_boundary([3, 0, 0])
841        domain.set_boundary({'wall': Br, 'wave': Bd})
842        yieldstep = 0.1
843        finaltime = 10
844        for t in domain.evolve(yieldstep, finaltime):
845            domain.write_time()
846
847    def test_create_mesh_from_regions_check_segs(self):
848        '''Test that create_mesh_from_regions fails
849        when an interior region is outside bounding polygon.
850        '''
851
852        # These are the absolute values
853        min_x = 10
854        min_y = 88
855        polygon = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
856        boundary_tags = {'walls': [0,1,3], 'bom': [2]}
857
858        # This one is inside bounding polygon - should pass
859        inner_polygon = [[800,400], [900,500], [800,600]]
860        interior_regions = [(inner_polygon, 5)]
861        m = create_mesh_from_regions(polygon,
862                                     boundary_tags,
863                                     10000000,
864                                     interior_regions=interior_regions)
865
866        boundary_tags = {'walls': [0,1,3,4], 'bom': [2]}
867        try:
868            m = create_mesh_from_regions(polygon,
869                                         boundary_tags,
870                                         10000000,
871                                         interior_regions=interior_regions)
872        except:
873            pass
874        else:
875            msg = 'Segment out of bounds not caught '
876            raise Exception, msg
877
878################################################################################
879
880if __name__ == "__main__":
881    suite = unittest.makeSuite(TestCase,'test')
882    runner = unittest.TextTestRunner() #verbosity=2)
883    runner.run(suite)
884
Note: See TracBrowser for help on using the repository browser.