source: anuga_core/source/anuga/pmesh/test_mesh_interface.py @ 7276

Last change on this file since 7276 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

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