source: production/karratha_2005/create_mesh.py @ 1796

"""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 pyvolution.coordinate_transforms.redfearn import redfearn


import project

#-------------------------------------------------------------
if __name__ == "__main__":


    
    #Basic geometry

    south = project.south 
    north = project.north
    west = project.west
    east = project.east

    refzone = project.refzone


    #SW
    zone, easting, northing  = redfearn(south, west)
    assert zone == refzone
    point_sw = [easting, northing]

    #SE
    zone, easting, northing  = redfearn(south, east)
    assert zone == refzone
    point_se = [easting, northing]    

    #NW
    zone, easting, northing  = redfearn(north, west)
    assert zone == refzone
    point_nw = [easting, northing]

    #NE
    zone, easting, northing  = redfearn(north, east)
    assert zone == refzone
    point_ne = [easting, northing]    
    

    geo = Geo_reference(xllcorner = 421544.35127423,  #from dem
                        yllcorner = 7677669.5257159,
                        zone = refzone)
    
                        
    print "***********************"
    print "geo ref", geo
    print "***********************"
    
    m = Mesh(geo_reference=geo)

    #Boundary
    dict = {}
    dict['points'] = [point_sw,
                      point_se,
                      point_ne,
                      point_nw]

    
    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
    dict['segment_tags'] = ['wall', 'wall', 'wall', 'wall']
        
    m.addVertsSegs(dict)


    base_resolution = 100000

    inner = m.addRegionEN(point_sw[0]+1, point_sw[1]+1)
    inner.setMaxArea(base_resolution)

    m.generateMesh('pzq28.0za1000000a')

    m.export_mesh_file(project.meshname + '.msh')