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

Last change on this file since 7244 was 6973, checked in by ole, 16 years ago

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