"""Create mesh for lwru2 validation of the Oshiri island tsunami """ #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # Assume that the root AnuGA dir is included in your pythonpath #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! from pmesh.mesh import * from pyvolution.coordinate_transforms.geo_reference import Geo_reference from conversion import convert_latlon_to_xycoords #------------------------------------------------------------- if __name__ == "__main__": #Basic geometry latlong_origin = 145, -24 coord_se = [0.0,0.0] coord_ne = convert_latlon_to_xycoords(145,-8,latlong_origin) coord_sw = convert_latlon_to_xycoords(-157,-24,latlong_origin) coord_nw = convert_latlon_to_xycoords(-157,-8,latlong_origin) # Sets origin which is needed by an instance of a mesh object geo = Geo_reference(xllcorner = coord_se[0], yllcorner = coord_se[1]) print "***********************" print "geo ref", geo print "***********************" m = Mesh(geo_reference=geo) #Boundary dict = {} dict['points'] = [coord_se, coord_ne, coord_nw, coord_sw] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] #The outer border dict['segment_tags'] = ['border', 'border', 'border', 'border'] m.addVertsSegs(dict) dict = {} # rupture zone # Tonga-Kermadec Subduction zone # located between 180E and -173W and 15S and 35S r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin) r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin) r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin) r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin) dict['points'] = [r0, r1, r2, r3] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] dict['segment_tags'] = ['', '', '', ''] m.addVertsSegs(dict) #Interior regions # Cairns lat and long (16.575S, 145.45E) c0 = convert_latlon_to_xycoords(149,-22,latlong_origin) c1 = convert_latlon_to_xycoords(146,-19,latlong_origin) c2 = convert_latlon_to_xycoords(145,-15,latlong_origin) c3 = convert_latlon_to_xycoords(145,-13,latlong_origin) c4 = convert_latlon_to_xycoords(146,-13,latlong_origin) c5 = convert_latlon_to_xycoords(147,-17,latlong_origin) c6 = convert_latlon_to_xycoords(153,-21,latlong_origin) c7 = convert_latlon_to_xycoords(153,-22,latlong_origin) dict['points'] = [c0, c1, c2, c3, c4, c5, c6, c7] dict['segments'] = [[0,1], [1,2], [2,3], [3,4], [4,5], [5,6], [6,7], [7,0]] dict['segment_tags'] = ['', '', '', '', '', '', '', ''] m.addVertsSegs(dict) # Fiji Islands f0 = convert_latlon_to_xycoords(177,-19,latlong_origin) f1 = convert_latlon_to_xycoords(177,-14,latlong_origin) f2 = convert_latlon_to_xycoords(-178,-14,latlong_origin) f3 = convert_latlon_to_xycoords(-178,-19,latlong_origin) dict['points'] = [f0, f1, f2, f3] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] dict['segment_tags'] = ['', '', '', ''] m.addVertsSegs(dict) # Vanuatu and Solomon Islands v0 = convert_latlon_to_xycoords(169.5,-20.5,latlong_origin) v1 = convert_latlon_to_xycoords(170,-20.5,latlong_origin) v2 = convert_latlon_to_xycoords(167,-10,latlong_origin) v3 = convert_latlon_to_xycoords(165,-10,latlong_origin) dict['points'] = [v0, v1, v2, v3] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] dict['segment_tags'] = ['', '', '', ''] m.addVertsSegs(dict) # Melanesia and New Caledonia m0 = convert_latlon_to_xycoords(163,-20.5,latlong_origin) m1 = convert_latlon_to_xycoords(167,-23,latlong_origin) m2 = convert_latlon_to_xycoords(168.5,-21,latlong_origin) m3 = convert_latlon_to_xycoords(163,-18,latlong_origin) dict['points'] = [m0, m1, m2, m3] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] dict['segment_tags'] = ['', '', '', ''] m.addVertsSegs(dict) #Exclude PNG because land mass has no effect on the wave that hits Cairns # Papua New Guinea p0 = convert_latlon_to_xycoords(145.5,-12.0,latlong_origin) p1 = convert_latlon_to_xycoords(145.5,-8.0,latlong_origin) p2 = convert_latlon_to_xycoords(154.5,-8.0,latlong_origin) p3 = convert_latlon_to_xycoords(154.5,-12.0,latlong_origin) dict['points'] = [p0, p1, p2, p3] dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] dict['segment_tags'] = ['', '', '', ''] m.addVertsSegs(dict) """ addregionEN - region is defined by a point in the region enclosed by segments setMaxArea - specifies the maximum area of a triangle within a region in m^2 """ # base resolution = 5e11 gives good size test grid #base_resolution = 5e11 base_resolution = 5e9 ocean = m.addRegionEN(coord_nw[0]-0.1, coord_nw[1]-0.1) ocean.setMaxArea(base_resolution) rupture = m.addRegionEN(r0[0]+0.1, r0[1]+0.1) rupture.setMaxArea(100*base_resolution) coast = m.addRegionEN(c0[0]+0.1, c0[1]+0.1) coast.setMaxArea(0.003*base_resolution) fiji = m.addRegionEN(f0[0]+0.1, f0[1]+0.1) fiji.setMaxArea(0.01*base_resolution) vanuatu = m.addRegionEN(v0[0]+0.1, v0[1]+0.1) vanuatu.setMaxArea(0.01*base_resolution) melanesia = m.addRegionEN(m0[0]+0.1, m0[1]+0.1) melanesia.setMaxArea(0.01*base_resolution) png = m.addRegionEN(p0[0]+0.1, p0[1]+0.1) png.setMaxArea(0.01*base_resolution) #a10000a refers to the max area passed to triangle #28.0 refers to the minimum angle of a triangle m.generateMesh('pzq28.0z') import project m.export_mesh_file(project.mesh_filename)