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

Last change on this file since 3741 was 3741, checked in by duncan, 18 years ago

fix bugs

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