source: anuga_work/development/Rudy_vandrie_2007/create_example_version1.py @ 4440

"""Create mesh and time boundary for Hydrograph example
"""


from Numeric import array, zeros, Float, allclose

from anuga.pmesh.mesh import *
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.coordinate_transforms.geo_reference import Geo_reference
from anuga.geospatial_data import Geospatial_data

import project


def prepare_timeboundary(filename, hydrograph, field):
    """Convert .out file to NetCDF tms file.
    This is a 'throw-away' code taylor made for this type of file

    hydrograph is an identifier e.g.
    HYDROGRAPHS_SUB17       

    field is one of the column headers.

    Time will always be read and converted into seconds.

    The file has contents like

    #####START_HYDROGRAPHS_SUB17       
     Time    Rain  Rainperv    Qtop      Qbot      Qper      Qimp    Qinto_OS   Qout_OS  Stage
      0.0    0.00    0.00     0.000     0.000     0.000     0.000     0.000     0.000     0.000
      5.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
     10.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
     15.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000

    
    """

    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import array

    assert filename[-4:] == '.tms'
   
    outfilename = filename
    infilename = filename[:-4] + '.out'

    print 'Creating', outfilename

    # Read the ascii (.txt) version of this file
    fid = open(infilename)

    # Read all lines and search for selected hydrograph
    lines = fid.readlines()
    fid.close()

    found = False
    search_string = 'START_'+hydrograph
    print 'looking for', search_string
    for i, line in enumerate(lines):
        if line.startswith('####'):
            if line.find(search_string) >= 0:
                print 'Found', line
                found = True
                break

    if found is False:
        msg = 'Did not find', hydrograph
        raise Exception, msg

    # Select data column
    headers = lines[i+1]
    fields = headers.split()
    index = fields.index(field)
    print 'Found header #%d: "%s"' %(index, fields[index])


    # Read data
    search_string = 'END_'+hydrograph
    time = []
    data = []    
    for line in lines[i+2:]:
        if line.startswith('####'):
            if line.find(search_string) >= 0:
                break        
        
        fields = line.split()
        time.append(float(fields[0]))
        data.append(float(fields[index]))


    # Convert to NetCDF
    N = len(time)
    T = array(time, Float)*60  # Time (seconds)
    Q = array(data, Float)     # Values (m^3/s)

    inflow_stage = 0           # Assumed stage at inflow
    boundary_width = 3.402     # Use width to convert flow into x-momentum 

    # Create tms NetCDF file

    fid = NetCDFFile(outfilename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for Benchmark 2'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', Float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', Float, ('number_of_timesteps',))
    fid.variables['stage'][:] = inflow_stage

    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] = Q/boundary_width

    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
    fid.variables['ymomentum'][:] = 0.0

    fid.close()


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


    # Prepare time boundary
    prepare_timeboundary(filename=project.boundary_filename,
                         hydrograph='HYDROGRAPHS_SUB1',
                         field='Qinto_OS')


    #--------------------------------------------------------------------------
    # Create the triangular mesh based on overall clipping polygon with a
    # tagged
    # boundary and interior regions defined in project.py along with
    # resolutions (maximal area of per triangle) for each polygon
    #--------------------------------------------------------------------------


    base_resolution = 1 # Use this to coarsen or refine entire mesh. 

    # Basic geometry and bounding polygon
    xleft   = 0
    xright  = 5.448
    ybottom = 0
    ytop    = 3.402

    point_sw = [xleft, ybottom]
    point_se = [xright, ybottom]
    point_nw = [xleft, ytop]    
    point_ne = [xright, ytop]

    bounding_polygon = [point_se,
                        point_ne,
                        point_nw,
                        point_sw]
    

    # Localised refined area for gulleys
    xl = 4.8
    xr = 5.3
    yb = 1.6
    yt = 2.3
    p0 = [xl, yb]
    p1 = [xr, yb]
    p2 = [xr, yt]
    p3 = [xl, yt]
    
    gulleys = [p0, p1, p2, p3]
    


    # Island area and drawdown region
    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
    island_7 = [xl-.01, yt]
 
    island = [island_0, island_1, island_2,
              island_3, island_4, island_5,
              island_6, island_7]



    # Region spanning half right hand side of domain just inside boundary
    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
    rhs_se = [xright-0.1, ybottom+0.2]
    rhs_ne = [xright-0.1, ytop-0.2]        

    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]

    
    # Interior regions and creation of mesh
    interior_regions = [[rhs_region, 0.0005],
                        [island, 0.0003*base_resolution],
                        [gulleys, 0.00003*base_resolution]]    

    meshname = project.mesh_filename + '.msh'
    m = create_mesh_from_regions(bounding_polygon,
                                 boundary_tags={'wall': [0, 1, 3],
                                                'wave': [2]},     
                                 maximum_triangle_area=0.1*base_resolution,
                                 interior_regions=interior_regions,
                                 filename=project.mesh_filename,
                                 use_cache=False,
                                 verbose=True)