[2276] | 1 | |
---|
| 2 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
[3943] | 3 | # Assume that the directory anuga_core/source is included in your pythonpath |
---|
[2276] | 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 5 | |
---|
[2402] | 6 | from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE |
---|
[3514] | 7 | from anuga.utilities.polygon import populate_polygon |
---|
| 8 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
[2402] | 9 | from Numeric import Float |
---|
[2397] | 10 | |
---|
| 11 | # This is due to pmesh being a package and a module and |
---|
[2282] | 12 | # the current dir being unknown |
---|
| 13 | try: |
---|
| 14 | from pmesh.mesh import Mesh |
---|
| 15 | except ImportError: |
---|
| 16 | from mesh import Mesh |
---|
[2402] | 17 | |
---|
| 18 | |
---|
[2276] | 19 | def create_mesh_from_regions(bounding_polygon, |
---|
| 20 | boundary_tags, |
---|
| 21 | maximum_triangle_area, |
---|
| 22 | filename=None, |
---|
| 23 | interior_regions=None, |
---|
| 24 | poly_geo_reference=None, |
---|
| 25 | mesh_geo_reference=None, |
---|
| 26 | minimum_triangle_angle=28.0): |
---|
| 27 | """Create mesh from bounding polygons, and resolutions. |
---|
| 28 | |
---|
[2402] | 29 | bounding_polygon is a list of points in Eastings and Northings, |
---|
| 30 | relative to the poly_geo_reference. |
---|
[2276] | 31 | |
---|
| 32 | Boundary tags is a dictionary of symbolic tags. For every tag there |
---|
[2391] | 33 | is a list of indices referring to segments associated with that tag |
---|
[2276] | 34 | |
---|
| 35 | Resolution is the maximal area per triangle for the bounding polygon |
---|
| 36 | (excluding interior regions, see later) |
---|
| 37 | |
---|
| 38 | Interior_regions is a list of tuples consisting of (polygon, resolution) |
---|
| 39 | for each region to be separately refined. |
---|
| 40 | |
---|
| 41 | This function does not allow segments to share points - use underlying |
---|
| 42 | pmesh functionality for that |
---|
| 43 | |
---|
| 44 | poly_geo_reference is the geo_reference of the polygons. |
---|
[2402] | 45 | If none, assume absolute. Please pass one though, since absolute |
---|
| 46 | references have a zone. |
---|
[2276] | 47 | |
---|
| 48 | mesh_geo_reference is the geo_reference of the mesh to be created. |
---|
[2402] | 49 | If none is given one will be automatically generated. It was use |
---|
| 50 | the lower left hand corner of bounding_polygon (absolute) |
---|
| 51 | as the x and y values for the geo_ref. |
---|
[2276] | 52 | |
---|
| 53 | Returns the mesh instance if no finename is given |
---|
[2430] | 54 | |
---|
| 55 | Note, interior regions should be fully nested, as overlaps may cause |
---|
| 56 | unintended resolutions. |
---|
[2276] | 57 | |
---|
| 58 | """ |
---|
[2402] | 59 | #FIXME (OLE-DSG) |
---|
[2276] | 60 | # To do make maximum_triangle_area optional? |
---|
| 61 | # check the segment indexes - throw an error if they are out of bounds |
---|
[2401] | 62 | # |
---|
[2391] | 63 | #In addition I reckon the polygons could be of class Geospatial_data |
---|
| 64 | |
---|
[2402] | 65 | if mesh_geo_reference is None: |
---|
| 66 | bounding_polygon = ensure_numeric(bounding_polygon, Float) |
---|
| 67 | xllcorner = min(bounding_polygon[:,0]) |
---|
| 68 | yllcorner = min(bounding_polygon[:,1]) |
---|
| 69 | # |
---|
| 70 | if poly_geo_reference is None: |
---|
| 71 | zone = DEFAULT_ZONE |
---|
| 72 | else: |
---|
| 73 | zone = poly_geo_reference.get_zone() |
---|
| 74 | [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \ |
---|
| 75 | [(xllcorner,yllcorner)]) |
---|
| 76 | # create a geo_ref, based on the llc of the bounding_polygon |
---|
| 77 | mesh_geo_reference = Geo_reference(xllcorner = xllcorner, |
---|
[2408] | 78 | yllcorner = yllcorner, |
---|
| 79 | zone = zone) |
---|
[2402] | 80 | |
---|
[2276] | 81 | m = Mesh(geo_reference=mesh_geo_reference) |
---|
| 82 | |
---|
| 83 | #Do bounding polygon |
---|
| 84 | m.add_region_from_polygon(bounding_polygon, |
---|
[2391] | 85 | tags=boundary_tags, |
---|
| 86 | geo_reference=poly_geo_reference) |
---|
[2276] | 87 | |
---|
| 88 | #Find one point inside region automatically |
---|
| 89 | if interior_regions is not None: |
---|
| 90 | excluded_polygons = [] |
---|
| 91 | for polygon, res in interior_regions: |
---|
| 92 | excluded_polygons.append( polygon ) |
---|
| 93 | else: |
---|
| 94 | excluded_polygons = None |
---|
| 95 | |
---|
| 96 | from Numeric import array |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | # convert bounding poly to absolute values |
---|
| 100 | # this sort of thing can be fixed with the geo_points class |
---|
| 101 | if poly_geo_reference is not None: |
---|
| 102 | bounding_polygon_absolute = \ |
---|
| 103 | poly_geo_reference.get_absolute(bounding_polygon) |
---|
| 104 | else: |
---|
| 105 | bounding_polygon_absolute = bounding_polygon |
---|
| 106 | |
---|
| 107 | [inner_point] = populate_polygon(bounding_polygon_absolute, |
---|
| 108 | 1, exclude = excluded_polygons) |
---|
| 109 | inner = m.add_region(inner_point[0], inner_point[1]) |
---|
| 110 | inner.setMaxArea(maximum_triangle_area) |
---|
| 111 | |
---|
| 112 | #Do interior regions |
---|
| 113 | if interior_regions is not None: |
---|
| 114 | for polygon, res in interior_regions: |
---|
| 115 | #polygon = convert_points_from_latlon_to_utm(P, refzone) |
---|
| 116 | m.add_region_from_polygon(polygon, |
---|
| 117 | geo_reference=poly_geo_reference) |
---|
| 118 | # convert bounding poly to absolute values |
---|
| 119 | if poly_geo_reference is not None: |
---|
| 120 | polygon_absolute = \ |
---|
| 121 | poly_geo_reference.get_absolute(polygon) |
---|
| 122 | else: |
---|
| 123 | polygon_absolute = polygon |
---|
| 124 | [inner_point] = populate_polygon(polygon_absolute, 1) |
---|
| 125 | region = m.add_region(inner_point[0], inner_point[1]) |
---|
| 126 | region.setMaxArea(res) |
---|
[2391] | 127 | |
---|
[2276] | 128 | m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle, |
---|
[2391] | 129 | maximum_triangle_area=maximum_triangle_area) |
---|
[2276] | 130 | |
---|
| 131 | if filename is None: |
---|
| 132 | return m |
---|
| 133 | else: |
---|
| 134 | m.export_mesh_file(filename) |
---|
[2402] | 135 | |
---|
| 136 | |
---|