source: inundation/pmesh/mesh_interface.py @ 3236

Last change on this file since 3236 was 3196, checked in by duncan, 19 years ago

adding tests and comments

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