Index: /anuga_work/development/landslide_4/benchmark_problem4.py
===================================================================
--- /anuga_work/development/landslide_4/benchmark_problem4.py	(revision 5567)
+++ /anuga_work/development/landslide_4/benchmark_problem4.py	(revision 5567)
@@ -0,0 +1,124 @@
+"""Script for running the 3rd Long Wave Workshop Benchmark Problem 4.
+
+Bridgette Lewis and Jane Sexton, July 2008
+"""
+
+#-------------------------------------------------------------------------------
+# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+from shutil import copy
+from os.path import dirname, basename
+from os import mkdir, access, F_OK, sep
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+
+#----------------------------------------------------------------------------
+# Create the triangular mesh based on the problem statement
+#-------------------------------------------------------------------------------
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+remainder_res = 250000.
+#local_res = 50000.
+#coast_res = 500.
+#interior_regions = [[project_slide.poly_syd1, local_res],
+#                    [project_slide.poly_coast, coast_res]]
+length = 10.
+width = 10.
+depth = 10.
+bounding_polygon = [[0,width],[0,0],[length,0],[length,width]]
+meshname = 'slide.msh'
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          bounding_polygon,
+           {'boundary_tags': {'top': [0], 'left': [1], 
+                              'bottom': [2], 'right': [3]},
+           'maximum_triangle_area': remainder_res,
+           'filename': meshname},
+           #'interior_regions': interior_regions},
+          verbose = True, evaluate=False)
+print 'created mesh'
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+domain = Domain(meshname, use_cache = True, verbose = True)
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name('slide')
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+from math import tan
+def topography(x,y):
+    z = -depth + tan(0.5)*x
+    return z
+domain.set_quantity('stage', 0.0)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation', topography,
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1
+                    )
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+from smf import slide_tsunami
+
+tsunami_source = slide_tsunami(length=2.,
+                               width=1.,
+                               depth=1.,
+                               slope=45,
+                               thickness=0.1, 
+                               x0=length/2, 
+                               y0=width/2, 
+                               alpha=0.,
+					 #dx=0.001,
+                               domain=domain)
+                               #verbose=True)
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions
+#-------------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+
+domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd} )
+
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+import time
+from Numeric import allclose
+from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 1, finaltime = 10): 
+    domain.write_time()
+    stagestep = domain.get_quantity('stage') 
+
+    if allclose(t, 1):
+        slide = Quantity(domain)
+        slide.set_values(tsunami_source)
+        domain.set_quantity('stage', slide + stagestep)
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
+print 'finished'
