Index: /anuga_work/development/attenuation_near_shore/create_mesh.py
===================================================================
--- /anuga_work/development/attenuation_near_shore/create_mesh.py	(revision 5160)
+++ /anuga_work/development/attenuation_near_shore/create_mesh.py	(revision 5161)
@@ -65,4 +65,27 @@
     print "mesh created"
 
+
+def generate_structured():
+    """
+    """
+    from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+    from anuga.shallow_water import Domain
+    #------------------------------------------------------------------------------
+    # Setup domain
+    #------------------------------------------------------------------------------
+    dx = 10.
+    dy = dx
+    L = 2080.
+    W = 2000.
+
+    # structured mesh
+    points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (-80.0, -1000.0))
+
+    domain = Domain(points, vertices, boundary) 
+
+    print "domain created"
+    return domain
+
+
     
 #-------------------------------------------------------------
Index: /anuga_work/development/attenuation_near_shore/run_beach_structured.py
===================================================================
--- /anuga_work/development/attenuation_near_shore/run_beach_structured.py	(revision 5161)
+++ /anuga_work/development/attenuation_near_shore/run_beach_structured.py	(revision 5161)
@@ -0,0 +1,198 @@
+"""
+
+Script for running a breaking wave simulation of Jon Hinwoods wave tank.
+Note: this is based on the frinction_ua_flume_2006 structure.
+
+
+Duncan Gray, GA - 2007 
+
+
+
+"""
+
+
+#----------------------------------------------------------------------------
+# Import necessary modules
+#----------------------------------------------------------------------------
+
+# Standard modules
+import time
+from time import localtime, strftime
+import sys
+from shutil import copy
+from os import path, sep
+from os.path import dirname  #, basename
+from math import sin, pi
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary,  Time_boundary, \
+                            File_boundary, \
+                            Transmissive_Momentum_Set_Stage_boundary
+from anuga.fit_interpolate.interpolate import interpolate_sww2csv
+from anuga.abstract_2d_finite_volumes.util import copy_code_files, \
+     file_function
+
+from anuga.shallow_water.data_manager import start_screen_catcher
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import File_boundary_time
+
+# Scenario specific imports
+import project                 # Definition of file names and polygons
+import create_mesh
+
+def elevation_function(x,y):
+    from Numeric import zeros, size, Float
+
+    z = zeros(size(x), Float)
+    for i in range(len(x)):
+        if x[i] <= 0:
+            # swash
+            z[i] = -0.05*x[i] 
+            #z[i] = -30
+        elif x[i] <= 1000:
+            # surf
+            z[i] = -0.15*(x[i]**(0.5))
+            #z[i] = -30
+        else:
+            # Nominal shoreface
+            z[i] = -0.01888*(x[i]**(0.8))
+            #z[i] = -30
+            #>>> -0.15*(1000**0.5)
+            #-4.7434164902525691
+
+            #>>> 0.01888*(1000**0.8)
+            #4.7424415826900885
+
+            #>>> -0.01888*10000**0.8
+            #-29.922783473665834
+    return z 
+
+def main(friction=0.01, outputdir_name=None, is_trial_run=False):
+    basename = 'zz' + str(friction)
+    if is_trial_run is True:
+        outputdir_name += '_test'
+        yieldstep = 1
+        finaltime = 100.
+        maximum_triangle_area=3000
+        thinner=True
+    else:
+        yieldstep = 10
+        finaltime = 1600
+        finaltime = 400
+        maximum_triangle_area = 100
+        thinner=True
+        
+        
+    pro_instance = project.Project(['results'],
+                                   outputdir_name=outputdir_name,
+                                   home='.')
+    print "The output dir is", pro_instance.outputdir
+    copy_code_files(pro_instance.outputdir,__file__,
+                    dirname(project.__file__) \
+                    + sep + project.__name__+'.py')
+    copy (pro_instance.codedir + 'run_beach.py',
+          pro_instance.outputdir + 'run_beach.py')
+    copy (pro_instance.codedir + 'create_mesh.py',
+          pro_instance.outputdir + 'create_mesh.py')
+    
+    #mesh_filename = pro_instance.meshdir + basename + '.tsh'
+    mesh_filename = pro_instance.meshdir + basename + '.msh'
+
+    #--------------------------------------------------------------------------
+    # Copy scripts to output directory and capture screen
+    # output to file
+    #--------------------------------------------------------------------------
+
+    # creates copy of code in output dir
+    if is_trial_run is False:
+        pass
+        #start_screen_catcher(pro_instance.outputdir) #, rank, pypar.size())
+
+    print 'USER:    ', pro_instance.user
+    #-------------------------------------------------------------------------
+    # Create the triangular mesh
+    #-------------------------------------------------------------------------
+
+    # this creates the mesh
+    #gate_position = 12.0
+    #create_mesh.generate(mesh_filename,
+    #                     maximum_triangle_area=maximum_triangle_area,
+    #                     thinner=thinner)
+
+    head,tail = path.split(mesh_filename)
+    copy (mesh_filename,
+          pro_instance.outputdir + tail )
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    #domain = Domain(mesh_filename, use_cache = False, verbose = True)
+    domain = create_mesh.generate_structured()
+   
+
+    print 'Number of triangles = ', len(domain)
+    print 'The extent is ', domain.get_extent()
+    print domain.statistics()
+
+    
+    domain.set_name(basename)
+    domain.set_datadir(pro_instance.outputdir)
+    domain.set_default_order(2) # Use second order spatial scheme
+    domain.set_timestepping_method('rk2')
+    domain.beta_w      = 2.0
+    domain.beta_uh     = 2.0
+    domain.beta_vh     = 2.0
+    domain.beta_h      = 2.0
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.set_minimum_storable_height(0.001)
+    domain.set_store_vertices_uniquely(True)  # for writting to sww
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+
+    domain.set_quantity('stage', 0.0)
+    domain.set_quantity('friction', friction) 
+    domain.set_quantity('elevation', elevation_function)
+
+    
+    print 'Available boundary tags', domain.get_boundary_tags()
+
+    # Create boundary function from timeseries provided in file
+    #function = file_function(project.boundary_file, domain, verbose=True)
+    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([10.,0,0])
+    Bwp = Transmissive_Momentum_Set_Stage_boundary(domain,
+                   lambda t: [(1.2*sin(2*t*pi/50)), 0.0 ,0.0])
+    Bwp_velocity = Time_boundary(domain,
+                   lambda t: [(1.2*sin(2*t*pi/10)),
+                              0.24*sin(2*t*pi/10),0.0])
+    Bws = Transmissive_Momentum_Set_Stage_boundary(domain,
+                   lambda t: [1.2*sin(2*t*pi/25),
+                              0.0, 0.0])
+    #domain.set_boundary( {'wall': Br, 'wave': Bwp_velocity} )
+    domain.set_boundary( {'top': Br, 'bottom': Br , 'left' : Br , 'right' : Bws } )
+
+    #-------------------------------------------------------------------------
+    # Evolve system through time
+    #-------------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep, finaltime):    
+        domain.write_time()
+    print 'That took %.2f seconds' %(time.time()-t0)
+    print 'finished'
+
+    #-------------------------------------------------------------------------
+    # Calculate gauge info
+    #-------------------------------------------------------------------------
+
+  
+    return pro_instance
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    main( is_trial_run = False, friction=0.0,
+         outputdir_name='Velocity_wave_2nd_order')
Index: /anuga_work/development/attenuation_near_shore/visualise_height.py
===================================================================
--- /anuga_work/development/attenuation_near_shore/visualise_height.py	(revision 5161)
+++ /anuga_work/development/attenuation_near_shore/visualise_height.py	(revision 5161)
@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+##########
+# Demonstration of the VTK sww Visualiser
+# Jack Kelly
+# September 2006
+##########
+
+# Import the offline visualiser
+from anuga.visualiser import OfflineVisualiser
+from vtk import vtkCubeAxesActor2D
+from Numeric import sqrt
+
+o = OfflineVisualiser("zz0.0.sww")
+
+
+# Specify the height-based-quantities to render.
+# Remember to set dynamic=True for time-varying quantities
+o.render_quantity_height("elevation", zScale=100.0, offset = 5.0, dynamic=False)
+o.render_quantity_height("stage", zScale=100.0, opacity=0.7, dynamic=True)
+
+# Colour the stage:
+o.colour_height_quantity('stage', (lambda q:(q['stage']), -1.5, 1.5))
+#o.colour_height_quantity('stage', (0,0,0.7))
+
+#o.colour_height_quantity('stage', (lambda q:sqrt(((q['xmomentum']/(q['stage']-q['elevation'])) ** 2) +
+#                                             ((q['ymomentum']/(q['stage']-q['elevation'])) ** 2)), 0, 0.5))
+
+# Start the visualiser (in its own thread).
+o.start()
+
+# Wait for the visualiser to terminate before shutting down.
+o.join()
