source: misc/inundation-numpy-branch/pmesh/mesh_interface.py @ 4151

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

comments change

File size: 5.3 KB
RevLine 
[2276]1
2#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[3943]3# Assume that the directory anuga_core/source is included in your pythonpath
[2276]4#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5
[2402]6from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
[3514]7from anuga.utilities.polygon import populate_polygon
8from anuga.utilities.numerical_tools import ensure_numeric
[2402]9from Numeric import Float
[2397]10
11# This is due to pmesh being a package and a module and
[2282]12# the current dir being unknown
13try:
14    from pmesh.mesh import Mesh
15except ImportError: 
16    from mesh import Mesh
[2402]17
18
[2276]19def create_mesh_from_regions(bounding_polygon,
20                             boundary_tags,
21                             maximum_triangle_area,
22                             filename=None,
23                             interior_regions=None,
24                             poly_geo_reference=None,
25                             mesh_geo_reference=None,
26                             minimum_triangle_angle=28.0):
27    """Create mesh from bounding polygons, and resolutions.
28
[2402]29    bounding_polygon is a list of points in Eastings and Northings,
30    relative to the poly_geo_reference.
[2276]31
32    Boundary tags is a dictionary of symbolic tags. For every tag there
[2391]33    is a list of indices referring to segments associated with that tag
[2276]34
35    Resolution is the maximal area per triangle for the bounding polygon
36    (excluding interior regions, see later)
37
38    Interior_regions is a list of tuples consisting of (polygon, resolution)
39    for each region to be separately refined.
40
41    This function does not allow segments to share points - use underlying
42    pmesh functionality for that
43
44    poly_geo_reference is the geo_reference of the polygons.
[2402]45    If none, assume absolute.  Please pass one though, since absolute
46    references have a zone.
[2276]47   
48    mesh_geo_reference is the geo_reference of the mesh to be created.
[2402]49    If none is given one will be automatically generated.  It was use
50    the lower left hand corner of  bounding_polygon (absolute)
51    as the x and y values for the geo_ref.
[2276]52   
53    Returns the mesh instance if no finename is given
[2430]54
55    Note, interior regions should be fully nested, as overlaps may cause
56    unintended resolutions.
[2276]57   
58    """
[2402]59    #FIXME (OLE-DSG)
[2276]60    # To do make maximum_triangle_area optional?
61    # check the segment indexes - throw an error if they are out of bounds
[2401]62    #
[2391]63    #In addition I reckon the polygons could be of class Geospatial_data
64
[2402]65    if mesh_geo_reference is None:
66        bounding_polygon = ensure_numeric(bounding_polygon, Float)
67        xllcorner = min(bounding_polygon[:,0])
68        yllcorner = min(bounding_polygon[:,1])   
69        #
70        if poly_geo_reference is None:
71            zone = DEFAULT_ZONE
72        else:
73            zone = poly_geo_reference.get_zone()
74            [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \
75            [(xllcorner,yllcorner)])
76        # create a geo_ref, based on the llc of the bounding_polygon
77        mesh_geo_reference = Geo_reference(xllcorner = xllcorner,
[2408]78                                           yllcorner = yllcorner,
79                                           zone = zone)
[2402]80
[2276]81    m = Mesh(geo_reference=mesh_geo_reference)
82
83    #Do bounding polygon
84    m.add_region_from_polygon(bounding_polygon,
[2391]85                              tags=boundary_tags,
86                              geo_reference=poly_geo_reference)
[2276]87
88    #Find one point inside region automatically
89    if interior_regions is not None:
90        excluded_polygons = []       
91        for polygon, res in interior_regions:
92            excluded_polygons.append( polygon )
93    else:
94        excluded_polygons = None
95
96    from Numeric import array
97
98
99    # convert bounding poly to absolute values
100    # this sort of thing can be fixed with the geo_points class
101    if poly_geo_reference is not None:
102        bounding_polygon_absolute = \
103            poly_geo_reference.get_absolute(bounding_polygon)
104    else:
105        bounding_polygon_absolute = bounding_polygon
106       
107    [inner_point] = populate_polygon(bounding_polygon_absolute,
108                                     1, exclude = excluded_polygons)
109    inner = m.add_region(inner_point[0], inner_point[1])
110    inner.setMaxArea(maximum_triangle_area)
111
112    #Do interior regions
113    if interior_regions is not None:   
114        for polygon, res in interior_regions:
115            #polygon = convert_points_from_latlon_to_utm(P, refzone)
116            m.add_region_from_polygon(polygon,
117                                      geo_reference=poly_geo_reference)
118            # convert bounding poly to absolute values
119            if poly_geo_reference is not None:
120                polygon_absolute = \
121                    poly_geo_reference.get_absolute(polygon)
122            else:
123                polygon_absolute = polygon
124            [inner_point] = populate_polygon(polygon_absolute, 1)
125            region = m.add_region(inner_point[0], inner_point[1])
126            region.setMaxArea(res)
[2391]127           
[2276]128    m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle,
[2391]129                    maximum_triangle_area=maximum_triangle_area)
[2276]130
131    if filename is None:
132        return m
133    else:
134        m.export_mesh_file(filename)
[2402]135
136
Note: See TracBrowser for help on using the repository browser.