source: trunk/anuga_core/source/anuga/pmesh/mesh_interface.py @ 8124

Last change on this file since 8124 was 8124, checked in by wilsonr, 12 years ago

Made changes described in ticket 359.

File size: 13.1 KB
Line 
1
2from anuga.coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
3from anuga.geometry.polygon import  point_in_polygon ,populate_polygon
4from anuga.utilities.numerical_tools import ensure_numeric
5import numpy as num
6from anuga.geometry.polygon import inside_polygon
7from anuga.geometry.polygon import polylist2points_verts
8import anuga.utilities.log as log
9
10
11# This is due to pmesh being a package and a module and
12# the current dir being unknown
13try:
14    from anuga.pmesh.mesh import Mesh
15except ImportError: 
16    from mesh import Mesh
17
18import exceptions
19class PolygonError(exceptions.Exception): pass
20class SegmentError(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                             hole_tags=None,
29                             poly_geo_reference=None,
30                             mesh_geo_reference=None,
31                             minimum_triangle_angle=28.0,
32                             fail_if_polygons_outside=True,
33                             breaklines=None,
34                             use_cache=False,
35                             verbose=True):
36    """Create mesh from bounding polygons, and resolutions.
37
38    bounding_polygon is a list of points in Eastings and Northings,
39    relative to the poly_geo_reference.
40
41    Boundary tags is a dictionary of symbolic tags. For every tag there
42    is a list of indices referring to segments associated with that tag.
43    If a segment is omitted an Exception will be raised.
44
45    maximum_triangle_area is the maximal area per triangle
46    for the bounding polygon, excluding the  interior regions.
47
48    Interior_regions is a list of tuples consisting of (polygon,
49    resolution) for each region to be separately refined. Do not have
50    polygon lines cross or be on-top of each other.  Also do not have
51    polygon close to each other.
52   
53    NOTE: If a interior_region is outside the bounding_polygon it should
54    throw an error
55   
56    Interior_holes is a list of polygons for each hole.
57    hole_tags is an optional list of boundary tags for the holes, see
58                boundary_tags parameter.
59
60    This function does not allow segments to share points - use underlying
61    pmesh functionality for that
62
63    poly_geo_reference is the geo_reference of the bounding polygon and
64    the interior polygons.
65    If none, assume absolute.  Please pass one though, since absolute
66    references have a zone.
67   
68    mesh_geo_reference is the geo_reference of the mesh to be created.
69    If none is given one will be automatically generated.  It was use
70    the lower left hand corner of  bounding_polygon (absolute)
71    as the x and y values for the geo_ref.
72
73    breaklines is a list of polygons. These lines will be preserved by the
74               triangulation algorithm - useful for coastlines, walls, etc.
75               The polygons are not closed.                       
76   
77    Returns the mesh instance if no filename is given
78
79    Note, interior regions should be fully nested, as overlaps may cause
80    unintended resolutions.
81
82    fail_if_polygons_outside: If True (the default) Exception in thrown
83    where interior polygons fall outside bounding polygon. If False, these
84    will be ignored and execution continued.
85       
86   
87    """
88   
89    # Build arguments and keyword arguments for use with caching or apply.
90    args = (bounding_polygon,
91            boundary_tags)
92   
93    kwargs = {'maximum_triangle_area': maximum_triangle_area,
94              'filename': filename,
95              'interior_regions': interior_regions,
96              'interior_holes': interior_holes,
97              'hole_tags': hole_tags,
98              'poly_geo_reference': poly_geo_reference,
99              'mesh_geo_reference': mesh_geo_reference,
100              'minimum_triangle_angle': minimum_triangle_angle,
101              'fail_if_polygons_outside': fail_if_polygons_outside,
102              'breaklines': breaklines,
103              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day. See ticket:14
104
105    # Call underlying engine with or without caching
106    if use_cache is True:
107        try:
108            from anuga.caching import cache
109        except:
110            msg = 'Caching was requested, but caching module'+\
111                  'could not be imported'
112            raise Exception(msg)
113
114
115        m = cache(_create_mesh_from_regions,
116                  args, kwargs,
117                  verbose=verbose,
118                  compression=False)
119    else:
120        m = apply(_create_mesh_from_regions,
121                  args, kwargs)
122
123    return m   
124       
125
126
127def _create_mesh_from_regions(bounding_polygon,
128                              boundary_tags,
129                              maximum_triangle_area=None,
130                              filename=None,                             
131                              interior_regions=None,
132                              interior_holes=None,
133                              hole_tags=None,
134                              poly_geo_reference=None,
135                              mesh_geo_reference=None,
136                              minimum_triangle_angle=28.0,
137                              fail_if_polygons_outside=True,
138                              breaklines=None,
139                              verbose=True):
140    """_create_mesh_from_regions - internal function.
141
142    See create_mesh_from_regions for documentation.
143    """
144
145   
146    # check the segment indexes - throw an error if they are out of bounds
147    if boundary_tags is not None:
148        max_points = len(bounding_polygon)
149        for key in boundary_tags.keys():
150            if len([x for x in boundary_tags[key] if x > max_points-1]) >= 1:
151                msg = 'Boundary tag %s has segment out of bounds. '\
152                      %(str(key))
153                msg += 'Number of points in bounding polygon = %d' % max_points
154                raise SegmentError(msg)
155
156        for i in range(max_points):
157            found = False
158            for tag in boundary_tags:
159                if i in boundary_tags[tag]:
160                    found = True
161            if found is False:
162                msg = 'Segment %d was not assigned a boundary_tag.' % i
163                msg +=  'Default tag "exterior" will be assigned to missing segment'
164                #raise Exception(msg)
165                # Fixme: Use proper Python warning
166                if verbose: log.critical('WARNING: %s' % msg)
167               
168
169   
170    #In addition I reckon the polygons could be of class Geospatial_data
171    #(DSG) If polygons were classes caching would break in places.
172
173    # Simple check
174    bounding_polygon = ensure_numeric(bounding_polygon, num.float)
175    msg = 'Bounding polygon must be a list of points or an Nx2 array'
176    assert len(bounding_polygon.shape) == 2, msg
177    assert bounding_polygon.shape[1] == 2, msg   
178
179    #
180    if interior_regions is not None:
181       
182        # Test that all the interior polygons are inside the
183        # bounding_poly and throw out those that aren't fully
184        # included.  #Note, Both poly's have the same geo_ref,
185        # therefore don't take into account # geo_ref
186
187
188        polygons_inside_boundary = []
189        for interior_polygon, res in interior_regions:
190            indices = inside_polygon(interior_polygon, bounding_polygon,
191                                     closed = True, verbose = False)
192   
193            if len(indices) <> len(interior_polygon): 
194                msg = 'Interior polygon %s is not fully inside'\
195                      %(str(interior_polygon))
196                msg += ' bounding polygon: %s.' %(str(bounding_polygon))
197
198                if fail_if_polygons_outside is True:
199                    raise PolygonError(msg)
200                else:
201                    msg += ' I will ignore it.'
202                    log.critical(msg)
203
204            else:
205                polygons_inside_boundary.append([interior_polygon, res])
206               
207        # Record only those that were fully contained       
208        interior_regions = polygons_inside_boundary
209
210           
211   
212# the following segment of code could be used to Test that all the
213# interior polygons are inside the bounding_poly... however it might need
214# to be change a bit   
215#
216#count = 0
217#for i in range(len(interior_regions)):
218#    region = interior_regions[i]
219#    interior_polygon = region[0]
220#    if len(inside_polygon(interior_polygon, bounding_polygon,
221#                   closed = True, verbose = False)) <> len(interior_polygon):
222#        print 'WARNING: interior polygon %d is outside bounding polygon' %(i)
223#        count += 1
224
225#if count == 0:
226#    print 'interior regions OK'
227#else:
228#    print 'check out your interior polygons'
229#    print 'check %s in production directory' %figname
230#    import sys; sys.exit()
231
232    if interior_holes is not None:       
233        # Test that all the interior polygons are inside the bounding_poly
234        for interior_polygon in interior_holes:
235            indices = inside_polygon(interior_polygon, bounding_polygon,
236                                     closed = True, verbose = False)
237   
238            if len(indices) <> len(interior_polygon): 
239                msg = 'Interior polygon %s is outside bounding polygon: %s'\
240                      %(str(interior_polygon), str(bounding_polygon))
241                raise PolygonError(msg)
242
243    # Resolve geo referencing       
244    if mesh_geo_reference is None:
245        xllcorner = min(bounding_polygon[:,0])
246        yllcorner = min(bounding_polygon[:,1])   
247        #
248        if poly_geo_reference is None:
249            zone = DEFAULT_ZONE
250        else:
251            zone = poly_geo_reference.get_zone()
252            [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \
253            [(xllcorner,yllcorner)])
254        # create a geo_ref, based on the llc of the bounding_polygon
255        mesh_geo_reference = Geo_reference(xllcorner = xllcorner,
256                                           yllcorner = yllcorner,
257                                           zone = zone)
258
259    m = Mesh(geo_reference=mesh_geo_reference)
260
261    # build a list of discrete segments from the breakline polygons
262    if breaklines is not None:
263        points, verts = polylist2points_verts(breaklines)
264        m.add_points_and_segments(points, verts)
265
266    # Do bounding polygon
267    m.add_region_from_polygon(bounding_polygon,
268                              segment_tags=boundary_tags,
269                              geo_reference=poly_geo_reference)
270
271    # Find one point inside region automatically
272    if interior_regions is not None:
273        excluded_polygons = []       
274        for polygon, res in interior_regions:
275            excluded_polygons.append( polygon )
276    else:
277        excluded_polygons = None
278
279    # Convert bounding poly to absolute values
280    # this sort of thing can be fixed with the geo_points class
281    if poly_geo_reference is not None:
282        bounding_polygon_absolute = \
283            poly_geo_reference.get_absolute(bounding_polygon)
284    else:
285        bounding_polygon_absolute = bounding_polygon
286   
287    inner_point = point_in_polygon(bounding_polygon_absolute)
288    inner = m.add_region(inner_point[0], inner_point[1])
289    inner.setMaxArea(maximum_triangle_area)
290
291    # Do interior regions
292#    if interior_regions is not None:   
293#        for polygon, res in interior_regions:
294#            m.add_region_from_polygon(polygon,
295#                                      geo_reference=poly_geo_reference)
296#            # convert bounding poly to absolute values
297#            if poly_geo_reference is not None:
298#                polygon_absolute = \
299#                    poly_geo_reference.get_absolute(polygon)
300#            else:
301#                polygon_absolute = polygon
302#            inner_point = point_in_polygon(polygon_absolute)
303#            region = m.add_region(inner_point[0], inner_point[1])
304#            region.setMaxArea(res)
305           
306           
307    if interior_regions is not None:   
308        for polygon, res in interior_regions:
309            m.add_region_from_polygon(polygon,
310                                      max_triangle_area=res,
311                                      geo_reference=poly_geo_reference)
312   
313           
314    # Do interior holes
315    if interior_holes is not None:   
316        for n, polygon in enumerate(interior_holes):
317            try:
318                tags = hole_tags[n]
319            except:
320                tags = {}
321            m.add_hole_from_polygon(polygon,
322                                    segment_tags=tags,
323                                    geo_reference=poly_geo_reference)
324       
325
326
327    # NOTE (Ole): This was moved here as it is annoying if mesh is always
328    # stored irrespective of whether the computation
329    # was cached or not. This caused Domain to
330    # recompute as it has meshfile as a dependency
331
332    # Decide whether to store this mesh or return it   
333    if filename is None:
334        return m
335    else:
336        if verbose: log.critical("Generating mesh to file '%s'" % filename)
337        m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle,
338                        verbose=verbose)
339        m.export_mesh_file(filename)
340
341
Note: See TracBrowser for help on using the repository browser.