[2286] | 1 | """Create mesh for lwru2 validation of the Oshiri island tsunami |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 5 | # Assume that the root AnuGA dir is included in your pythonpath |
---|
| 6 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 7 | |
---|
| 8 | from pmesh.mesh import * |
---|
| 9 | from pyvolution.coordinate_transforms.geo_reference import Geo_reference |
---|
| 10 | from conversion import convert_latlon_to_xycoords |
---|
| 11 | |
---|
| 12 | #------------------------------------------------------------- |
---|
| 13 | if __name__ == "__main__": |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | #Basic geometry |
---|
| 17 | |
---|
| 18 | latlong_origin = 145, -24 |
---|
| 19 | |
---|
| 20 | coord_se = [0.0,0.0] |
---|
| 21 | coord_ne = convert_latlon_to_xycoords(145,-8,latlong_origin) |
---|
| 22 | coord_sw = convert_latlon_to_xycoords(-157,-24,latlong_origin) |
---|
| 23 | coord_nw = convert_latlon_to_xycoords(-157,-8,latlong_origin) |
---|
| 24 | |
---|
| 25 | # Sets origin which is needed by an instance of a mesh object |
---|
| 26 | geo = Geo_reference(xllcorner = coord_se[0], |
---|
| 27 | yllcorner = coord_se[1]) |
---|
| 28 | |
---|
| 29 | print "***********************" |
---|
| 30 | print "geo ref", geo |
---|
| 31 | print "***********************" |
---|
| 32 | |
---|
| 33 | m = Mesh(geo_reference=geo) |
---|
| 34 | |
---|
| 35 | #Boundary |
---|
| 36 | dict = {} |
---|
| 37 | |
---|
| 38 | dict['points'] = [coord_se, |
---|
| 39 | coord_ne, |
---|
| 40 | coord_nw, |
---|
| 41 | coord_sw] |
---|
| 42 | |
---|
| 43 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] #The outer border |
---|
| 44 | |
---|
| 45 | dict['segment_tags'] = ['border', |
---|
| 46 | 'border', |
---|
| 47 | 'border', |
---|
| 48 | 'border'] |
---|
| 49 | |
---|
| 50 | m.addVertsSegs(dict) |
---|
| 51 | |
---|
[2312] | 52 | |
---|
[2286] | 53 | dict = {} |
---|
| 54 | # rupture zone |
---|
| 55 | # Tonga-Kermadec Subduction zone |
---|
| 56 | # located between 180E and -173W and 15S and 35S |
---|
[2312] | 57 | r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin) |
---|
| 58 | r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin) |
---|
| 59 | r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin) |
---|
| 60 | r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin) |
---|
[2286] | 61 | |
---|
[2312] | 62 | dict['points'] = [r0, r1, r2, r3] |
---|
[2286] | 63 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] |
---|
| 64 | dict['segment_tags'] = ['', '', '', ''] |
---|
| 65 | m.addVertsSegs(dict) |
---|
| 66 | |
---|
| 67 | #Interior regions |
---|
| 68 | # Cairns lat and long (16.575S, 145.45E) |
---|
[2312] | 69 | c0 = convert_latlon_to_xycoords(149,-22,latlong_origin) |
---|
| 70 | c1 = convert_latlon_to_xycoords(146,-19,latlong_origin) |
---|
| 71 | c2 = convert_latlon_to_xycoords(145,-15,latlong_origin) |
---|
| 72 | c3 = convert_latlon_to_xycoords(145,-13,latlong_origin) |
---|
| 73 | c4 = convert_latlon_to_xycoords(146,-13,latlong_origin) |
---|
| 74 | c5 = convert_latlon_to_xycoords(147,-17,latlong_origin) |
---|
| 75 | c6 = convert_latlon_to_xycoords(153,-21,latlong_origin) |
---|
| 76 | c7 = convert_latlon_to_xycoords(153,-22,latlong_origin) |
---|
[2286] | 77 | |
---|
[2312] | 78 | dict['points'] = [c0, c1, c2, c3, c4, c5, c6, c7] |
---|
| 79 | dict['segments'] = [[0,1], [1,2], [2,3], [3,4], [4,5], [5,6], [6,7], [7,0]] |
---|
| 80 | dict['segment_tags'] = ['', '', '', '', '', '', '', ''] |
---|
[2286] | 81 | m.addVertsSegs(dict) |
---|
| 82 | |
---|
| 83 | # Fiji Islands |
---|
| 84 | f0 = convert_latlon_to_xycoords(177,-19,latlong_origin) |
---|
| 85 | f1 = convert_latlon_to_xycoords(177,-14,latlong_origin) |
---|
| 86 | f2 = convert_latlon_to_xycoords(-178,-14,latlong_origin) |
---|
| 87 | f3 = convert_latlon_to_xycoords(-178,-19,latlong_origin) |
---|
| 88 | |
---|
| 89 | dict['points'] = [f0, f1, f2, f3] |
---|
| 90 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] |
---|
| 91 | dict['segment_tags'] = ['', '', '', ''] |
---|
| 92 | m.addVertsSegs(dict) |
---|
| 93 | |
---|
| 94 | # Vanuatu and Solomon Islands |
---|
| 95 | v0 = convert_latlon_to_xycoords(169.5,-20.5,latlong_origin) |
---|
| 96 | v1 = convert_latlon_to_xycoords(170,-20.5,latlong_origin) |
---|
| 97 | v2 = convert_latlon_to_xycoords(167,-10,latlong_origin) |
---|
| 98 | v3 = convert_latlon_to_xycoords(165,-10,latlong_origin) |
---|
| 99 | |
---|
| 100 | dict['points'] = [v0, v1, v2, v3] |
---|
| 101 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] |
---|
| 102 | dict['segment_tags'] = ['', '', '', ''] |
---|
| 103 | m.addVertsSegs(dict) |
---|
| 104 | |
---|
| 105 | # Melanesia and New Caledonia |
---|
| 106 | m0 = convert_latlon_to_xycoords(163,-20.5,latlong_origin) |
---|
| 107 | m1 = convert_latlon_to_xycoords(167,-23,latlong_origin) |
---|
| 108 | m2 = convert_latlon_to_xycoords(168.5,-21,latlong_origin) |
---|
| 109 | m3 = convert_latlon_to_xycoords(163,-18,latlong_origin) |
---|
| 110 | |
---|
| 111 | dict['points'] = [m0, m1, m2, m3] |
---|
| 112 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] |
---|
| 113 | dict['segment_tags'] = ['', '', '', ''] |
---|
| 114 | m.addVertsSegs(dict) |
---|
| 115 | |
---|
| 116 | |
---|
[2312] | 117 | #Exclude PNG because land mass has no effect on the wave that hits Cairns |
---|
| 118 | # Papua New Guinea |
---|
| 119 | p0 = convert_latlon_to_xycoords(145.5,-12.0,latlong_origin) |
---|
| 120 | p1 = convert_latlon_to_xycoords(145.5,-8.0,latlong_origin) |
---|
| 121 | p2 = convert_latlon_to_xycoords(154.5,-8.0,latlong_origin) |
---|
| 122 | p3 = convert_latlon_to_xycoords(154.5,-12.0,latlong_origin) |
---|
[2286] | 123 | |
---|
[2312] | 124 | dict['points'] = [p0, p1, p2, p3] |
---|
[2286] | 125 | dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] |
---|
| 126 | dict['segment_tags'] = ['', '', '', ''] |
---|
| 127 | m.addVertsSegs(dict) |
---|
| 128 | |
---|
[2312] | 129 | |
---|
[2286] | 130 | """ |
---|
| 131 | addregionEN - region is defined by a point in the region enclosed |
---|
| 132 | by segments |
---|
| 133 | |
---|
| 134 | setMaxArea - specifies the maximum area of a triangle within a region in |
---|
| 135 | m^2 |
---|
| 136 | """ |
---|
[2312] | 137 | |
---|
| 138 | # base resolution = 5e11 gives good size test grid |
---|
| 139 | #base_resolution = 5e11 |
---|
| 140 | base_resolution = 5e9 |
---|
| 141 | |
---|
[2286] | 142 | ocean = m.addRegionEN(coord_nw[0]-0.1, coord_nw[1]-0.1) |
---|
| 143 | ocean.setMaxArea(base_resolution) |
---|
[2312] | 144 | |
---|
[2286] | 145 | |
---|
[2312] | 146 | rupture = m.addRegionEN(r0[0]+0.1, r0[1]+0.1) |
---|
| 147 | rupture.setMaxArea(100*base_resolution) |
---|
[2286] | 148 | |
---|
| 149 | coast = m.addRegionEN(c0[0]+0.1, c0[1]+0.1) |
---|
[2312] | 150 | coast.setMaxArea(0.003*base_resolution) |
---|
[2286] | 151 | |
---|
| 152 | fiji = m.addRegionEN(f0[0]+0.1, f0[1]+0.1) |
---|
| 153 | fiji.setMaxArea(0.01*base_resolution) |
---|
| 154 | |
---|
| 155 | vanuatu = m.addRegionEN(v0[0]+0.1, v0[1]+0.1) |
---|
| 156 | vanuatu.setMaxArea(0.01*base_resolution) |
---|
| 157 | |
---|
| 158 | melanesia = m.addRegionEN(m0[0]+0.1, m0[1]+0.1) |
---|
[2312] | 159 | melanesia.setMaxArea(0.01*base_resolution) |
---|
[2286] | 160 | |
---|
[2312] | 161 | png = m.addRegionEN(p0[0]+0.1, p0[1]+0.1) |
---|
| 162 | png.setMaxArea(0.01*base_resolution) |
---|
| 163 | |
---|
[2286] | 164 | #a10000a refers to the max area passed to triangle |
---|
| 165 | #28.0 refers to the minimum angle of a triangle |
---|
| 166 | m.generateMesh('pzq28.0z') |
---|
| 167 | |
---|
| 168 | import project |
---|
| 169 | m.export_mesh_file(project.mesh_filename) |
---|
| 170 | |
---|