source: inundation-numpy-branch/pmesh/create_mesh.py @ 3514

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 4.9 KB
Line 
1"""Create mesh for lwru2 validation of the Oshiri island tsunami
2"""
3
4# Assume that the root AnuGA dir (inundation) is included in your pythonpath
5
6from coordinate_transforms.redfearn import redfearn
7
8def convert_points_from_latlon_to_utm(polygon, refzone=None):
9    points = []
10    for point in polygon:
11        zone, easting, northing  = redfearn(point[0], point[1])
12        #FIXME: Use point.latitude etc once we have a proper point set
13
14        if refzone is None:
15            refzone = zone
16        else:
17            assert zone == refzone
18       
19        points.append([easting, northing])
20
21    return points, zone   
22
23
24
25
26def create_region(polygon, tags, refzone):
27    """Create dictionary representation of region bounded by given polygon for use with pmesh
28    """
29
30
31    #
32    dict = {}
33    dict['points'] = polygon
34
35    #Create segments
36    #E.g. [[0,1], [1,2], [2,3], [3,0]]
37    segments = []
38    N = len(polygon)
39    for i in range(N):
40        lo = i
41        hi = (lo + 1) % N
42        segments.append( [lo, hi] ) 
43       
44       
45    dict['segments'] = segments
46
47
48    #Create tags
49    #E.g. ['wall', 'wall', 'ocean', 'wall']
50
51    if tags is not None:
52        segment_tags = [0]*N
53        for key in tags:
54            indices = tags[key]
55            for i in indices:
56                segment_tags[i] = key
57   
58        dict['segment_tags'] = segment_tags
59
60
61
62    return dict
63
64
65
66def create_mesh_from_regions(bounding_polygon,
67                             boundary_tags,
68                             resolution,
69                             filename = None,
70                             interior_regions = None,
71                             UTM = False,
72                             refzone = None):
73    """Create mesh from bounding polygon, tags for all segments and resolution.
74
75    Polygon is a list of points in latitude and longitudes - decimal degrees
76
77    Boundary tags is a dictionary of symbolic tags. For every tag there is a list of
78    indices referring to segments associated with that tag
79
80    Resolution is the maximal area per triangle for the bounding polygon
81    (excluding interior regions, see later)
82
83    Interior_regions is a list of tuples consisting of (polygon, resolution)
84    for each region to be separately refined.
85
86    This function does not allow segments to share points - use underlying
87    pmesh functionality for that
88
89    FIXME - FUTURE. Something like this.
90      Use class Point_set
91
92      Accept deg, min, sec, e.g.
93      [ [(-20,30,55), (116, 20, 00)], ...]
94
95
96
97   
98    """
99
100    from mesh import Mesh
101    from coordinate_transforms.redfearn import redfearn
102    from anuga.utilities.polygon import populate_polygon
103    from anuga.utilities.numerical_tools import ensure_numeric
104    from coordinate_transforms.geo_reference import Geo_reference
105
106
107    bounding_polygon = ensure_numeric(bounding_polygon)
108    #print 'refzone', refzone
109
110    #Convert to UTM
111    if not UTM:
112        bounding_polygon, refzone = convert_points_from_latlon_to_utm(bounding_polygon)
113    else:
114       assert refzone is not None, 'Refzone must be specified, got %s' %refzone 
115
116
117    #FIXME: Must give georeference here.
118    #Why can't Mesh work out the ll corners?
119
120    xllcorner = min(bounding_polygon[:,0])
121    yllcorner = min(bounding_polygon[:,1])   
122   
123    geo = Geo_reference(xllcorner = xllcorner,
124                        yllcorner = yllcorner,
125                        zone = refzone)
126
127    print 'geo reference derived from bounding polygon', geo
128    m = Mesh(geo_reference=geo)
129   
130   
131   
132    #Do bounding polygon
133    region_dict = create_region(bounding_polygon, boundary_tags, refzone)   
134    m.addVertsSegs(region_dict)
135
136
137    #Find one point inside region automatically
138    if interior_regions is not None:
139        excluded_polygons = []       
140        for P, res in interior_regions:
141            if not UTM:           
142                polygon, _ = convert_points_from_latlon_to_utm(P, refzone)
143            else:
144                polygon = P
145            excluded_polygons.append( polygon )
146    else:
147        excluded_polygons = None
148
149    from Numeric import array
150   
151    [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons)
152    print 'I:', inner_point, resolution
153    inner = m.addRegionEN(inner_point[0], inner_point[1])
154    inner.setMaxArea(resolution)
155
156
157
158    #Do interior regions
159    if interior_regions is not None:   
160        for P, res in interior_regions:
161            if not UTM:                       
162                polygon, _ = convert_points_from_latlon_to_utm(P, refzone)
163            else:
164                polygon = P
165               
166            region_dict = create_region(polygon, None, refzone)
167
168            m.addVertsSegs(region_dict)
169
170            [inner_point] = populate_polygon(polygon, 1)
171            print 'I', inner_point, res
172            X = m.addRegionEN(inner_point[0], inner_point[1])
173            X.setMaxArea(res)
174           
175
176   
177    m.generateMesh('pzq28.0za1000000a')
178   
179    m.export_mesh_file(filename)
180   
181
182
183
184
185   
Note: See TracBrowser for help on using the repository browser.