source: anuga_work/production/gippsland_2005/create_mesh2.py @ 5587

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

change imports so reflect the new structure

File size: 2.9 KB
Line 
1"""Create a mesh of the Gippsland Lakes coast based on 100m
2   bathymetric data.
3"""
4
5# add inundation dir to your pythonpath
6from anuga.pmesh.mesh import *
7from anuga.coordinate_transforms.geo_reference import Geo_reference
8
9   
10def create_mesh(max_area, mesh_file=None, triangles_in_name=False,
11                factor = 1000000):
12
13    #Based on a photocopy of the Gippsland area
14    p1 = [546000, 5775000]
15    p2 = [555000, 5762000]
16    p3 = [603000, 5796000]
17    p4 = [597000, 5806000]
18    p5 = [594000, 5812000]
19    p6 = [556000, 5813000]
20    p7 = [553000, 5798000]
21    p8 = [534000, 5793000]
22   
23   
24    #create geo reference object for use with mesh
25    geo = Geo_reference(xllcorner = p2[0],
26                        yllcorner = p8[1],
27                        zone = 56)
28    print "***********************"
29    print "geo ref", geo
30    print "***********************"
31   
32    m = Mesh(geo_reference=geo)
33
34    #Boundary
35    dict = {}
36    dict['points'] = [p1, p2, p3, p4, p5, p6, p7, p8]
37   
38    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
39                        [4,5], [5,6], [6,7], [7,0]] # the outer border
40    dict['segment_tags'] = ['side',
41                            'ocean',
42                            'side',                           
43                            'side',
44                            'back',
45                            'back',
46                            'back',                           
47                            'side']
48    m.addVertsSegs(dict)
49
50     #Based on a photocopy of the Gippsland area
51    i1 = [564000, 5794000]
52    i2 = [592000, 5807000]
53    i3 = [590000, 5809000]
54    i4 = [557000, 5805000]
55
56    #Boundary
57    dict = {}
58    dict['points'] = [i1, i2, i3, i4]
59   
60    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] # the inner border
61    m.addVertsSegs(dict)
62
63
64   
65    #factor = 10000 #low res 10000, high res 1000
66    low = m.addRegionEN(p1[0]+0.5, p1[1]+0.5)
67    maxArea = 100*factor
68    print "low region maxArea",maxArea
69    low.setMaxArea(maxArea)
70   
71    medium = m.addRegionEN(i1[0], i1[1]+0.5)
72    maxArea = 10*factor
73    print "medium region maxArea",maxArea
74    medium.setMaxArea(maxArea)
75    medium.setTag('medium')
76   
77    m.generateMesh('q28.0z', maxArea = max_area)
78    triangle_count = len(m.getTriangulation())
79    if mesh_file is None:   
80        return m, triangle_count
81    else:
82        if triangles_in_name is True:
83            mesh_file = mesh_file[:-4] + '_' + str(triangle_count) \
84                        + mesh_file[-4:]
85        m.export_mesh_file(mesh_file)
86        return mesh_file, triangle_count
87
88
89#-------------------------------------------------------------
90if __name__ == "__main__":
91    _, triangle_count = create_mesh(1000000, mesh_file = 'newtest.tsh',
92                                    triangles_in_name=True)
93    print "triangle_count",triangle_count
Note: See TracBrowser for help on using the repository browser.