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

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

Minor bug in test
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,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
674        segments are listed repeatedly in boundary_tags.
675        """
676       
677       
678       
679
680        # These are the absolute values
681        min_x = 10 
682        min_y = 88
683        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
684
685       
686        boundary_tags = {'walls':[0,1],'bom':[1,2]}
687
688        # This one is inside bounding polygon - should pass
689        inner_polygon = [[800,400],[900,500],[800,600]]
690       
691        interior_regions = [(inner_polygon, 5)]
692        m = create_mesh_from_regions(polygon,
693                                     boundary_tags,
694                                     10000000,
695                                     interior_regions=interior_regions,verbose=False)
696
697
698        # This one sticks outside bounding polygon - should fail
699        inner_polygon = [[800,400],[900,500],[800,600]]
700        interior_regions = [(inner_polygon, 5)]
701
702
703
704        try:
705            m = create_mesh_from_regions(polygon,
706                                         boundary_tags,
707                                         10000000,
708                                         interior_regions=interior_regions)
709        except:
710            pass
711        else:
712            msg = 'Tags are listed repeatedly, but create mesh from regions '
713            msg += 'does not cause an Exception to be raised'
714            raise msg
715
716       
717
718
719    def test_create_mesh_from_regions_with_duplicate_verts(self):
720
721        # These are the absolute values
722       
723        polygon_absolute = [[0.0,0.0],
724                            [0,4.0],
725                            [4.0,4.0],
726                            [4.0,0.0],
727                            [4.0,0.0]]
728       
729        x_p = -10
730        y_p = -40
731        zone = 808
732        geo_ref_poly = Geo_reference(zone, x_p, y_p)
733        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
734
735        boundary_tags = {'50':[0],
736                         '40':[1],
737                         '30':[2],
738                         'no where seg':[3],
739                         '20':[4]
740                         }
741       
742        m = create_mesh_from_regions(polygon,
743                                     boundary_tags,
744                                     10000000,
745                                     poly_geo_reference=geo_ref_poly,verbose=False)
746       
747
748        fileName = "badmesh.tsh"
749        #m.export_mesh_file(fileName)
750 
751       
752    def concept_create_mesh_from_regions_with_ungenerate(self):
753        x=0
754        y=0
755        mesh_geo = geo_reference=Geo_reference(56,x,y)
756       
757        # These are the absolute values
758        polygon_absolute = [[0,0],[100,0],[100,100],[0,100]]
759       
760        x_p = -10
761        y_p = -40
762        geo_ref_poly = Geo_reference(56, x_p, y_p)
763        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
764
765        boundary_tags = {'walls':[0,1],'bom':[2]}
766       
767        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
768        inner1_polygon = geo_ref_poly. \
769                         change_points_geo_ref(inner1_polygon_absolute)
770
771        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
772        inner2_polygon = geo_ref_poly. \
773                         change_points_geo_ref(inner2_polygon_absolute)
774       
775        max_area = 10000000
776        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 10)]
777        m = create_mesh_from_regions(polygon,
778                                     boundary_tags,
779                                     max_area,
780                                     interior_regions=interior_regions,
781                                     poly_geo_reference=geo_ref_poly,
782                                     mesh_geo_reference=mesh_geo)
783                   
784        m.export_mesh_file('a_test_mesh_iknterface.tsh')                 
785       
786        fileName = tempfile.mktemp(".txt")
787        file = open(fileName,"w")
788        file.write("         1       ??      ??\n\
789       90.0       90.0\n\
790       81.0       90.0\n\
791       81.0       81.0\n\
792       90.0       81.0\n\
793       90.0       90.0\n\
794END\n\
795         2      ?? ??\n\
796       10.0       80.0\n\
797       10.0       90.0\n\
798       20.0       90.0\n\
799       10.0       80.0\n\
800END\n\
801END\n")
802        file.close() 
803       
804        m.import_ungenerate_file(fileName, tag='wall')
805        os.remove(fileName)
806        m.generate_mesh(maximum_triangle_area=max_area,verbose=False)
807        m.export_mesh_file('b_test_mesh_iknterface.tsh')
808
809    def concept_ungenerateII(self):
810       
811        from anuga.shallow_water import Domain, Reflective_boundary, \
812                            Dirichlet_boundary
813        x=0
814        y=0
815        mesh_geo = geo_reference=Geo_reference(56,x,y)
816       
817        # These are the absolute values
818        polygon_absolute = [[0,0],[100,0],[100,100],[0,100]]
819       
820        x_p = -10
821        y_p = -40
822        geo_ref_poly = Geo_reference(56, x_p, y_p)
823        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
824
825        boundary_tags = {'wall':[0,1,3],'wave':[2]}
826       
827        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
828        inner1_polygon = geo_ref_poly. \
829                         change_points_geo_ref(inner1_polygon_absolute)
830
831        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
832        inner2_polygon = geo_ref_poly. \
833                         change_points_geo_ref(inner2_polygon_absolute)
834       
835        max_area = 1
836        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 10)]
837        m = create_mesh_from_regions(polygon,
838                                     boundary_tags,
839                                     max_area,
840                                     interior_regions=interior_regions,
841                                     poly_geo_reference=geo_ref_poly,
842                                     mesh_geo_reference=mesh_geo)
843                   
844        m.export_mesh_file('a_test_mesh_iknterface.tsh')                 
845       
846        fileName = tempfile.mktemp(".txt")
847        file = open(fileName,"w")
848        file.write("         1       ??      ??\n\
849       90.0       90.0\n\
850       81.0       90.0\n\
851       81.0       81.0\n\
852       90.0       81.0\n\
853       90.0       90.0\n\
854END\n\
855         2      ?? ??\n\
856       10.0       80.0\n\
857       10.0       90.0\n\
858       20.0       90.0\n\
859       10.0       80.0\n\
860END\n\
861END\n")
862        file.close() 
863       
864        m.import_ungenerate_file(fileName) #, tag='wall')
865        os.remove(fileName)
866        m.generate_mesh(maximum_triangle_area=max_area,verbose=False)
867        mesh_filename = "bento_b.tsh"
868        m.export_mesh_file(mesh_filename)
869
870        domain = Domain(mesh_filename, use_cache = False)
871       
872        Br = Reflective_boundary(domain)
873        Bd = Dirichlet_boundary([3,0,0]) 
874        domain.set_boundary( {'wall': Br, 'wave': Bd} )
875        yieldstep = 0.1
876        finaltime = 10
877        for t in domain.evolve(yieldstep, finaltime):   
878            domain.write_time()
879   
880    def concept_ungenerateIII(self):
881       
882        from anuga.shallow_water import Domain, Reflective_boundary, \
883                            Dirichlet_boundary
884       
885        from anuga.pmesh.mesh_interface import create_mesh_from_regions
886       
887        # These are the absolute values
888        polygon = [[0,0],[100,0],[100,100],[0,100]]
889       
890        boundary_tags = {'wall':[0,1,3],'wave':[2]}
891       
892        inner1_polygon = [[10,10],[20,10],[20,20],[10,20]]
893       
894
895        inner2_polygon = [[30,30],[40,30],[40,40],[30,40]]
896       
897       
898        max_area = 1
899        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 10)]
900        m = create_mesh_from_regions(polygon,
901                                     boundary_tags,
902                                     max_area,
903                                     interior_regions=interior_regions)
904                   
905        fileName = tempfile.mktemp(".txt")
906        file = open(fileName,"w")
907        file.write("         1       ??      ??\n\
908       90.0       90.0\n\
909       81.0       90.0\n\
910       81.0       81.0\n\
911       90.0       81.0\n\
912       90.0       90.0\n\
913END\n\
914         2      ?? ??\n\
915       10.0       80.0\n\
916       10.0       90.0\n\
917       20.0       90.0\n\
918       10.0       80.0\n\
919END\n\
920END\n")
921        file.close() 
922       
923        m.import_ungenerate_file(fileName) 
924        os.remove(fileName)
925        m.generate_mesh(maximum_triangle_area=max_area,verbose=False)
926        mesh_filename = "mesh.tsh"
927        m.export_mesh_file(mesh_filename)
928
929        domain = Domain(mesh_filename, use_cache = False)
930       
931        Br = Reflective_boundary(domain)
932        Bd = Dirichlet_boundary([3,0,0]) 
933        domain.set_boundary( {'wall': Br, 'wave': Bd} )
934        yieldstep = 0.1
935        finaltime = 10
936        for t in domain.evolve(yieldstep, finaltime):   
937            domain.write_time()
938
939           
940       
941    def test_create_mesh_from_regions_check_segs(self):
942        """Test that create_mesh_from_regions fails when an interior region is
943         outside bounding polygon.       """
944       
945
946        # These are the absolute values
947        min_x = 10 
948        min_y = 88
949        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
950       
951        boundary_tags = {'walls':[0,1,3],'bom':[2]}
952#        boundary_tags = {'walls':[0,1]}
953        # This one is inside bounding polygon - should pass
954        inner_polygon = [[800,400],[900,500],[800,600]]
955
956        interior_regions = [(inner_polygon, 5)]
957        m = create_mesh_from_regions(polygon,
958                                     boundary_tags,
959                                     10000000,
960                                     interior_regions=interior_regions)
961
962        boundary_tags = {'walls':[0,1,3,4],'bom':[2]}
963       
964        try:
965            m = create_mesh_from_regions(polygon,
966                                         boundary_tags,
967                                         10000000,
968                                         interior_regions=interior_regions)
969        except:
970            pass
971        else:
972            msg = 'segment out of bounds not caught '
973            raise msg
974
975       
976#-------------------------------------------------------------
977if __name__ == "__main__":
978    suite = unittest.makeSuite(TestCase,'test')
979    #suite = unittest.makeSuite(TestCase,'test_create_mesh_from_regions_check_segs')
980    runner = unittest.TextTestRunner() #verbosity=2)
981    runner.run(suite)
982   
Note: See TracBrowser for help on using the repository browser.