source: development/cairns_2006/create_mesh.py @ 3434

"""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 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 cairns_project as project
    m.export_mesh_file(project.mesh_filename)