Index: /production/sydney_2006/project.py
===================================================================
--- /production/sydney_2006/project.py	(revision 2181)
+++ /production/sydney_2006/project.py	(revision 2181)
@@ -0,0 +1,69 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+
+from os import sep
+from os.path import expanduser
+import sys
+
+#Making assumptions about the location of scenario data
+scenario_dir_name = 'sydney_tsunami_scenario_2006'
+basename = 'sydney_2006_100' # get from Neil (DEM or topo data)
+
+if sys.platform == 'win32':
+    home = '..\..\..\..'     #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
+#boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
+
+meshname = meshdir + basename
+demname = datadir + basename
+#boundaryname = boundarydir + boundary_basename
+outputname = outputdir + basename  #Used by post processing
+
+
+#Georeferencing
+from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
+
+#Origin of existing dem (FIXME: Temporary measure)
+mesh_origin = (56, 314036, 6224951)  # input from Neil's data
+
+# define clipping polygon
+#south = degminsec2decimal_degrees(-20,45,0)
+#north = degminsec2decimal_degrees(-20,15,0)
+#west = degminsec2decimal_degrees(116,30,0)
+#east = degminsec2decimal_degrees(117,0,0)
+
+south = degminsec2decimal_degrees(-34,05,0)
+north = degminsec2decimal_degrees(-33,45,0)
+west = degminsec2decimal_degrees(151,0,0)
+east = degminsec2decimal_degrees(151,30,0)
+
+p0 = [south, west]
+p1 = [south, east]
+p2 = [north, east]
+p3 = [north, west]
+    
+polygon = [p0, p1, p2, p3]
+refzone = 56
+
+#Interior regions
+harbour_south = degminsec2decimal_degrees(-33,53,0)
+harbour_north = degminsec2decimal_degrees(-33,47,0)
+harbour_west = degminsec2decimal_degrees(151,5,0)
+harbour_east = degminsec2decimal_degrees(151,15,0)
+
+k0 = [harbour_south, harbour_west]
+k1 = [harbour_south, harbour_east]
+k2 = [harbour_north, harbour_east]
+k3 = [harbour_north, harbour_west]    
+
+harbour_polygon = [k0, k1, k2, k3]
+
+
Index: /production/sydney_2006/run_sydney.py
===================================================================
--- /production/sydney_2006/run_sydney.py	(revision 2181)
+++ /production/sydney_2006/run_sydney.py	(revision 2181)
@@ -0,0 +1,113 @@
+"""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 boundary data obtained from a tsunami simulation done with MOST.
+
+Ole Nielsen, GA - 2005
+"""
+
+tide = 0       #Australian Height Datum (mean sea level)
+
+import os
+import time
+
+
+from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
+     Dirichlet_boundary, Time_boundary, Transmissive_boundary
+from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
+     dem2pts
+from pyvolution.pmesh2domain import pmesh_to_domain_instance 
+from caching import cache
+import project
+from math import pi, sin
+
+#Data preparation
+#Convert ASC 2 DEM 2 PTS using source data and store result in source data
+demname = project.demname
+meshname = project.meshname+'.msh'
+
+cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
+      dependencies = [demname + '.asc'],
+      verbose = True)
+      #evaluate = True)
+
+cache(dem2pts, demname, {'verbose': True},
+      dependencies = [demname + '.dem'],      
+      verbose = True)
+
+south = project.south 
+north = project.north
+west = project.west
+east = project.east
+
+#Create Triangular Mesh
+from pmesh.create_mesh import create_mesh_from_regions
+
+interior_regions = [[project.harbour_polygon, 25000]] # maximal area of per triangle
+print project.harbour_polygon
+                   
+m = cache(create_mesh_from_regions,
+          project.polygon,
+          {'boundary_tags': {'bottom': [0], 'top': [2],
+                             'right': [1], 'left': [3]},
+           'resolution': 100000,
+           'filename': meshname,           
+           'interior_regions': interior_regions},     
+          verbose = True)
+
+
+#Setup domain
+
+domain = cache(pmesh_to_domain_instance, (meshname, Domain),
+               dependencies = [meshname],                     
+               verbose = True)               
+
+
+domain.set_name(project.basename)
+domain.set_datadir(project.outputdir)
+domain.store = True
+
+#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+domain.quantities_to_be_stored = ['stage']
+        
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+
+
+
+#Setup Initial Conditions
+domain.set_quantity('friction', 0)
+domain.set_quantity('stage', tide)
+domain.set_quantity('elevation',
+                    filename = demname + '.pts',
+                    use_cache = True,
+                    verbose = True)
+
+
+
+#Setup Boundary Conditions
+print domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bt = Transmissive_boundary(domain)
+Bd = Dirichlet_boundary([0,0,0])
+Bw = Time_boundary(domain=domain,
+                   f=lambda t: [4*sin(2*pi*t/100)+4, 0, 0])
+
+
+domain.set_boundary( {'top': Bd, 'bottom': Bd, 'right': Bd, 'left': Bd} )
+
+
+#Evolve
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 1, finaltime = 4):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'right') #quantities = 'stage')       
+
+print 'That took %.2f seconds' %(time.time()-t0)
