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

Last change on this file since 3719 was 3719, checked in by steve, 17 years ago

Adding holes to create_mesh_from_regions
File size: 25.4 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
8from anuga.pmesh.mesh_interface import _create_mesh_from_regions
9from load_mesh.loadASCII import *
10from anuga.utilities.polygon import is_inside_polygon
11from anuga.coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
12
13class TestCase(unittest.TestCase):
14    def setUp(self):
15        pass
16   
17    def tearDown(self):
18        pass
19
20    def test_create_mesh_from_regions(self):
21        x=-500
22        y=-1000
23        mesh_geo = geo_reference=Geo_reference(56,x,y)
24       
25        # These are the absolute values
26        polygon_absolute = [[0,0],[100,0],[100,100],[0,100]]
27       
28        x_p = -10
29        y_p = -40
30        geo_ref_poly = Geo_reference(56, x_p, y_p)
31        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
32
33        boundary_tags = {'walls':[0,1],'bom':[2]}
34       
35        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
36        inner1_polygon = geo_ref_poly. \
37                         change_points_geo_ref(inner1_polygon_absolute)
38
39        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
40        inner2_polygon = geo_ref_poly. \
41                         change_points_geo_ref(inner2_polygon_absolute)
42       
43        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
44        m = create_mesh_from_regions(polygon,
45                                     boundary_tags,
46                                     10000000,
47                                     interior_regions=interior_regions,
48                                     poly_geo_reference=geo_ref_poly,
49                                     mesh_geo_reference=mesh_geo)
50
51        # Test the mesh instance
52        self.failUnless(len(m.regions)==3,
53                        'FAILED!')
54        segs = m.getUserSegments()
55        self.failUnless(len(segs)==12,
56                        'FAILED!')
57        self.failUnless(len(m.userVertices)==12,
58                        'FAILED!') 
59        self.failUnless(segs[0].tag=='walls',
60                        'FAILED!') 
61        self.failUnless(segs[1].tag=='walls',
62                        'FAILED!') 
63         
64        self.failUnless(segs[2].tag=='bom',
65                        'FAILED!') 
66        self.failUnless(segs[3].tag=='',
67                        'FAILED!')
68
69        # Assuming the order of the region points is known.
70        # (This isn't true, if you consider create_mesh_from_regions
71        # a black box)
72        poly_point = m.getRegions()[0]
73       
74        #print "poly_point", poly_point
75        #print "polygon_absolute",polygon_absolute
76         
77        # poly_point values are relative to the mesh geo-ref
78        # make them absolute
79        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
80                                          polygon_absolute, closed = False),
81                        'FAILED!')
82       
83        # Assuming the order of the region points is known.
84        # (This isn't true, if you consider create_mesh_from_regions
85        # a black box)
86        poly_point = m.getRegions()[1]
87       
88        #print "poly_point", poly_point
89        #print "polygon_absolute",polygon_absolute
90         
91        # poly_point values are relative to the mesh geo-ref
92        # make them absolute
93        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
94                                          inner1_polygon_absolute,
95                                          closed = False),
96                        'FAILED!')
97       
98        # Assuming the order of the region points is known.
99        # (This isn't true, if you consider create_mesh_from_regions
100        # a black box)
101        poly_point = m.getRegions()[2]
102       
103        #print "poly_point", poly_point
104        #print "polygon_absolute",polygon_absolute
105         
106        # poly_point values are relative to the mesh geo-ref
107        # make them absolute
108        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
109                                          inner2_polygon_absolute,
110                                          closed = False),
111                        'FAILED!')
112
113
114    def test_create_mesh_from_regions_with_caching(self):
115        x=-500
116        y=-1000
117        mesh_geo = geo_reference=Geo_reference(56,x,y)
118       
119        # These are the absolute values
120        polygon_absolute = [[0,0],[100,0],[100,100],[0,100]]
121       
122        x_p = -10
123        y_p = -40
124        geo_ref_poly = Geo_reference(56, x_p, y_p)
125        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
126
127        boundary_tags = {'walls':[0,1],'bom':[2]}
128       
129        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
130        inner1_polygon = geo_ref_poly. \
131                         change_points_geo_ref(inner1_polygon_absolute)
132
133        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
134        inner2_polygon = geo_ref_poly. \
135                         change_points_geo_ref(inner2_polygon_absolute)
136       
137        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
138
139        interior_holes = None
140
141        # Clear cache first
142        from anuga.caching import cache
143        cache(_create_mesh_from_regions,
144              (polygon, boundary_tags),
145              {'minimum_triangle_angle': 28.0,
146               'maximum_triangle_area': 10000000,
147               'filename': None,
148               'interior_regions': interior_regions,
149               'interior_holes': interior_holes,
150               'poly_geo_reference': geo_ref_poly,
151               'mesh_geo_reference': mesh_geo,
152               'verbose': False},
153              verbose=False,
154              clear=1)
155
156       
157        m = create_mesh_from_regions(polygon,
158                                     boundary_tags,
159                                     maximum_triangle_area=10000000,
160                                     interior_regions=interior_regions,
161                                     poly_geo_reference=geo_ref_poly,
162                                     mesh_geo_reference=mesh_geo,
163                                     verbose=False,
164                                     use_cache=True)
165
166
167        # Test the mesh instance
168        self.failUnless(len(m.regions)==3,
169                        'FAILED!')
170        segs = m.getUserSegments()
171        self.failUnless(len(segs)==12,
172                        'FAILED!')
173        self.failUnless(len(m.userVertices)==12,
174                        'FAILED!') 
175        self.failUnless(segs[0].tag=='walls',
176                        'FAILED!') 
177        self.failUnless(segs[1].tag=='walls',
178                        'FAILED!') 
179         
180        self.failUnless(segs[2].tag=='bom',
181                        'FAILED!') 
182        self.failUnless(segs[3].tag=='',
183                        'FAILED!')
184
185        # Assuming the order of the region points is known.
186        # (This isn't true, if you consider create_mesh_from_regions
187        # a black box)
188        poly_point = m.getRegions()[0]
189       
190        #print "poly_point", poly_point
191        #print "polygon_absolute",polygon_absolute
192         
193        # poly_point values are relative to the mesh geo-ref
194        # make them absolute
195        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
196                                          polygon_absolute, closed = False),
197                        'FAILED!')
198       
199        # Assuming the order of the region points is known.
200        # (This isn't true, if you consider create_mesh_from_regions
201        # a black box)
202        poly_point = m.getRegions()[1]
203       
204        #print "poly_point", poly_point
205        #print "polygon_absolute",polygon_absolute
206         
207        # poly_point values are relative to the mesh geo-ref
208        # make them absolute
209        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
210                                          inner1_polygon_absolute,
211                                          closed = False),
212                        'FAILED!')
213       
214        # Assuming the order of the region points is known.
215        # (This isn't true, if you consider create_mesh_from_regions
216        # a black box)
217        poly_point = m.getRegions()[2]
218       
219        #print "poly_point", poly_point
220        #print "polygon_absolute",polygon_absolute
221         
222        # poly_point values are relative to the mesh geo-ref
223        # make them absolute
224        self.failUnless(is_inside_polygon([poly_point.x+x,poly_point.y+y],
225                                          inner2_polygon_absolute,
226                                          closed = False),
227                        'FAILED!')
228
229
230        # Now create m using cached values
231        m_cache = create_mesh_from_regions(polygon,
232                                           boundary_tags,
233                                           10000000,
234                                           interior_regions=interior_regions,
235                                           poly_geo_reference=geo_ref_poly,
236                                           mesh_geo_reference=mesh_geo,
237                                           verbose=False,
238                                           use_cache=True)
239
240
241
242       
243    def test_create_mesh_from_regions2(self):
244
245        # These are the absolute values
246        min_x = -10 
247        min_y = -88
248        polygon_absolute = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
249       
250        x_p = -10
251        y_p = -40
252        zone = 808
253        geo_ref_poly = Geo_reference(zone, x_p, y_p)
254        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
255
256        boundary_tags = {'walls':[0,1],'bom':[2]}
257       
258        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
259        inner1_polygon = geo_ref_poly. \
260                         change_points_geo_ref(inner1_polygon_absolute)
261
262        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
263        inner2_polygon = geo_ref_poly. \
264                         change_points_geo_ref(inner2_polygon_absolute)
265       
266        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
267        m = create_mesh_from_regions(polygon,
268                                     boundary_tags,
269                                     10000000,
270                                     interior_regions=interior_regions,
271                                     poly_geo_reference=geo_ref_poly)
272       
273
274        # Test the mesh instance
275        self.failUnless(len(m.regions)==3,
276                        'FAILED!')
277        segs = m.getUserSegments()
278        self.failUnless(len(segs)==12,
279                        'FAILED!')
280        self.failUnless(len(m.userVertices)==12,
281                        'FAILED!') 
282        self.failUnless(segs[0].tag=='walls',
283                        'FAILED!') 
284        self.failUnless(segs[1].tag=='walls',
285                        'FAILED!') 
286         
287        self.failUnless(segs[2].tag=='bom',
288                        'FAILED!') 
289        self.failUnless(segs[3].tag=='',
290                        'FAILED!')
291       
292        self.failUnless(m.geo_reference.get_zone()==zone,
293                        'FAILED!')
294        self.failUnless(m.geo_reference.get_xllcorner()==min_x,
295                        'FAILED!')
296        self.failUnless(m.geo_reference.get_yllcorner()==min_y,
297                        'FAILED!')
298
299       
300    def test_create_mesh_from_regions3(self):
301
302        # These are the absolute values
303        min_x = -10 
304        min_y = -88
305        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
306       
307
308        x_p = -10
309        y_p = -40
310        geo_ref_poly = Geo_reference(56, x_p, y_p)
311       
312        boundary_tags = {'walls':[0,1],'bom':[2]}
313       
314        inner1_polygon_absolute = [[10,10],[20,10],[20,20],[10,20]]
315        inner1_polygon = geo_ref_poly. \
316                         change_points_geo_ref(inner1_polygon_absolute)
317
318        inner2_polygon_absolute = [[30,30],[40,30],[40,40],[30,40]]
319        inner2_polygon = geo_ref_poly. \
320                         change_points_geo_ref(inner2_polygon_absolute)
321       
322        interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
323        m = create_mesh_from_regions(polygon,
324                                     boundary_tags,
325                                     10000000,
326                                     interior_regions=interior_regions)
327       
328
329        # Test the mesh instance
330        self.failUnless(len(m.regions)==3,
331                        'FAILED!')
332        segs = m.getUserSegments()
333        self.failUnless(len(segs)==12,
334                        'FAILED!')
335        self.failUnless(len(m.userVertices)==12,
336                        'FAILED!') 
337        self.failUnless(segs[0].tag=='walls',
338                        'FAILED!') 
339        self.failUnless(segs[1].tag=='walls',
340                        'FAILED!') 
341         
342        self.failUnless(segs[2].tag=='bom',
343                        'FAILED!') 
344        self.failUnless(segs[3].tag=='',
345                        'FAILED!')
346       
347        self.failUnless(m.geo_reference.get_zone()==DEFAULT_ZONE,
348                        'FAILED!')
349        self.failUnless(m.geo_reference.get_xllcorner()==min_x,
350                        'FAILED!')
351        self.failUnless(m.geo_reference.get_yllcorner()==min_y,
352                        'FAILED!')
353
354    def test_create_mesh_from_regions4(self):
355
356        file_name = tempfile.mktemp(".tsh")
357       
358        # These are the absolute values
359        density_outer = 1000
360        min_outer = 0 
361        max_outer = 1000
362        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
363                   [max_outer,max_outer],[min_outer,max_outer]]
364       
365        density_inner1 = 10000000
366        inner_buffer = 100
367        min_inner1 = min_outer + inner_buffer
368        max_inner1 = max_outer - inner_buffer
369        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
370                   [max_inner1,max_inner1],[min_inner1,max_inner1]]
371     
372       
373        boundary_tags = {'walls':[0,1],'bom':[2]}
374       
375        interior_regions = [(inner1_polygon, density_inner1)]
376        create_mesh_from_regions(polygon_outer
377                                     , boundary_tags
378                                     , density_outer
379                                     , interior_regions=interior_regions
380                                     ,filename=file_name
381                                     #,verbose=True
382                                     ,verbose=False
383                                     )
384       
385        m = importMeshFromFile(file_name)
386       
387        #print "file_name",file_name
388        self.failUnless(len(m.meshTriangles) <= 900,
389                        'Test mesh interface failed!')
390        self.failUnless(len(m.meshTriangles) >= 200,
391                        'Test mesh interface failed!')
392       
393        create_mesh_from_regions(polygon_outer
394                                     , boundary_tags
395                                     , interior_regions=interior_regions
396                                     ,filename=file_name
397                                     #,verbose=True
398                                     ,verbose=False
399                                     )
400       
401        m = importMeshFromFile(file_name)
402       
403        #print "len(m.meshTriangles)",len(m.meshTriangles)
404        self.failUnless(len(m.meshTriangles) <= 100,
405                        'Test mesh interface failed!')
406
407        os.remove(file_name)
408       
409    def test_create_mesh_from_regions5(self):
410
411        file_name = tempfile.mktemp(".tsh")
412       
413        # These are the absolute values
414        density_outer = 10000000 
415        min_outer = 0 
416        max_outer = 1000
417        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
418                   [max_outer,max_outer],[min_outer,max_outer]]
419       
420        density_inner1 = 1000
421        inner_buffer = 100
422        min_inner1 = min_outer + inner_buffer
423        max_inner1 = max_outer - inner_buffer
424        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
425                   [max_inner1,max_inner1],[min_inner1,max_inner1]]
426     
427       
428        boundary_tags = {'walls':[0,1],'bom':[2]}
429       
430        interior_regions = [(inner1_polygon, density_inner1)]
431        create_mesh_from_regions(polygon_outer
432                                     , boundary_tags
433                                     , density_outer
434                                     , interior_regions=interior_regions
435                                     ,filename=file_name
436                                     #,verbose=True
437                                     ,verbose=False
438                                     )
439       
440        m = importMeshFromFile(file_name)
441        #print "file_name",file_name
442        #print "len(m.meshTriangles",len(m.meshTriangles)
443        self.failUnless(len(m.meshTriangles) <= 2000, 
444                        'Test mesh interface failed!')
445 
446        self.failUnless(len(m.meshTriangles) >= 900,
447                        'Test mesh interface failed!')
448
449        os.remove(file_name)
450       
451    def test_create_mesh_from_regions6(self):
452
453        file_name = tempfile.mktemp(".tsh")
454       
455        # These are the absolute values
456        density_outer = 1000
457        min_outer = 0 
458        max_outer = 1000
459        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
460                         [max_outer,max_outer],[min_outer,max_outer]]
461
462        delta = 10
463        density_inner1 = 1000
464        min_inner1 = min_outer + delta
465        max_inner1 = max_outer - delta
466        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
467                          [max_inner1,max_inner1],[min_inner1,max_inner1]]
468     
469       
470        density_inner2 = 10000000 
471        min_inner2 = min_outer +  2*delta
472        max_inner2 = max_outer -  2*delta
473        inner2_polygon = [[min_inner2,min_inner2],[max_inner2,min_inner2],
474                          [max_inner2,max_inner2],[min_inner2,max_inner2]]
475       
476        boundary_tags = {'walls':[0,1],'bom':[2]}
477       
478        interior_regions = [(inner1_polygon, density_inner1),
479                            (inner2_polygon, density_inner2)]
480        create_mesh_from_regions(polygon_outer,
481                                 boundary_tags,
482                                 density_outer,
483                                 interior_regions=interior_regions,
484                                 filename=file_name,
485                                 verbose=False)
486       
487        m = importMeshFromFile(file_name)
488        #print "file_name",file_name
489        #print "len(m.meshTriangles",len(m.meshTriangles)
490        self.failUnless(len(m.meshTriangles) <= 2000, 
491                        'Test mesh interface failed!')
492 
493        self.failUnless(len(m.meshTriangles) >= 900,
494                        'Test mesh interface failed!')
495
496        os.remove(file_name)
497       
498    def test_create_mesh_from_regions7(self):
499
500        file_name = tempfile.mktemp(".tsh")
501       
502        # These are the absolute values
503        density_outer = 1001
504        min_outer = 0 
505        max_outer = 1000
506        polygon_outer = [[min_outer,min_outer],[max_outer,min_outer],
507                   [max_outer,max_outer],[min_outer,max_outer]]
508
509        delta = 10
510        density_inner1 = 100000000
511        min_inner1 = min_outer + delta
512        max_inner1 = max_outer - delta
513        inner1_polygon = [[min_inner1,min_inner1],[max_inner1,min_inner1],
514                   [max_inner1,max_inner1],[min_inner1,max_inner1]]
515     
516       
517        density_inner2 = 1000 
518        min_inner2 = min_outer +  2*delta
519        max_inner2 = max_outer -  2*delta
520        inner2_polygon = [[min_inner2,min_inner2],[max_inner2,min_inner2],
521                   [max_inner2,max_inner2],[min_inner2,max_inner2]]
522       
523        boundary_tags = {'walls':[0,1],'bom':[2]}
524
525        #Note the list order is important
526        # The last region added will be the region triangle uses,
527        # if two regions points are in the same bounded area.
528        interior_regions = [(inner2_polygon, density_inner2),(inner1_polygon, density_inner1)]
529        create_mesh_from_regions(polygon_outer,
530                                 boundary_tags,
531                                 density_outer,
532                                 interior_regions=interior_regions,
533                                 filename=file_name,
534                                 verbose=False)
535       
536        m = importMeshFromFile(file_name)
537        #print "file_name",file_name
538        #print "len(m.meshTriangles",len(m.meshTriangles)
539        self.failUnless(len(m.meshTriangles) <= 3000, 
540                        'Test mesh interface failed!')
541 
542        self.failUnless(len(m.meshTriangles) >= 2000,
543                        'Test mesh interface failed!')
544
545        os.remove(file_name)
546
547       
548    def test_create_mesh_from_regions_interior_regions(self):
549        """Test that create_mesh_from_regions fails when an interior region is
550         outside bounding polygon.       """
551       
552
553        # These are the absolute values
554        min_x = 10 
555        min_y = 88
556        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
557
558       
559        boundary_tags = {'walls':[0,1],'bom':[2]}
560
561        # This one is inside bounding polygon - should pass
562        inner_polygon = [[800,400],[900,500],[800,600]]
563       
564        interior_regions = [(inner_polygon, 5)]
565        m = create_mesh_from_regions(polygon,
566                                     boundary_tags,
567                                     10000000,
568                                     interior_regions=interior_regions)
569
570
571        # This one sticks outside bounding polygon - should fail
572        inner_polygon = [[800,400],[1100,500],[800,600]]
573        interior_regions = [(inner_polygon, 5)]
574
575
576       
577        try:
578            m = create_mesh_from_regions(polygon,
579                                         boundary_tags,
580                                         10000000,
581                                         interior_regions=interior_regions,
582                                         verbose=False)
583        except:
584            pass
585        else:
586            msg = 'Interior polygon sticking outside bounding polygon should '
587            msg += 'cause an Exception to be raised'
588            raise msg
589
590
591
592    def FIXME_test_create_mesh_with_multiply_tagged_segments(self):
593        """Test that create_mesh_from_regions fails when
594        segments are listed repeatedly in boundary_tags.
595        """
596       
597       
598       
599
600        # These are the absolute values
601        min_x = 10 
602        min_y = 88
603        polygon = [[min_x,min_y],[1000,100],[1000,1000],[100,1000]]
604
605       
606        boundary_tags = {'walls':[0,1],'bom':[1,2]}
607
608        # This one is inside bounding polygon - should pass
609        inner_polygon = [[800,400],[900,500],[800,600]]
610       
611        interior_regions = [(inner_polygon, 5)]
612        m = create_mesh_from_regions(polygon,
613                                     boundary_tags,
614                                     10000000,
615                                     interior_regions=interior_regions,verbose=False)
616
617
618        # This one sticks outside bounding polygon - should fail
619        inner_polygon = [[800,400],[900,500],[800,600]]
620        interior_regions = [(inner_polygon, 5)]
621
622
623
624        try:
625            m = create_mesh_from_regions(polygon,
626                                         boundary_tags,
627                                         10000000,
628                                         interior_regions=interior_regions)
629        except:
630            pass
631        else:
632            msg = 'Tags are listed repeatedly, but create mesh from regions '
633            msg += 'does not cause an Exception to be raised'
634            raise msg
635
636       
637
638
639    def test_create_mesh_from_regions_with_duplicate_verts(self):
640
641        # These are the absolute values
642       
643        polygon_absolute = [[0.0,0.0],
644                            [0,4.0],
645                            [4.0,4.0],
646                            [4.0,0.0],
647                            [4.0,0.0]]
648       
649        x_p = -10
650        y_p = -40
651        zone = 808
652        geo_ref_poly = Geo_reference(zone, x_p, y_p)
653        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
654
655        boundary_tags = {'50':[0],
656                         '40':[1],
657                         '30':[2],
658                         'no where seg':[3],
659                         '20':[4]
660                         }
661       
662        m = create_mesh_from_regions(polygon,
663                                     boundary_tags,
664                                     10000000,
665                                     poly_geo_reference=geo_ref_poly,verbose=False)
666       
667
668        fileName = "badmesh.tsh"
669        #m.export_mesh_file(fileName)
670       
671#-------------------------------------------------------------
672if __name__ == "__main__":
673    suite = unittest.makeSuite(TestCase,'test')
674    #suite = unittest.makeSuite(meshTestCase,'test_asciiFile')
675    runner = unittest.TextTestRunner() #verbosity=2)
676    runner.run(suite)
677   
Note: See TracBrowser for help on using the repository browser.