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

Last change on this file since 8279 was 8279, checked in by neweyv, 11 years ago

Additional information written to the log file. This information will assist in creating an algorithm to predict the runtime of ANUGA given various parameters.

File size: 13.4 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
9import datetime
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   
90    if verbose: log.resource_usage_timing(log.logging.INFO, "start_")
91    if verbose: log.timingInfo("maximum_triangle_area, " + str(maximum_triangle_area))
92    if verbose: log.timingInfo("minimum_triangle_angle, " + str(minimum_triangle_angle))
93    if verbose: log.timingInfo("startMesh, '%s'" % log.CurrentDateTime())
94   
95    # Build arguments and keyword arguments for use with caching or apply.
96    args = (bounding_polygon,
97            boundary_tags)
98   
99    kwargs = {'maximum_triangle_area': maximum_triangle_area,
100              'filename': filename,
101              'interior_regions': interior_regions,
102              'interior_holes': interior_holes,
103              'hole_tags': hole_tags,
104              'poly_geo_reference': poly_geo_reference,
105              'mesh_geo_reference': mesh_geo_reference,
106              'minimum_triangle_angle': minimum_triangle_angle,
107              'fail_if_polygons_outside': fail_if_polygons_outside,
108              'breaklines': breaklines,
109              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day. See ticket:14
110
111    # Call underlying engine with or without caching
112    if use_cache is True:
113        try:
114            from anuga.caching import cache
115        except:
116            msg = 'Caching was requested, but caching module'+\
117                  'could not be imported'
118            raise Exception(msg)
119
120
121        m = cache(_create_mesh_from_regions,
122                  args, kwargs,
123                  verbose=verbose,
124                  compression=False)
125    else:
126        m = apply(_create_mesh_from_regions,
127                  args, kwargs)
128
129
130   
131    return m   
132       
133
134
135def _create_mesh_from_regions(bounding_polygon,
136                              boundary_tags,
137                              maximum_triangle_area=None,
138                              filename=None,                             
139                              interior_regions=None,
140                              interior_holes=None,
141                              hole_tags=None,
142                              poly_geo_reference=None,
143                              mesh_geo_reference=None,
144                              minimum_triangle_angle=28.0,
145                              fail_if_polygons_outside=True,
146                              breaklines=None,
147                              verbose=True):
148    """_create_mesh_from_regions - internal function.
149
150    See create_mesh_from_regions for documentation.
151    """
152
153   
154    # check the segment indexes - throw an error if they are out of bounds
155    if boundary_tags is not None:
156        max_points = len(bounding_polygon)
157        for key in boundary_tags.keys():
158            if len([x for x in boundary_tags[key] if x > max_points-1]) >= 1:
159                msg = 'Boundary tag %s has segment out of bounds. '\
160                      %(str(key))
161                msg += 'Number of points in bounding polygon = %d' % max_points
162                raise SegmentError(msg)
163
164        for i in range(max_points):
165            found = False
166            for tag in boundary_tags:
167                if i in boundary_tags[tag]:
168                    found = True
169            if found is False:
170                msg = 'Segment %d was not assigned a boundary_tag.' % i
171                msg +=  'Default tag "exterior" will be assigned to missing segment'
172                #raise Exception(msg)
173                # Fixme: Use proper Python warning
174                if verbose: log.critical('WARNING: %s' % msg)
175               
176
177   
178    #In addition I reckon the polygons could be of class Geospatial_data
179    #(DSG) If polygons were classes caching would break in places.
180
181    # Simple check
182    bounding_polygon = ensure_numeric(bounding_polygon, num.float)
183    msg = 'Bounding polygon must be a list of points or an Nx2 array'
184    assert len(bounding_polygon.shape) == 2, msg
185    assert bounding_polygon.shape[1] == 2, msg   
186
187    #
188    if interior_regions is not None:
189       
190        # Test that all the interior polygons are inside the
191        # bounding_poly and throw out those that aren't fully
192        # included.  #Note, Both poly's have the same geo_ref,
193        # therefore don't take into account # geo_ref
194
195
196        polygons_inside_boundary = []
197        for interior_polygon, res in interior_regions:
198            indices = inside_polygon(interior_polygon, bounding_polygon,
199                                     closed = True, verbose = False)
200   
201            if len(indices) <> len(interior_polygon): 
202                msg = 'Interior polygon %s is not fully inside'\
203                      %(str(interior_polygon))
204                msg += ' bounding polygon: %s.' %(str(bounding_polygon))
205
206                if fail_if_polygons_outside is True:
207                    raise PolygonError(msg)
208                else:
209                    msg += ' I will ignore it.'
210                    log.critical(msg)
211
212            else:
213                polygons_inside_boundary.append([interior_polygon, res])
214               
215        # Record only those that were fully contained       
216        interior_regions = polygons_inside_boundary
217
218           
219   
220# the following segment of code could be used to Test that all the
221# interior polygons are inside the bounding_poly... however it might need
222# to be change a bit   
223#
224#count = 0
225#for i in range(len(interior_regions)):
226#    region = interior_regions[i]
227#    interior_polygon = region[0]
228#    if len(inside_polygon(interior_polygon, bounding_polygon,
229#                   closed = True, verbose = False)) <> len(interior_polygon):
230#        print 'WARNING: interior polygon %d is outside bounding polygon' %(i)
231#        count += 1
232
233#if count == 0:
234#    print 'interior regions OK'
235#else:
236#    print 'check out your interior polygons'
237#    print 'check %s in production directory' %figname
238#    import sys; sys.exit()
239
240    if interior_holes is not None:       
241        # Test that all the interior polygons are inside the bounding_poly
242        for interior_polygon in interior_holes:
243            indices = inside_polygon(interior_polygon, bounding_polygon,
244                                     closed = True, verbose = False)
245   
246            if len(indices) <> len(interior_polygon): 
247                msg = 'Interior polygon %s is outside bounding polygon: %s'\
248                      %(str(interior_polygon), str(bounding_polygon))
249                raise PolygonError(msg)
250
251    # Resolve geo referencing       
252    if mesh_geo_reference is None:
253        xllcorner = min(bounding_polygon[:,0])
254        yllcorner = min(bounding_polygon[:,1])   
255        #
256        if poly_geo_reference is None:
257            zone = DEFAULT_ZONE
258        else:
259            zone = poly_geo_reference.get_zone()
260            [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \
261            [(xllcorner,yllcorner)])
262        # create a geo_ref, based on the llc of the bounding_polygon
263        mesh_geo_reference = Geo_reference(xllcorner = xllcorner,
264                                           yllcorner = yllcorner,
265                                           zone = zone)
266
267    m = Mesh(geo_reference=mesh_geo_reference)
268
269    # build a list of discrete segments from the breakline polygons
270    if breaklines is not None:
271        points, verts = polylist2points_verts(breaklines)
272        m.add_points_and_segments(points, verts)
273
274    # Do bounding polygon
275    m.add_region_from_polygon(bounding_polygon,
276                              segment_tags=boundary_tags,
277                              geo_reference=poly_geo_reference)
278
279    # Find one point inside region automatically
280    if interior_regions is not None:
281        excluded_polygons = []       
282        for polygon, res in interior_regions:
283            excluded_polygons.append( polygon )
284    else:
285        excluded_polygons = None
286
287    # Convert bounding poly to absolute values
288    # this sort of thing can be fixed with the geo_points class
289    if poly_geo_reference is not None:
290        bounding_polygon_absolute = \
291            poly_geo_reference.get_absolute(bounding_polygon)
292    else:
293        bounding_polygon_absolute = bounding_polygon
294   
295    inner_point = point_in_polygon(bounding_polygon_absolute)
296    inner = m.add_region(inner_point[0], inner_point[1])
297    inner.setMaxArea(maximum_triangle_area)
298
299    # Do interior regions
300#    if interior_regions is not None:   
301#        for polygon, res in interior_regions:
302#            m.add_region_from_polygon(polygon,
303#                                      geo_reference=poly_geo_reference)
304#            # convert bounding poly to absolute values
305#            if poly_geo_reference is not None:
306#                polygon_absolute = \
307#                    poly_geo_reference.get_absolute(polygon)
308#            else:
309#                polygon_absolute = polygon
310#            inner_point = point_in_polygon(polygon_absolute)
311#            region = m.add_region(inner_point[0], inner_point[1])
312#            region.setMaxArea(res)
313           
314           
315    if interior_regions is not None:   
316        for polygon, res in interior_regions:
317            m.add_region_from_polygon(polygon,
318                                      max_triangle_area=res,
319                                      geo_reference=poly_geo_reference)
320   
321           
322    # Do interior holes
323    if interior_holes is not None:   
324        for n, polygon in enumerate(interior_holes):
325            try:
326                tags = hole_tags[n]
327            except:
328                tags = {}
329            m.add_hole_from_polygon(polygon,
330                                    segment_tags=tags,
331                                    geo_reference=poly_geo_reference)
332       
333
334
335    # NOTE (Ole): This was moved here as it is annoying if mesh is always
336    # stored irrespective of whether the computation
337    # was cached or not. This caused Domain to
338    # recompute as it has meshfile as a dependency
339
340    # Decide whether to store this mesh or return it
341   
342   
343    if filename is None:
344        return m
345    else:
346        if verbose: log.critical("Generating mesh to file '%s'" % filename)
347     
348        m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle,
349                        verbose=verbose)
350        m.export_mesh_file(filename)
351       
Note: See TracBrowser for help on using the repository browser.