Index: /documentation/user_manual/examples/bedslope.py
===================================================================
--- /documentation/user_manual/examples/bedslope.py	(revision 2481)
+++ /documentation/user_manual/examples/bedslope.py	(revision 2481)
@@ -0,0 +1,58 @@
+"""Simple example of shallow water wave equation using Pyvolution
+
+Water driven by linear slope and Dirichlet boundary
+
+"""
+
+######################
+# Module imports 
+#
+from pyvolution.mesh_factory import rectangular
+from pyvolution.shallow_water import Domain, Reflective_boundary,\
+     Dirichlet_boundary, Time_boundary, Transmissive_boundary
+
+#Create basic triangular mesh
+points, vertices, boundary = rectangular(2, 1)
+
+#Create shallow water domain
+domain = Domain(points, vertices, boundary)
+domain.set_name('bedslope')
+domain.set_datadir('.')                      #Use current directory for output
+domain.set_quantities_to_be_stored('stage')  #See shallow_water.py
+
+
+#######################
+# Initial conditions
+def f(x,y):
+    return -x/2
+
+domain.set_quantity('elevation', f)
+domain.set_quantity('friction', 0.1)
+
+h = 0.05  # Constant depth
+domain.set_quantity('stage', expression = 'elevation + %f' %h)
+
+
+######################
+# Boundary conditions
+from math import sin, pi
+Br = Reflective_boundary(domain)
+Bt = Transmissive_boundary(domain)
+Bd = Dirichlet_boundary([0.2,0.,0.])
+Bw = Time_boundary(domain=domain,
+                   f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
+
+
+print 'Tags are ', domain.get_boundary_tags()
+
+domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+#domain.set_boundary({'left': Bw, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+######################
+#Evolution
+domain.check_integrity()
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 0.2):
+    domain.write_time()
+
Index: /documentation/user_manual/examples/project.py
===================================================================
--- /documentation/user_manual/examples/project.py	(revision 2481)
+++ /documentation/user_manual/examples/project.py	(revision 2481)
@@ -0,0 +1,175 @@
+"""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
+
+# We'll do this for the final thing 
+home = '.'
+
+
+if sys.platform == 'win32':
+    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
+else:    
+    home = expanduser('~')
+    
+
+basename = 'Sydney'
+
+#Setting the location of scenario data
+scenario_dir_name = 'sydney_tsunami_example_2006'
+# revised 100m data
+coarsename = 'bathyland100' # get from Neil/Ingo (DEM or topo data)
+# revised 25m data
+finename = 'bathy_dem25' # get from Neil/Ingo (DEM or topo data) Wed 25 Jan
+
+
+#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 + 'sydneytopo'
+
+#Georeferencing
+from coordinate_transforms.redfearn import degminsec2decimal_degrees
+
+refzone = 56
+print 'refzone', refzone
+
+dsouth = degminsec2decimal_degrees(-34,05,0)
+dnorth = degminsec2decimal_degrees(-33,33,0)
+dnorth1 = degminsec2decimal_degrees(-33,40,0)
+dnorth2 = degminsec2decimal_degrees(-33,58,30)
+dnorth3 = degminsec2decimal_degrees(-33,46,0)
+dwest = degminsec2decimal_degrees(151,2,20)
+deast1 = degminsec2decimal_degrees(151,20,0)
+deast2 = degminsec2decimal_degrees(151,48,0)
+deast3 = degminsec2decimal_degrees(151,10,0)
+deast4 = degminsec2decimal_degrees(151,9,0)
+
+dp0 = [dsouth, dwest]
+dp1 = [dsouth, deast1]
+dp2 = [dnorth2, deast2]
+dp3 = [dnorth1, deast2]
+dp4 = [dnorth, deast1]
+dp5 = [dnorth, deast4]
+dp6 = [dnorth3, deast3]
+dp7 = [dnorth3, dwest]
+dp8 = [dnorth, dwest]
+    
+diffpolygonall, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7])
+
+qp0 = [343494.0, 6270650.0]
+qp1 = [343337.0, 6270303.0]
+qp2 = [343466.0, 6270228.0]
+qp3 = [343139.0, 6269901.0]
+qp4 = [342472.0, 6268573.0]
+qp5 = [342111.0, 6267115.0]
+qp6 = [342479.0, 6266121.0]
+qp7 = [342860.0, 6266386.0]
+qp8 = [342635.0, 6267438.0]
+qp9 = [343105.0, 6269070.0]
+qp10 = [343548.0, 6269567.0]
+qp11 = [343487.0, 6269928.0]
+qp12 = [343991.0, 6270269.0]
+
+#Interior regions - the Harbour
+harbour_1x = degminsec2decimal_degrees(-33,51,0)
+harbour_1y = degminsec2decimal_degrees(151,2,30)
+harbour_12x = degminsec2decimal_degrees(-33,51,0)
+harbour_12y = degminsec2decimal_degrees(151,5,0)
+harbour_13x = degminsec2decimal_degrees(-33,52,15)
+harbour_13y = degminsec2decimal_degrees(151,5,0)
+harbour_2x = degminsec2decimal_degrees(-33,53,0)
+harbour_2y = degminsec2decimal_degrees(151,17,20)
+harbour_3x = degminsec2decimal_degrees(-33,47,0)
+harbour_3y = degminsec2decimal_degrees(151,20,30)
+harbour_4x = degminsec2decimal_degrees(-33,47,50)
+harbour_4y = degminsec2decimal_degrees(151,8,10)
+harbour_5x = degminsec2decimal_degrees(-33,48,10)
+harbour_5y = degminsec2decimal_degrees(151,8,0)
+harbour_6x = degminsec2decimal_degrees(-33,49,0)
+harbour_6y = degminsec2decimal_degrees(151,2,30)
+harbour_7x = degminsec2decimal_degrees(-33,34,30)
+harbour_7y = degminsec2decimal_degrees(151,20,20)
+harbour_8x = degminsec2decimal_degrees(-33,33,30)
+harbour_8y = degminsec2decimal_degrees(151,17,0)
+harbour_9x = degminsec2decimal_degrees(-33,45,30)
+harbour_9y = degminsec2decimal_degrees(151,17,0)
+harbour_10x = degminsec2decimal_degrees(-33,45,10)
+harbour_10y = degminsec2decimal_degrees(151,11,40)
+harbour_11x = degminsec2decimal_degrees(-33,45,10)
+harbour_11y = degminsec2decimal_degrees(151,11,40)
+harbour_14x = degminsec2decimal_degrees(-33,49,10)
+harbour_14y = degminsec2decimal_degrees(151,11,40)
+harbour_15x = degminsec2decimal_degrees(-33,48,55)
+harbour_15y = degminsec2decimal_degrees(151,2,30)
+
+k02 = [harbour_1x, harbour_1y]
+k12 = [harbour_2x, harbour_2y]
+k22 = [harbour_3x, harbour_3y]
+k32 = [harbour_4x, harbour_4y]
+k42 = [harbour_5x, harbour_5y]
+k52 = [harbour_6x, harbour_6y]
+k62 = [harbour_7x, harbour_7y]
+k72 = [harbour_8x, harbour_8y]
+k82 = [harbour_9x, harbour_9y]
+k92 = [harbour_10x, harbour_10y]
+k102 = [harbour_11x, harbour_11y]
+k112 = [harbour_12x, harbour_12y]
+k122 = [harbour_13x, harbour_13y]
+k132 = [harbour_14x, harbour_14y]
+k142 = [harbour_15x, harbour_15y]
+
+harbour_polygon_2, zone = convert_points_from_latlon_to_utm([k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52]) 
+assert zone == refzone
+
+
+#Interior region - Botany Bay
+bb_1x = degminsec2decimal_degrees(-34,3,0)
+bb_1y = degminsec2decimal_degrees(151,2,30)
+bb_10x = degminsec2decimal_degrees(-34,3,0)
+bb_10y = degminsec2decimal_degrees(151,8,0)
+bb_2x = degminsec2decimal_degrees(-34,3,0)
+bb_2y = degminsec2decimal_degrees(151,14,0)
+bb_3x = degminsec2decimal_degrees(-33,53,30)
+bb_3y = degminsec2decimal_degrees(151,17,20)
+bb_4x = degminsec2decimal_degrees(-33,53,0)
+bb_4y = degminsec2decimal_degrees(151,8,0)
+bb_5x = degminsec2decimal_degrees(-33,57,30)
+bb_5y = degminsec2decimal_degrees(151,8,0)
+bb_6x = degminsec2decimal_degrees(-33,57,30)
+bb_6y = degminsec2decimal_degrees(151,2,30)  
+bb_7x = degminsec2decimal_degrees(-33,53,30)
+bb_7y = degminsec2decimal_degrees(151,12,30)  
+bb_8x = degminsec2decimal_degrees(-33,55,20)
+bb_8y = degminsec2decimal_degrees(151,8,0)  
+bb_9x = degminsec2decimal_degrees(-33,55,20)
+bb_9y = degminsec2decimal_degrees(151,12,30)  
+
+j02 = [bb_1x, bb_1y]
+j12 = [bb_2x, bb_2y]
+j22 = [bb_3x, bb_3y]
+j32 = [bb_4x, bb_4y]
+j42 = [bb_5x, bb_5y]
+j52 = [bb_6x, bb_6y]
+j62 = [bb_7x, bb_7y]
+j72 = [bb_8x, bb_8y]
+j82 = [bb_9x, bb_9y]
+j92 = [bb_10x, bb_10y]
+
+botanybay_polygon_2, zone = convert_points_from_latlon_to_utm([j92, j12, j22, j62, j82, j72, j42]) 
+
+
+
+slump_origin = [385000.0, 6255000.0] #Absolute UTM
Index: /documentation/user_manual/examples/runsydney.py
===================================================================
--- /documentation/user_manual/examples/runsydney.py	(revision 2481)
+++ /documentation/user_manual/examples/runsydney.py	(revision 2481)
@@ -0,0 +1,141 @@
+"""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
+from pmesh.mesh_interface import create_mesh_from_regions
+
+# Application specific imports
+import project                           # Definition of file names and polygons
+from pyvolution.smf import slump_tsunami # Function for submarine mudslide 
+
+
+
+#-------------------------------------------------------------------------------   
+# Prepare bathymetric data
+#-------------------------------------------------------------------------------
+
+# filenames
+coarsedemname = project.coarsedemname + '.pts'
+finedemname = project.finedemname + '.pts'
+combineddemname = project.combineddemname + '.pts'
+meshname = project.meshname+'.msh'
+
+
+# combining the coarse and fine data
+combine_rectangular_points_files(finedemname,
+                                 coarsedemname,
+                                 combineddemname)
+
+
+
+#-------------------------------------------------------------------------------   
+# Setup computational domain
+#-------------------------------------------------------------------------------
+
+# Interior regions and mesh resolutions
+interior_res = 5000
+interior_regions = [[project.harbour_polygon_2, interior_res],
+                    [project.botanybay_polygon_2, interior_res]]
+
+
+create_mesh_from_regions(project.diffpolygonall,
+                         boundary_tags= {'bottom': [0], 
+                                         'right1': [1],
+                                         'right0': [2],
+                                         'right2': [3],
+                                         'top': [4],
+                                         'left1': [5],
+                                         'left2': [6],
+                                         'left3': [7]},
+                         maximum_triangle_area=100000,
+                         filename=meshname,           
+                         interior_regions=interior_regions)
+
+
+
+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', tsunami_source)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation',
+                    filename = combineddemname,
+                    use_cache = True,
+                    verbose = True)
+
+
+#-------------------------------------------------------------------------------
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
+                      'right2': Br, 'top': Br, 'left1': Br,
+                      'left2': Br, 'left3': Br} )
+
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 120, finaltime = 18000): 
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')     
+    
+print 'That took %.2f seconds' %(time.time()-t0)
Index: cumentation/user_manual/project.py
===================================================================
--- /documentation/user_manual/project.py	(revision 2480)
+++ 	(revision )
@@ -1,175 +1,0 @@
-"""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
-
-# We'll do this for the final thing 
-home = '.'
-
-
-if sys.platform == 'win32':
-    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
-else:    
-    home = expanduser('~')
-    
-
-basename = 'Sydney'
-
-#Setting the location of scenario data
-scenario_dir_name = 'sydney_tsunami_example_2006'
-# revised 100m data
-coarsename = 'bathyland100' # get from Neil/Ingo (DEM or topo data)
-# revised 25m data
-finename = 'bathy_dem25' # get from Neil/Ingo (DEM or topo data) Wed 25 Jan
-
-
-#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 + 'sydneytopo'
-
-#Georeferencing
-from coordinate_transforms.redfearn import degminsec2decimal_degrees
-
-refzone = 56
-print 'refzone', refzone
-
-dsouth = degminsec2decimal_degrees(-34,05,0)
-dnorth = degminsec2decimal_degrees(-33,33,0)
-dnorth1 = degminsec2decimal_degrees(-33,40,0)
-dnorth2 = degminsec2decimal_degrees(-33,58,30)
-dnorth3 = degminsec2decimal_degrees(-33,46,0)
-dwest = degminsec2decimal_degrees(151,2,20)
-deast1 = degminsec2decimal_degrees(151,20,0)
-deast2 = degminsec2decimal_degrees(151,48,0)
-deast3 = degminsec2decimal_degrees(151,10,0)
-deast4 = degminsec2decimal_degrees(151,9,0)
-
-dp0 = [dsouth, dwest]
-dp1 = [dsouth, deast1]
-dp2 = [dnorth2, deast2]
-dp3 = [dnorth1, deast2]
-dp4 = [dnorth, deast1]
-dp5 = [dnorth, deast4]
-dp6 = [dnorth3, deast3]
-dp7 = [dnorth3, dwest]
-dp8 = [dnorth, dwest]
-    
-diffpolygonall, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7])
-
-qp0 = [343494.0, 6270650.0]
-qp1 = [343337.0, 6270303.0]
-qp2 = [343466.0, 6270228.0]
-qp3 = [343139.0, 6269901.0]
-qp4 = [342472.0, 6268573.0]
-qp5 = [342111.0, 6267115.0]
-qp6 = [342479.0, 6266121.0]
-qp7 = [342860.0, 6266386.0]
-qp8 = [342635.0, 6267438.0]
-qp9 = [343105.0, 6269070.0]
-qp10 = [343548.0, 6269567.0]
-qp11 = [343487.0, 6269928.0]
-qp12 = [343991.0, 6270269.0]
-
-#Interior regions - the Harbour
-harbour_1x = degminsec2decimal_degrees(-33,51,0)
-harbour_1y = degminsec2decimal_degrees(151,2,30)
-harbour_12x = degminsec2decimal_degrees(-33,51,0)
-harbour_12y = degminsec2decimal_degrees(151,5,0)
-harbour_13x = degminsec2decimal_degrees(-33,52,15)
-harbour_13y = degminsec2decimal_degrees(151,5,0)
-harbour_2x = degminsec2decimal_degrees(-33,53,0)
-harbour_2y = degminsec2decimal_degrees(151,17,20)
-harbour_3x = degminsec2decimal_degrees(-33,47,0)
-harbour_3y = degminsec2decimal_degrees(151,20,30)
-harbour_4x = degminsec2decimal_degrees(-33,47,50)
-harbour_4y = degminsec2decimal_degrees(151,8,10)
-harbour_5x = degminsec2decimal_degrees(-33,48,10)
-harbour_5y = degminsec2decimal_degrees(151,8,0)
-harbour_6x = degminsec2decimal_degrees(-33,49,0)
-harbour_6y = degminsec2decimal_degrees(151,2,30)
-harbour_7x = degminsec2decimal_degrees(-33,34,30)
-harbour_7y = degminsec2decimal_degrees(151,20,20)
-harbour_8x = degminsec2decimal_degrees(-33,33,30)
-harbour_8y = degminsec2decimal_degrees(151,17,0)
-harbour_9x = degminsec2decimal_degrees(-33,45,30)
-harbour_9y = degminsec2decimal_degrees(151,17,0)
-harbour_10x = degminsec2decimal_degrees(-33,45,10)
-harbour_10y = degminsec2decimal_degrees(151,11,40)
-harbour_11x = degminsec2decimal_degrees(-33,45,10)
-harbour_11y = degminsec2decimal_degrees(151,11,40)
-harbour_14x = degminsec2decimal_degrees(-33,49,10)
-harbour_14y = degminsec2decimal_degrees(151,11,40)
-harbour_15x = degminsec2decimal_degrees(-33,48,55)
-harbour_15y = degminsec2decimal_degrees(151,2,30)
-
-k02 = [harbour_1x, harbour_1y]
-k12 = [harbour_2x, harbour_2y]
-k22 = [harbour_3x, harbour_3y]
-k32 = [harbour_4x, harbour_4y]
-k42 = [harbour_5x, harbour_5y]
-k52 = [harbour_6x, harbour_6y]
-k62 = [harbour_7x, harbour_7y]
-k72 = [harbour_8x, harbour_8y]
-k82 = [harbour_9x, harbour_9y]
-k92 = [harbour_10x, harbour_10y]
-k102 = [harbour_11x, harbour_11y]
-k112 = [harbour_12x, harbour_12y]
-k122 = [harbour_13x, harbour_13y]
-k132 = [harbour_14x, harbour_14y]
-k142 = [harbour_15x, harbour_15y]
-
-harbour_polygon_2, zone = convert_points_from_latlon_to_utm([k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52]) 
-assert zone == refzone
-
-
-#Interior region - Botany Bay
-bb_1x = degminsec2decimal_degrees(-34,3,0)
-bb_1y = degminsec2decimal_degrees(151,2,30)
-bb_10x = degminsec2decimal_degrees(-34,3,0)
-bb_10y = degminsec2decimal_degrees(151,8,0)
-bb_2x = degminsec2decimal_degrees(-34,3,0)
-bb_2y = degminsec2decimal_degrees(151,14,0)
-bb_3x = degminsec2decimal_degrees(-33,53,30)
-bb_3y = degminsec2decimal_degrees(151,17,20)
-bb_4x = degminsec2decimal_degrees(-33,53,0)
-bb_4y = degminsec2decimal_degrees(151,8,0)
-bb_5x = degminsec2decimal_degrees(-33,57,30)
-bb_5y = degminsec2decimal_degrees(151,8,0)
-bb_6x = degminsec2decimal_degrees(-33,57,30)
-bb_6y = degminsec2decimal_degrees(151,2,30)  
-bb_7x = degminsec2decimal_degrees(-33,53,30)
-bb_7y = degminsec2decimal_degrees(151,12,30)  
-bb_8x = degminsec2decimal_degrees(-33,55,20)
-bb_8y = degminsec2decimal_degrees(151,8,0)  
-bb_9x = degminsec2decimal_degrees(-33,55,20)
-bb_9y = degminsec2decimal_degrees(151,12,30)  
-
-j02 = [bb_1x, bb_1y]
-j12 = [bb_2x, bb_2y]
-j22 = [bb_3x, bb_3y]
-j32 = [bb_4x, bb_4y]
-j42 = [bb_5x, bb_5y]
-j52 = [bb_6x, bb_6y]
-j62 = [bb_7x, bb_7y]
-j72 = [bb_8x, bb_8y]
-j82 = [bb_9x, bb_9y]
-j92 = [bb_10x, bb_10y]
-
-botanybay_polygon_2, zone = convert_points_from_latlon_to_utm([j92, j12, j22, j62, j82, j72, j42]) 
-
-
-
-slump_origin = [385000.0, 6255000.0] #Absolute UTM
Index: cumentation/user_manual/runsydney.py
===================================================================
--- /documentation/user_manual/runsydney.py	(revision 2480)
+++ 	(revision )
@@ -1,141 +1,0 @@
-"""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
-from pmesh.mesh_interface import create_mesh_from_regions
-
-# Application specific imports
-import project                           # Definition of file names and polygons
-from pyvolution.smf import slump_tsunami # Function for submarine mudslide 
-
-
-
-#-------------------------------------------------------------------------------   
-# Prepare bathymetric data
-#-------------------------------------------------------------------------------
-
-# filenames
-coarsedemname = project.coarsedemname + '.pts'
-finedemname = project.finedemname + '.pts'
-combineddemname = project.combineddemname + '.pts'
-meshname = project.meshname+'.msh'
-
-
-# combining the coarse and fine data
-combine_rectangular_points_files(finedemname,
-                                 coarsedemname,
-                                 combineddemname)
-
-
-
-#-------------------------------------------------------------------------------   
-# Setup computational domain
-#-------------------------------------------------------------------------------
-
-# Interior regions and mesh resolutions
-interior_res = 5000
-interior_regions = [[project.harbour_polygon_2, interior_res],
-                    [project.botanybay_polygon_2, interior_res]]
-
-
-create_mesh_from_regions(project.diffpolygonall,
-                         boundary_tags= {'bottom': [0], 
-                                         'right1': [1],
-                                         'right0': [2],
-                                         'right2': [3],
-                                         'top': [4],
-                                         'left1': [5],
-                                         'left2': [6],
-                                         'left3': [7]},
-                         maximum_triangle_area=100000,
-                         filename=meshname,           
-                         interior_regions=interior_regions)
-
-
-
-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', tsunami_source)
-domain.set_quantity('friction', 0.0) 
-domain.set_quantity('elevation',
-                    filename = combineddemname,
-                    use_cache = True,
-                    verbose = True)
-
-
-#-------------------------------------------------------------------------------
-# Setup boundary conditions (all reflective)
-#-------------------------------------------------------------------------------
-
-print 'Available boundary tags', domain.get_boundary_tags()
-
-Br = Reflective_boundary(domain)
-domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
-                      'right2': Br, 'top': Br, 'left1': Br,
-                      'left2': Br, 'left3': Br} )
-
-
-#-------------------------------------------------------------------------------
-# Evolve system through time
-#-------------------------------------------------------------------------------
-
-import time
-t0 = time.time()
-
-for t in domain.evolve(yieldstep = 120, finaltime = 18000): 
-    domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom')     
-    
-print 'That took %.2f seconds' %(time.time()-t0)
