source: anuga_core/source/anuga/pmesh/mesh_interface.py @ 3689

Last change on this file since 3689 was 3689, checked in by ole, 17 years ago

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.

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