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

Last change on this file since 6671 was 6360, checked in by rwilson, 16 years ago

Ongoing conversion changes.

File size: 32.8 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        try:
592            m = create_mesh_from_regions(polygon,
593                                         boundary_tags,
594                                         10000000,
595                                         interior_regions=interior_regions)
596        except:
597            pass
598        else:
599            msg = 'Tags are listed repeatedly, but create mesh from regions '
600            msg += 'does not cause an Exception to be raised'
601            raise Exception, msg
602
603    def test_create_mesh_from_regions_with_duplicate_verts(self):
604        # These are the absolute values
605        polygon_absolute = [[0.0, 0.0],
606                            [0, 4.0],
607                            [4.0, 4.0],
608                            [4.0, 0.0],
609                            [4.0, 0.0]]
610        x_p = -10
611        y_p = -40
612        zone = 808
613        geo_ref_poly = Geo_reference(zone, x_p, y_p)
614        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
615        boundary_tags = {'50': [0],
616                         '40': [1],
617                         '30': [2],
618                         'no where seg': [3],
619                         '20': [4]}
620        m = create_mesh_from_regions(polygon,
621                                     boundary_tags,
622                                     10000000,
623                                     poly_geo_reference=geo_ref_poly,
624                                     verbose=False)
625
626        fileName = 'badmesh.tsh'
627        #m.export_mesh_file(fileName)
628
629    def concept_create_mesh_from_regions_with_ungenerate(self):
630        x=0
631        y=0
632        mesh_geo = geo_reference=Geo_reference(56, x, y)
633
634        # These are the absolute values
635        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
636        x_p = -10
637        y_p = -40
638        geo_ref_poly = Geo_reference(56, x_p, y_p)
639        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
640
641        boundary_tags = {'walls': [0,1], 'bom': [2]}
642
643        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
644        inner1_polygon = geo_ref_poly.\
645                            change_points_geo_ref(inner1_polygon_absolute)
646
647        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
648        inner2_polygon = geo_ref_poly.\
649                            change_points_geo_ref(inner2_polygon_absolute)
650
651        max_area = 10000000
652        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
653        m = create_mesh_from_regions(polygon,
654                                     boundary_tags,
655                                     max_area,
656                                     interior_regions=interior_regions,
657                                     poly_geo_reference=geo_ref_poly,
658                                     mesh_geo_reference=mesh_geo)
659
660        m.export_mesh_file('a_test_mesh_iknterface.tsh')
661
662        fileName = tempfile.mktemp('.txt')
663        file = open(fileName, 'w')
664        file.write('         1       ??      ??\n\
665       90.0       90.0\n\
666       81.0       90.0\n\
667       81.0       81.0\n\
668       90.0       81.0\n\
669       90.0       90.0\n\
670END\n\
671         2      ?? ??\n\
672       10.0       80.0\n\
673       10.0       90.0\n\
674       20.0       90.0\n\
675       10.0       80.0\n\
676END\n\
677END\n')
678        file.close()
679
680        m.import_ungenerate_file(fileName, tag='wall')
681        os.remove(fileName)
682        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
683        m.export_mesh_file('b_test_mesh_iknterface.tsh')
684
685    def concept_ungenerateII(self):
686        from anuga.shallow_water import Domain, Reflective_boundary, \
687                            Dirichlet_boundary
688
689        x=0
690        y=0
691        mesh_geo = geo_reference=Geo_reference(56, x, y)
692
693        # These are the absolute values
694        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
695        x_p = -10
696        y_p = -40
697        geo_ref_poly = Geo_reference(56, x_p, y_p)
698        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
699
700        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
701
702        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
703        inner1_polygon = geo_ref_poly.\
704                            change_points_geo_ref(inner1_polygon_absolute)
705
706        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
707        inner2_polygon = geo_ref_poly.\
708                            change_points_geo_ref(inner2_polygon_absolute)
709
710        max_area = 1
711        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
712        m = create_mesh_from_regions(polygon,
713                                     boundary_tags,
714                                     max_area,
715                                     interior_regions=interior_regions,
716                                     poly_geo_reference=geo_ref_poly,
717                                     mesh_geo_reference=mesh_geo)
718
719        m.export_mesh_file('a_test_mesh_iknterface.tsh')
720
721        fileName = tempfile.mktemp('.txt')
722        file = open(fileName, 'w')
723        file.write('         1       ??      ??\n\
724       90.0       90.0\n\
725       81.0       90.0\n\
726       81.0       81.0\n\
727       90.0       81.0\n\
728       90.0       90.0\n\
729END\n\
730         2      ?? ??\n\
731       10.0       80.0\n\
732       10.0       90.0\n\
733       20.0       90.0\n\
734       10.0       80.0\n\
735END\n\
736END\n')
737        file.close()
738
739        m.import_ungenerate_file(fileName)      #, tag='wall')
740        os.remove(fileName)
741        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
742        mesh_filename = 'bento_b.tsh'
743        m.export_mesh_file(mesh_filename)
744
745        domain = Domain(mesh_filename, use_cache = False)
746
747        Br = Reflective_boundary(domain)
748        Bd = Dirichlet_boundary([3, 0, 0])
749        domain.set_boundary({'wall': Br, 'wave': Bd})
750        yieldstep = 0.1
751        finaltime = 10
752        for t in domain.evolve(yieldstep, finaltime):
753            domain.write_time()
754
755    def concept_ungenerateIII(self):
756        from anuga.shallow_water import Domain, Reflective_boundary, \
757                            Dirichlet_boundary
758        from anuga.pmesh.mesh_interface import create_mesh_from_regions
759
760        # These are the absolute values
761        polygon = [[0,0], [100,0], [100,100], [0,100]]
762
763        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
764        inner1_polygon = [[10,10], [20,10], [20,20], [10,20]]
765        inner2_polygon = [[30,30], [40,30], [40,40], [30,40]]
766
767        max_area = 1
768        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
769        m = create_mesh_from_regions(polygon,
770                                     boundary_tags,
771                                     max_area,
772                                     interior_regions=interior_regions)
773
774        fileName = tempfile.mktemp('.txt')
775        file = open(fileName, 'w')
776        file.write('         1       ??      ??\n\
777       90.0       90.0\n\
778       81.0       90.0\n\
779       81.0       81.0\n\
780       90.0       81.0\n\
781       90.0       90.0\n\
782END\n\
783         2      ?? ??\n\
784       10.0       80.0\n\
785       10.0       90.0\n\
786       20.0       90.0\n\
787       10.0       80.0\n\
788END\n\
789END\n')
790        file.close()
791
792        m.import_ungenerate_file(fileName)
793        os.remove(fileName)
794        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
795        mesh_filename = 'mesh.tsh'
796        m.export_mesh_file(mesh_filename)
797
798        domain = Domain(mesh_filename, use_cache=False)
799
800        Br = Reflective_boundary(domain)
801        Bd = Dirichlet_boundary([3, 0, 0])
802        domain.set_boundary({'wall': Br, 'wave': Bd})
803        yieldstep = 0.1
804        finaltime = 10
805        for t in domain.evolve(yieldstep, finaltime):
806            domain.write_time()
807
808    def test_create_mesh_from_regions_check_segs(self):
809        '''Test that create_mesh_from_regions fails
810        when an interior region is outside bounding polygon.
811        '''
812
813        # These are the absolute values
814        min_x = 10
815        min_y = 88
816        polygon = [[min_x,min_y], [1000,100], [1000,1000], [100,1000]]
817        boundary_tags = {'walls': [0,1,3], 'bom': [2]}
818
819        # This one is inside bounding polygon - should pass
820        inner_polygon = [[800,400], [900,500], [800,600]]
821        interior_regions = [(inner_polygon, 5)]
822        m = create_mesh_from_regions(polygon,
823                                     boundary_tags,
824                                     10000000,
825                                     interior_regions=interior_regions)
826
827        boundary_tags = {'walls': [0,1,3,4], 'bom': [2]}
828        try:
829            m = create_mesh_from_regions(polygon,
830                                         boundary_tags,
831                                         10000000,
832                                         interior_regions=interior_regions)
833        except:
834            pass
835        else:
836            msg = 'Segment out of bounds not caught '
837            raise Exception, msg
838
839################################################################################
840
841if __name__ == "__main__":
842    suite = unittest.makeSuite(TestCase,'test')
843    runner = unittest.TextTestRunner() #verbosity=2)
844    runner.run(suite)
845
Note: See TracBrowser for help on using the repository browser.