Index: /anuga_work/development/near_shore_PMD/create_mesh.py
===================================================================
--- /anuga_work/development/near_shore_PMD/create_mesh.py	(revision 5093)
+++ /anuga_work/development/near_shore_PMD/create_mesh.py	(revision 5093)
@@ -0,0 +1,59 @@
+"""Create mesh for near shore PMD
+"""
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+
+xslope = 5.0 # Distance between the boundary ande the start of the slope
+
+
+def generate(mesh_filename, maximum_triangle_area=300000):
+    """
+    """
+    # Basic geometry and bounding polygon
+    xleft   = -100  # Beach
+    xright  = 10000  # wave generated here
+    ytop    = 15000
+    ybottom = -15000
+
+    point_sw = [xleft, ybottom]
+    point_se = [xright, ybottom]
+    point_nw = [xleft, ytop]    
+    point_ne = [xright, ytop]
+
+    bounding_polygon = [point_se,
+                        point_ne,
+                        point_nw,
+                        point_sw]
+
+    # local refinement
+    xleft   = -80
+    xright  = 500
+    ytop    = 1000
+    ybottom = -1000
+
+    point_sw = [xleft, ybottom]
+    point_se = [xright, ybottom]
+    point_nw = [xleft, ytop]    
+    point_ne = [xright, ytop]
+    
+    polygon_2 = [point_se,
+                        point_ne,
+                        point_nw,
+                        point_sw]
+
+    interior_regions = [[polygon_2, maximum_triangle_area/10.]]
+    
+    m = create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags={'wall': [1, 2, 3],
+                                                'wave': [0]},     
+                                 maximum_triangle_area=maximum_triangle_area,
+                                 interior_regions=interior_regions,
+                                 filename=mesh_filename,
+                                 verbose=True)
+    print "mesh created"
+
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    generate("aa.tsh")
Index: /anuga_work/development/near_shore_PMD/project.py
===================================================================
--- /anuga_work/development/near_shore_PMD/project.py	(revision 5093)
+++ /anuga_work/development/near_shore_PMD/project.py	(revision 5093)
@@ -0,0 +1,52 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+Also includes origin for slump scenario.
+"""
+
+from os import sep, environ, getenv, getcwd
+import sys
+from time import localtime, strftime
+from os import mkdir, access, F_OK
+from anuga.abstract_2d_finite_volumes.util import add_directories, \
+     copy_code_files
+
+from anuga.utilities.system_tools import get_user_name
+
+boundary_file = 'boundary.tsm'
+
+class Project:    
+    def __init__(self,
+                 trunk,
+                 outputdir_name = None,
+                 home = None):
+        
+        self.user = get_user_name()
+        if home is None:
+            home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+        self.home = home
+        # Create the structure of where the output directories will go
+        scenariodir=add_directories(home, trunk)
+
+        #Derive subdirectories and filenames
+        if outputdir_name is None:
+            #gets time for dir
+            outputdir_name = strftime('%Y%m%d_%H%M%S',localtime()) 
+        general_outputdir = scenariodir+sep+'output'+sep
+        self.outputdir = general_outputdir+outputdir_name+sep
+        self.meshdir = scenariodir+sep+'meshes'+sep
+
+        # creates copy of output dir structure, if it doesn't exist
+        if not access(self.meshdir,F_OK):
+            mkdir (self.meshdir)       
+        if not access(general_outputdir,F_OK):
+            mkdir (general_outputdir)
+        if not access(self.outputdir,F_OK):
+            mkdir (self.outputdir)
+
+        self.codedir = getcwd()+sep
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    p = Project(['eagle'], 'lego','.')
+    print p.outputdir
+    print p.meshdir
Index: /anuga_work/development/near_shore_PMD/run_beach.py
===================================================================
--- /anuga_work/development/near_shore_PMD/run_beach.py	(revision 5093)
+++ /anuga_work/development/near_shore_PMD/run_beach.py	(revision 5093)
@@ -0,0 +1,183 @@
+"""
+
+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
+
+
+# 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 start_screen_catcher, \
+     copy_code_files, file_function
+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] 
+        elif x[i] <= 1000:
+            # surf
+            z[i] = -0.15*(x[i]**(0.5))
+        else:
+            # Nominal shoreface
+            z[i] = -0.01888*(x[i]**(0.8))
+            #>>> -0.15*(1000**0.5)
+            #-4.7434164902525691
+
+            #>>> 0.01888*(1000**0.8)
+            #4.7424415826900885
+    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 = 0.1
+        finaltime = 15.
+        maximum_triangle_area=300000
+    else:
+        yieldstep = 0.02
+        finaltime = 15.1
+        maximum_triangle_area=30000
+        
+        
+    pro_instance = project.Project(['data','near_shore_PMD'],
+                                   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_dam.py',
+          pro_instance.outputdir + 'run_dam.py')
+    copy (pro_instance.codedir + 'create_mesh.py',
+          pro_instance.outputdir + 'create_mesh.py')
+    
+    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:
+        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)
+
+    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)
+   
+
+    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_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([1000.,0,0]) 
+    domain.set_boundary( {'wall': Br, 'wave': Bd} )
+
+    #-------------------------------------------------------------------------
+    # 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'
+
+    points = [[2.8,0.5],  #-1.8m from SWL
+              [5.1,0.5],  #0.5m from SWL
+              [6.6,0.5],  #2m from SWL
+              [6.95,0.5], #2.35m from SWL
+              [7.6,0.5],  #3m from SWL
+              [8.2,0.5],  #3.5m from SWL
+              [9.2,0.5]  #4.5m from SWL
+              ]
+
+
+    #-------------------------------------------------------------------------
+    # Calculate gauge info
+    #-------------------------------------------------------------------------
+
+    if False:
+        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
+                            points,
+                            pro_instance.outputdir + "depth_manning_"+str(friction)+".csv",
+                            pro_instance.outputdir + "velocity_x.csv",
+                            pro_instance.outputdir + "velocity_y.csv")
+  
+    return pro_instance
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    main( is_trial_run = True,
+         outputdir_name='Hinwood_low_stage_low_velocity_draft')
