Index: /production/onslow_2006/project.py
===================================================================
--- /production/onslow_2006/project.py	(revision 2436)
+++ /production/onslow_2006/project.py	(revision 2436)
@@ -0,0 +1,77 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+Also includes origin for slump scenario.
+"""
+
+from os import sep, environ
+from os.path import expanduser
+from utilities.polygon import read_polygon
+import sys
+
+from pmesh.create_mesh import convert_points_from_latlon_to_utm
+                
+#Making assumptions about the location of scenario data
+scenario_dir_name = 'onslow_tsunami_scenario_2006'
+
+"""
+# 250m data to be provided
+coarsename = 'onshore_250' # get from Neil/Ingo (DEM or topo data)
+# 30m data to be provided
+finename = 'onshore_30' # get from Neil/Ingo (DEM or topo data) 
+"""
+
+# clipping region for fine elevation data
+eastingmin = 250000
+eastingmax = 305000
+northingmin = 7580000
+northingmax = 7635000
+
+#swollen/ all data output
+basename = 'source'
+
+if sys.platform == 'win32':
+    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
+else:    
+    home = expanduser('~')
+
+#Derive subdirectories and filenames
+meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
+datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
+outputdir = home+sep+scenario_dir_name+sep+'output'+sep
+polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
+
+meshname = meshdir + basename
+"""
+coarsedemname = datadir + coarsename
+finedemname = datadir + finename
+combineddemname = datadir + 'onslow_combined_elevation'
+"""
+
+outputname = outputdir + basename  #Used by post processing
+
+#!gauge_filename = outputdir + 'onslow_gauges.xya'
+#!gauge_outname = outputdir + 'gauges_max_output.xya'
+
+#Georeferencing
+from coordinate_transforms.redfearn import degminsec2decimal_degrees
+
+refzone = 50
+
+#Main Domain of Onslow: first run NB 21/2/06
+d0 = [305000, 7635000]
+d1 = [280000, 7635000]
+d2 = [250000, 7615000]
+d3 = [250000, 7590000]
+d4 = [310000, 7580000]
+d5 = [330000, 7610000]
+
+polyAll = [d0, d1, d2, d3, d4, d5]
+
+#Interior region - Onslow - take 1
+i0 = [310000, 7610000]
+i1 = [300000, 7610000]
+i2 = [300000, 7600000]
+i3 = [310000, 7600000]
+
+polyInterior = [i0, i1, i2, i3]
+
+#!slump_origin = [385000.0, 6255000.0] #Absolute UTM
Index: /production/onslow_2006/run_onslow.py
===================================================================
--- /production/onslow_2006/run_onslow.py	(revision 2436)
+++ /production/onslow_2006/run_onslow.py	(revision 2436)
@@ -0,0 +1,166 @@
+"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.outputdir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
+"""
+
+
+#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+
+# Related major packages
+from pyvolution.shallow_water import Domain, Reflective_boundary
+from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from pyvolution.combine_pts import combine_rectangular_points_files 
+from pyvolution.pmesh2domain import pmesh_to_domain_instance
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+from smf import slump_tsunami  # Function for submarine mudslide 
+
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+# Do for coarse and fine data
+# Fine pts file to be clipped to area of interest
+#-------------------------------------------------------------------------------
+'''
+# filenames
+coarsedemname = project.coarsedemname
+#finedemname = project.finedemname
+'''
+meshname = project.meshname+'.msh'
+'''
+# coarse data
+convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
+dem2pts(coarsedemname, use_cache=True, verbose=True)
+
+# fine data (clipping the points file to smaller area)
+convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
+dem2pts(finedemname,
+        easting_min=project.eastingmin,
+        easting_max=project.eastingmax,
+        northing_min=project.northingmin,
+        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
+
+
+# combining the coarse and fine data
+combine_rectangular_points_files(project.finedemname + '.pts',
+                                 project.coarsedemname + '.pts',
+                                 project.combineddemname + '.pts')
+'''
+
+#-------------------------------------------------------------------------------                                 
+# 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
+#-------------------------------------------------------------------------------
+
+from pmesh.create_mesh import create_mesh_from_regions
+
+# original
+interior_res = 5000
+interior_regions = [[project.polyInterior, interior_res]]
+
+#FIXME: Fix caching of this one once the mesh_interface is ready
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.polyAll,
+          {'boundary_tags': {'top': [0], 'topleft': [1],
+                             'left': [2], 'bottom': [3],
+                             'bottomright': [4], 'topright': [5]},
+           'resolution': 100000,
+           'filename': meshname,           
+           'interior_regions': interior_regions,
+           'UTM': True,
+           'refzone': project.refzone},
+          verbose = True)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+
+domain = pmesh_to_domain_instance(meshname, Domain,
+                                  use_cache = True,
+                                  verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+
+domain.set_name(project.basename)
+domain.set_datadir(project.outputdir)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+'''
+tsunami_source = slump_tsunami(length=30000.0,
+                               depth=400.0,
+                               slope=6.0,
+                               thickness=176.0, 
+                               radius=3330,
+                               dphi=0.23,
+                               x0=project.slump_origin[0], 
+                               y0=project.slump_origin[1], 
+                               alpha=0.0, 
+                               domain=domain)
+
+'''
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+domain.set_quantity('stage', 0.)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation', 0.
+#                    filename = project.combineddemname + '.pts',
+#                    filename = project.coarsedemname + '.pts',
+#                    use_cache = True,
+#                    verbose = True
+                    )
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+# 10 min square wave starting at 1 min, 6m high
+Bw = Time_boundary(domain=domain,
+                   f=lambda t: [(20<t<200)*6, 0, 0])
+domain.set_boundary( {'top': Bw, 'topleft': Br,
+                             'left': Br, 'bottom': Br,
+                             'bottomright': Br, 'topright': Br} )
+
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 10, finaltime = 300): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'top')     
+    
+print 'That took %.2f seconds' %(time.time()-t0)
