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

Last change on this file since 3734 was 3734, checked in by ole, 18 years ago

Small fix ensuring that mesh file is always stored irrespective of caching.

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, res 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, res 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.