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

Last change on this file since 3734 was 3734, checked in by ole, 18 years ago

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