source: anuga_work/production/gippsland_2005/create_mesh2.py @ 4072

"""Create a mesh of the Gippsland Lakes coast based on 100m
   bathymetric data.
"""

# add inundation dir to your pythonpath
from anuga.pmesh.mesh import *
from anuga.coordinate_transforms.geo_reference import Geo_reference

    
def create_mesh(max_area, mesh_file=None, triangles_in_name=False,
                factor = 1000000):

    #Based on a photocopy of the Gippsland area
    p1 = [546000, 5775000]
    p2 = [555000, 5762000]
    p3 = [603000, 5796000]
    p4 = [597000, 5806000]
    p5 = [594000, 5812000]
    p6 = [556000, 5813000]
    p7 = [553000, 5798000]
    p8 = [534000, 5793000]
    
    
    #create geo reference object for use with mesh
    geo = Geo_reference(xllcorner = p2[0],
                        yllcorner = p8[1],
                        zone = 56)
    print "***********************"
    print "geo ref", geo
    print "***********************"
    
    m = Mesh(geo_reference=geo)

    #Boundary
    dict = {}
    dict['points'] = [p1, p2, p3, p4, p5, p6, p7, p8]
    
    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
                        [4,5], [5,6], [6,7], [7,0]] # the outer border
    dict['segment_tags'] = ['side',
                            'ocean',
                            'side',                           
                            'side',
                            'back',
                            'back',
                            'back',                            
                            'side']
    m.addVertsSegs(dict)

     #Based on a photocopy of the Gippsland area
    i1 = [564000, 5794000]
    i2 = [592000, 5807000]
    i3 = [590000, 5809000]
    i4 = [557000, 5805000]

    #Boundary
    dict = {}
    dict['points'] = [i1, i2, i3, i4]
    
    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]] # the inner border
    m.addVertsSegs(dict)


    
    #factor = 10000 #low res 10000, high res 1000
    low = m.addRegionEN(p1[0]+0.5, p1[1]+0.5)
    maxArea = 100*factor
    print "low region maxArea",maxArea 
    low.setMaxArea(maxArea)
    
    medium = m.addRegionEN(i1[0], i1[1]+0.5)
    maxArea = 10*factor
    print "medium region maxArea",maxArea 
    medium.setMaxArea(maxArea)
    medium.setTag('medium')
    
    m.generateMesh('q28.0z', maxArea = max_area)
    triangle_count = len(m.getTriangulation())
    if mesh_file is None:    
        return m, triangle_count
    else:
        if triangles_in_name is True:
            mesh_file = mesh_file[:-4] + '_' + str(triangle_count) \
                        + mesh_file[-4:]
        m.export_mesh_file(mesh_file)
        return mesh_file, triangle_count


#-------------------------------------------------------------
if __name__ == "__main__":
    _, triangle_count = create_mesh(1000000, mesh_file = 'newtest.tsh',
                                    triangles_in_name=True)
    print "triangle_count",triangle_count