#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # Assume that the directory anuga_core/source is included in your pythonpath #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE from anuga.utilities.polygon import populate_polygon from anuga.utilities.numerical_tools import ensure_numeric from Numeric import Float # This is due to pmesh being a package and a module and # the current dir being unknown try: from pmesh.mesh import Mesh except ImportError: from mesh import Mesh def create_mesh_from_regions(bounding_polygon, boundary_tags, maximum_triangle_area, filename=None, interior_regions=None, poly_geo_reference=None, mesh_geo_reference=None, minimum_triangle_angle=28.0): """Create mesh from bounding polygons, and resolutions. bounding_polygon is a list of points in Eastings and Northings, relative to the poly_geo_reference. Boundary tags is a dictionary of symbolic tags. For every tag there is a list of indices referring to segments associated with that tag Resolution is the maximal area per triangle for the bounding polygon (excluding interior regions, see later) Interior_regions is a list of tuples consisting of (polygon, resolution) for each region to be separately refined. This function does not allow segments to share points - use underlying pmesh functionality for that poly_geo_reference is the geo_reference of the polygons. If none, assume absolute. Please pass one though, since absolute references have a zone. mesh_geo_reference is the geo_reference of the mesh to be created. If none is given one will be automatically generated. It was use the lower left hand corner of bounding_polygon (absolute) as the x and y values for the geo_ref. Returns the mesh instance if no finename is given Note, interior regions should be fully nested, as overlaps may cause unintended resolutions. """ #FIXME (OLE-DSG) # To do make maximum_triangle_area optional? # check the segment indexes - throw an error if they are out of bounds # #In addition I reckon the polygons could be of class Geospatial_data if mesh_geo_reference is None: bounding_polygon = ensure_numeric(bounding_polygon, Float) xllcorner = min(bounding_polygon[:,0]) yllcorner = min(bounding_polygon[:,1]) # if poly_geo_reference is None: zone = DEFAULT_ZONE else: zone = poly_geo_reference.get_zone() [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \ [(xllcorner,yllcorner)]) # create a geo_ref, based on the llc of the bounding_polygon mesh_geo_reference = Geo_reference(xllcorner = xllcorner, yllcorner = yllcorner, zone = zone) m = Mesh(geo_reference=mesh_geo_reference) #Do bounding polygon m.add_region_from_polygon(bounding_polygon, tags=boundary_tags, geo_reference=poly_geo_reference) #Find one point inside region automatically if interior_regions is not None: excluded_polygons = [] for polygon, res in interior_regions: excluded_polygons.append( polygon ) else: excluded_polygons = None from Numeric import array # convert bounding poly to absolute values # this sort of thing can be fixed with the geo_points class if poly_geo_reference is not None: bounding_polygon_absolute = \ poly_geo_reference.get_absolute(bounding_polygon) else: bounding_polygon_absolute = bounding_polygon [inner_point] = populate_polygon(bounding_polygon_absolute, 1, exclude = excluded_polygons) inner = m.add_region(inner_point[0], inner_point[1]) inner.setMaxArea(maximum_triangle_area) #Do interior regions if interior_regions is not None: for polygon, res in interior_regions: #polygon = convert_points_from_latlon_to_utm(P, refzone) m.add_region_from_polygon(polygon, geo_reference=poly_geo_reference) # convert bounding poly to absolute values if poly_geo_reference is not None: polygon_absolute = \ poly_geo_reference.get_absolute(polygon) else: polygon_absolute = polygon [inner_point] = populate_polygon(polygon_absolute, 1) region = m.add_region(inner_point[0], inner_point[1]) region.setMaxArea(res) m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle, maximum_triangle_area=maximum_triangle_area) if filename is None: return m else: m.export_mesh_file(filename)