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

Last change on this file since 6179 was 6179, checked in by ole, 10 years ago

Fixed up tests reflecting changes in changeset:6177

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