| 1 | |
|---|
| 2 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 3 | # Assume that the root AnuGA dir (inundation) is included in your pythonpath |
|---|
| 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 5 | |
|---|
| 6 | from coordinate_transforms.geo_reference import Geo_reference |
|---|
| 7 | # This is due to pmesh not having a class mesh and |
|---|
| 8 | # the current dir being unknown |
|---|
| 9 | try: |
|---|
| 10 | from pmesh.mesh import Mesh |
|---|
| 11 | except ImportError: |
|---|
| 12 | from mesh import Mesh |
|---|
| 13 | from coordinate_transforms.redfearn import redfearn |
|---|
| 14 | from utilities.polygon import populate_polygon |
|---|
| 15 | |
|---|
| 16 | def create_mesh_from_regions(bounding_polygon, |
|---|
| 17 | boundary_tags, |
|---|
| 18 | maximum_triangle_area, |
|---|
| 19 | filename=None, |
|---|
| 20 | interior_regions=None, |
|---|
| 21 | poly_geo_reference=None, |
|---|
| 22 | mesh_geo_reference=None, |
|---|
| 23 | minimum_triangle_angle=28.0): |
|---|
| 24 | """Create mesh from bounding polygons, and resolutions. |
|---|
| 25 | |
|---|
| 26 | Polygon is a list of points in Eastings and Northings, absolute |
|---|
| 27 | |
|---|
| 28 | Boundary tags is a dictionary of symbolic tags. For every tag there |
|---|
| 29 | is a list of |
|---|
| 30 | indices referring to segments associated with that tag |
|---|
| 31 | |
|---|
| 32 | Resolution is the maximal area per triangle for the bounding polygon |
|---|
| 33 | (excluding interior regions, see later) |
|---|
| 34 | |
|---|
| 35 | Interior_regions is a list of tuples consisting of (polygon, resolution) |
|---|
| 36 | for each region to be separately refined. |
|---|
| 37 | |
|---|
| 38 | This function does not allow segments to share points - use underlying |
|---|
| 39 | pmesh functionality for that |
|---|
| 40 | |
|---|
| 41 | poly_geo_reference is the geo_reference of the polygons. |
|---|
| 42 | If none, assume absolute. |
|---|
| 43 | |
|---|
| 44 | mesh_geo_reference is the geo_reference of the mesh to be created. |
|---|
| 45 | If none, assume absolute. |
|---|
| 46 | |
|---|
| 47 | Returns the mesh instance if no finename is given |
|---|
| 48 | |
|---|
| 49 | """ |
|---|
| 50 | # To do make maximum_triangle_area optional? |
|---|
| 51 | # check the segment indexes - throw an error if they are out of bounds |
|---|
| 52 | |
|---|
| 53 | m = Mesh(geo_reference=mesh_geo_reference) |
|---|
| 54 | |
|---|
| 55 | #Do bounding polygon |
|---|
| 56 | m.add_region_from_polygon(bounding_polygon, |
|---|
| 57 | tags=boundary_tags, |
|---|
| 58 | geo_reference=poly_geo_reference) |
|---|
| 59 | |
|---|
| 60 | #Find one point inside region automatically |
|---|
| 61 | if interior_regions is not None: |
|---|
| 62 | excluded_polygons = [] |
|---|
| 63 | for polygon, res in interior_regions: |
|---|
| 64 | #polygon = convert_points_from_latlon_to_utm(P, refzone) |
|---|
| 65 | excluded_polygons.append( polygon ) |
|---|
| 66 | else: |
|---|
| 67 | excluded_polygons = None |
|---|
| 68 | |
|---|
| 69 | from Numeric import array |
|---|
| 70 | |
|---|
| 71 | |
|---|
| 72 | # convert bounding poly to absolute values |
|---|
| 73 | # this sort of thing can be fixed with the geo_points class |
|---|
| 74 | if poly_geo_reference is not None: |
|---|
| 75 | bounding_polygon_absolute = \ |
|---|
| 76 | poly_geo_reference.get_absolute(bounding_polygon) |
|---|
| 77 | else: |
|---|
| 78 | bounding_polygon_absolute = bounding_polygon |
|---|
| 79 | |
|---|
| 80 | [inner_point] = populate_polygon(bounding_polygon_absolute, |
|---|
| 81 | 1, exclude = excluded_polygons) |
|---|
| 82 | inner = m.add_region(inner_point[0], inner_point[1]) |
|---|
| 83 | inner.setMaxArea(maximum_triangle_area) |
|---|
| 84 | |
|---|
| 85 | #Do interior regions |
|---|
| 86 | if interior_regions is not None: |
|---|
| 87 | for polygon, res in interior_regions: |
|---|
| 88 | #polygon = convert_points_from_latlon_to_utm(P, refzone) |
|---|
| 89 | m.add_region_from_polygon(polygon, |
|---|
| 90 | geo_reference=poly_geo_reference) |
|---|
| 91 | # convert bounding poly to absolute values |
|---|
| 92 | if poly_geo_reference is not None: |
|---|
| 93 | polygon_absolute = \ |
|---|
| 94 | poly_geo_reference.get_absolute(polygon) |
|---|
| 95 | else: |
|---|
| 96 | polygon_absolute = polygon |
|---|
| 97 | [inner_point] = populate_polygon(polygon_absolute, 1) |
|---|
| 98 | region = m.add_region(inner_point[0], inner_point[1]) |
|---|
| 99 | region.setMaxArea(res) |
|---|
| 100 | m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle, |
|---|
| 101 | maximum_triangle_area=maximum_triangle_area) |
|---|
| 102 | |
|---|
| 103 | if filename is None: |
|---|
| 104 | return m |
|---|
| 105 | else: |
|---|
| 106 | m.export_mesh_file(filename) |
|---|
