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

Last change on this file since 7317 was 7317, checked in by rwilson, 15 years ago

Replaced 'print' statements with log.critical() calls.

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