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

Last change on this file since 3737 was 3737, checked in by steve, 17 years ago

In island.py the creap is very thin!

File size: 9.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                             interior_holes=None,
28                             poly_geo_reference=None,
29                             mesh_geo_reference=None,
30                             minimum_triangle_angle=28.0,
31                             use_cache=False,
32                             verbose=True):
33    """Create mesh from bounding polygons, and resolutions.
34
35    bounding_polygon is a list of points in Eastings and Northings,
36    relative to the poly_geo_reference.
37
38    Boundary tags is a dictionary of symbolic tags. For every tag there
39    is a list of indices referring to segments associated with that tag
40
41    maximum_triangle_area is the maximal area per triangle
42    for the bounding polygon, excluding the  interior regions.
43
44    Interior_regions is a list of tuples consisting of (polygon, resolution)
45    for each region to be separately refined.
46   
47    Interior_holes is a list of ploygons for each hole.
48
49    This function does not allow segments to share points - use underlying
50    pmesh functionality for that
51
52    poly_geo_reference is the geo_reference of the bounding polygon and
53    the interior polygons.
54    If none, assume absolute.  Please pass one though, since absolute
55    references have a zone.
56   
57    mesh_geo_reference is the geo_reference of the mesh to be created.
58    If none is given one will be automatically generated.  It was use
59    the lower left hand corner of  bounding_polygon (absolute)
60    as the x and y values for the geo_ref.
61   
62    Returns the mesh instance if no finename is given
63
64    Note, interior regions should be fully nested, as overlaps may cause
65    unintended resolutions.
66   
67    """
68
69    # Build arguments and keyword arguments for use with caching or apply.
70    args = (bounding_polygon,
71            boundary_tags)
72   
73    kwargs = {'maximum_triangle_area': maximum_triangle_area,
74              'interior_regions': interior_regions,
75              'interior_holes': interior_holes,
76              'poly_geo_reference': poly_geo_reference,
77              'mesh_geo_reference': mesh_geo_reference,
78              'minimum_triangle_angle': minimum_triangle_angle,
79              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day
80
81
82    #print 'kwargs', kwargs
83
84    # Call underlying engine with or without caching
85    if use_cache is True:
86        try:
87            from caching import cache
88        except:
89            msg = 'Caching was requested, but caching module'+\
90                  'could not be imported'
91            raise msg
92
93
94        m = cache(_create_mesh_from_regions,
95                  args, kwargs,
96                  verbose=verbose,
97                  compression=False)
98    else:
99        m = apply(_create_mesh_from_regions,
100                  args, kwargs)
101
102
103    # Decide whether to store this mesh or return it   
104    if filename is None:
105        return m
106    else:
107        if verbose: print 'Generating mesh to file "%s"' %filename
108        m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle,
109                        verbose=verbose)
110        m.export_mesh_file(filename)
111       
112
113
114def _create_mesh_from_regions(bounding_polygon,
115                              boundary_tags,
116                              maximum_triangle_area=None,
117                              interior_regions=None,
118                              interior_holes=None,
119                              poly_geo_reference=None,
120                              mesh_geo_reference=None,
121                              minimum_triangle_angle=28.0,
122                              verbose=True):
123    """_create_mesh_from_regions - internal function.
124
125    See create_mesh_from_regions for documentation.
126    """
127    #FIXME (OLE-DSG)
128    # check the segment indexes - throw an error if they are out of bounds
129    #(DSG) Yes!
130
131   
132    #In addition I reckon the polygons could be of class Geospatial_data
133    #(DSG) Yes!
134
135    # First check that interior polygons are fully contained in bounding
136    # polygon
137    #Note, Both poly's have the same geo_ref, therefore don't take into account
138    # geo_ref
139
140    # Simple check
141    bounding_polygon = ensure_numeric(bounding_polygon, Float)
142    msg = 'Bounding polygon must be a list of points or an Nx2 array'
143    assert len(bounding_polygon.shape) == 2, msg
144    assert bounding_polygon.shape[1] == 2, msg   
145
146
147    #
148    if interior_regions is not None:       
149        # Test that all the interior polygons are inside the bounding_poly
150        for interior_polygon, res in interior_regions:
151            indices = inside_polygon(interior_polygon, bounding_polygon,
152                                     closed = True, verbose = False)
153   
154            if len(indices) <> len(interior_polygon): 
155                msg = 'Interior polygon %s is outside bounding polygon: %s'\
156                      %(str(interior_polygon), str(bounding_polygon))
157                raise PolygonError, msg
158
159    if interior_holes is not None:       
160        # Test that all the interior polygons are inside the bounding_poly
161        for interior_polygon in interior_holes:
162            indices = inside_polygon(interior_polygon, bounding_polygon,
163                                     closed = True, verbose = False)
164   
165            if len(indices) <> len(interior_polygon): 
166                msg = 'Interior polygon %s is outside bounding polygon: %s'\
167                      %(str(interior_polygon), str(bounding_polygon))
168                raise PolygonError, msg
169
170    # Resolve geo referencing       
171    if mesh_geo_reference is None:
172        xllcorner = min(bounding_polygon[:,0])
173        yllcorner = min(bounding_polygon[:,1])   
174        #
175        if poly_geo_reference is None:
176            zone = DEFAULT_ZONE
177        else:
178            zone = poly_geo_reference.get_zone()
179            [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \
180            [(xllcorner,yllcorner)])
181        # create a geo_ref, based on the llc of the bounding_polygon
182        mesh_geo_reference = Geo_reference(xllcorner = xllcorner,
183                                           yllcorner = yllcorner,
184                                           zone = zone)
185
186    m = Mesh(geo_reference=mesh_geo_reference)
187
188    # Do bounding polygon
189    m.add_region_from_polygon(bounding_polygon,
190                              tags=boundary_tags,
191                              geo_reference=poly_geo_reference)
192
193    # Find one point inside region automatically
194    if interior_regions is not None:
195        excluded_polygons = []       
196        for polygon, res in interior_regions:
197            excluded_polygons.append( polygon )
198    else:
199        excluded_polygons = None
200
201    # Convert bounding poly to absolute values
202    # this sort of thing can be fixed with the geo_points class
203    if poly_geo_reference is not None:
204        bounding_polygon_absolute = \
205            poly_geo_reference.get_absolute(bounding_polygon)
206    else:
207        bounding_polygon_absolute = bounding_polygon
208   
209    inner_point = point_in_polygon(bounding_polygon_absolute)
210    inner = m.add_region(inner_point[0], inner_point[1])
211    inner.setMaxArea(maximum_triangle_area)
212
213    # Do interior regions
214#    if interior_regions is not None:   
215#        for polygon, res in interior_regions:
216#            m.add_region_from_polygon(polygon,
217#                                      geo_reference=poly_geo_reference)
218#            # convert bounding poly to absolute values
219#            if poly_geo_reference is not None:
220#                polygon_absolute = \
221#                    poly_geo_reference.get_absolute(polygon)
222#            else:
223#                polygon_absolute = polygon
224#            inner_point = point_in_polygon(polygon_absolute)
225#            region = m.add_region(inner_point[0], inner_point[1])
226#            region.setMaxArea(res)
227           
228           
229    if interior_regions is not None:   
230        for polygon, res in interior_regions:
231            m.add_region_from_polygon(polygon,
232                                      max_triangle_area = res,
233                                      geo_reference=poly_geo_reference)
234   
235           
236    # Do interior holes
237    if interior_holes is not None:   
238        for polygon in interior_holes:
239            m.add_hole_from_polygon(polygon,
240                                      geo_reference=poly_geo_reference)
241       
242           
243    return m       
244
Note: See TracBrowser for help on using the repository browser.