Index: /anuga_work/development/saman_flume/create_mesh.py
===================================================================
--- /anuga_work/development/saman_flume/create_mesh.py	(revision 4444)
+++ /anuga_work/development/saman_flume/create_mesh.py	(revision 4444)
@@ -0,0 +1,143 @@
+"""Create mesh for University of Queensland dam break flume
+"""
+
+from math import cos, radians
+from anuga.pmesh.mesh import *
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+
+
+#Basic geometry
+xslope = 2.5
+
+def generate(mesh_filename, beach_angle, toe_to_nails_distance,
+             is_coarse = True):
+    """
+    Generate mesh for Saman's flume flume.
+    The gate position is the distance from the right wall of the wave tank
+    to the gate. 
+    
+    """
+    #Basic geometry
+    
+    xright  = 5.25 
+    ybottom = 0
+    ytop    = 0.30
+    xdam = 0.0
+    xleft = -0.5
+    #xslope = self.xslope
+    
+    nails_length_x = 0.15 
+
+    xnailsleft = xslope + toe_to_nails_distance*cos(radians(beach_angle))
+    xnailsright = xslope + (toe_to_nails_distance+ nails_length_x) \
+                  *cos(radians(beach_angle))
+
+    
+
+    #Outline 
+    point_sw = [xleft, ybottom]
+    point_se = [xright, ybottom]
+    point_nw = [xleft, ytop]    
+    point_ne = [xright, ytop]
+
+    # Dam seperation (left)
+    point_dam_top = [xdam, ytop]
+    point_dam_bottom = [xdam, ybottom]    
+
+    # slope seperation (middle)
+    point_slope_top = [xslope, ytop]
+    point_slope_bottom = [xslope, ybottom]
+    
+    # nail seperation (after slope, left)
+    point_nails_left_top = [xnailsleft, ytop]
+    point_nails_left_bottom = [xnailsleft, ybottom]
+    
+    # nail seperation (after slope, right)
+    point_nails_right_top = [xnailsright, ytop]
+    point_nails_right_bottom = [xnailsright, ybottom]
+     
+
+    m = Mesh()
+
+    #Boundary
+    points = [point_sw,   #se
+              point_nw,
+              point_dam_top,
+              point_slope_top,
+              point_nails_left_top,
+              point_nails_right_top,
+              point_ne, 
+              point_se,
+              point_nails_right_bottom,
+              point_nails_left_bottom,
+              point_slope_bottom,
+              point_dam_bottom
+              ]
+    
+    segments = [[0,1], [1,2], [2,3],
+                [3,4 ],[4,5], [5,6],[6,7],[7,8],
+                [8,9], [9,10],[10,11],[11,0],#The outer border
+                [2,11],    #dam Separator
+                [3,10],        # slope separator
+                [4,9],        # nails left separator
+                [5,8]]        # nails right separator
+    
+    segment_tags = {'wall':[0,1,2,3,5,6,7,8,9,10,11],'edge':[4]} 
+        
+    m.add_points_and_segments(points, segments, segment_tags)
+
+    
+    if is_coarse == True:
+        #nails_max_area = 0.001
+        general_max_area = 0.01
+        segment_count = 8
+        nails_x = 2
+        nails_y = 5
+        nail_spacing = 0.060  # distance between nail centers
+        nail_radius = 0.01
+    else:
+        #nails_max_area = 0.000001
+        general_max_area = 0.0001
+        segment_count = 20
+        nails_x = 8
+        nails_y = 19
+        nail_spacing = 0.015  # distance between nail centers
+        nail_radius = 0.0025
+
+    #Add the nails
+    for x_i in range(1,nails_x + 1):
+        for y_i in range(1,nails_y + 1):
+            center_point = (xnailsleft+(x_i*nail_spacing), y_i*nail_spacing)
+            m.add_circle(center_point, radius = nail_radius,
+                         segment_count = segment_count,
+                         tag = 'wall',
+                         hole = True)
+             
+    
+    # this is the location of the reservoir region.
+    dam = m.add_region(-0.0000001,(ytop - ybottom)/2)
+    dam.setTag("dam")
+    
+    # this is the location of the slope region, before the nails.
+    slope = m.add_region(xslope + 0.0000001,(ytop - ybottom)/2)
+    slope.setTag("slope")
+    
+    # this is the location of the slope region, after the nails.
+    slope = m.add_region(xnailsright + 0.0000001,(ytop - ybottom)/2)
+    slope.setTag("slope")
+
+    # this is the location of the nails region, after the nails.
+    nails = m.add_region(xnailsleft + 0.0000001,(ytop - ybottom)/2)
+    nails.setTag("nails")
+    #nails.setMaxArea(nails_max_area)
+    
+    
+    m.generate_mesh(maximum_triangle_area=general_max_area)
+
+    m.export_mesh_file(mesh_filename)
+    print "mesh created"
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    generate("aa.tsh", beach_angle = 15, toe_to_nails_distance = 1.0,
+             is_coarse = False)
Index: /anuga_work/development/saman_flume/project.py
===================================================================
--- /anuga_work/development/saman_flume/project.py	(revision 4444)
+++ /anuga_work/development/saman_flume/project.py	(revision 4444)
@@ -0,0 +1,50 @@
+"""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
+
+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/saman_flume/run_dam.py
===================================================================
--- /anuga_work/development/saman_flume/run_dam.py	(revision 4444)
+++ /anuga_work/development/saman_flume/run_dam.py	(revision 4444)
@@ -0,0 +1,178 @@
+"""Script for running a dam break simulation of UQ's dam break tank.
+
+
+Ole Nielsen and Duncan Gray, GA - 2006 
+
+
+Issues
+If running this is hand-set-up parallel, the python files are overwritten.
+
+"""
+
+
+#----------------------------------------------------------------------------
+# 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 tan, radians
+
+# Related major packages
+from anuga.shallow_water import Domain, Reflective_boundary, \
+                            Dirichlet_boundary, Time_boundary, File_boundary
+from anuga.abstract_2d_finite_volumes.region import Set_region 
+from anuga.fit_interpolate.interpolate import interpolate_sww2csv
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
+     copy_code_files
+
+
+# Scenario specific imports
+import project                 # Definition of file names and polygons
+import create_mesh
+
+
+
+def main(friciton,beach_angle, toe_to_nails_distance,
+         outputdir_name=None, is_trial_run=False):
+    """
+    beach angle is the angle in degrees of the sloped section
+    toe_to_nails_distance is the distance from the beginning if the sloped
+    section to the beginning of the nails section
+    """
+    basename = 'zz' + str(friction)
+    if is_trial_run is True:
+        outputdir_name += '_test'
+        yieldstep = 0.1
+        finaltime = 4
+    else:
+        yieldstep = 0.01
+        finaltime = 31
+    pro_instance = project.Project(['data','flumes','saman_2007'],
+                                       outputdir_name=outputdir_name)
+    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
+    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' + str(friction)+ '.py')
+    copy (pro_instance.codedir + 'create_mesh.py',
+          pro_instance.outputdir + 'create_mesh.py')
+    copy (pro_instance.codedir + 'project.py',
+          pro_instance.outputdir + 'project.py')
+    if is_trial_run is False:
+        start_screen_catcher(pro_instance.outputdir, int(friction*100))
+
+    print 'USER:    ', pro_instance.user
+    #-------------------------------------------------------------------------
+    # Create the triangular mesh
+    #-------------------------------------------------------------------------
+    
+    #beach_angle = 15.0
+    #toe_to_nails_distance = .5
+    create_mesh.generate(mesh_filename,
+                         beach_angle,
+                         toe_to_nails_distance,
+                         is_coarse=is_trial_run) # this creates the mesh
+
+    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.01)
+    #domain.set_store_vertices_uniquely(True)  # for writting to sww
+
+
+    def elevation_function(x,y):
+        from Numeric import zeros, size, Float
+        slope = create_mesh.xslope 
+    
+        z = zeros(size(x), Float)
+        for i in range(len(x)):
+            if x[i] < slope:
+                z[i] = 0.0
+            else:
+                z[i] = (x[i]-slope)*tan(radians(beach_angle))
+        return z
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+        
+    domain.set_quantity('stage', elevation_function)
+    domain.set_quantity('friction', friciton) 
+    domain.set_quantity('elevation', elevation_function)
+
+    print 'Available boundary tags', domain.get_boundary_tags()
+    domain.set_region('dam','stage',0.30,
+                                 location = 'unique vertices') 
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
+    domain.set_boundary( {'wall': Br, 'edge': 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'
+
+
+    #-------------------------------------------------------------------------
+    # Calculate gauge info
+    #-------------------------------------------------------------------------
+    if not is_trial_run:
+        points = [[gate_position - 0.65,0.2],
+                  [gate_position - 0.55,0.2],
+                  [gate_position - 0.45,0.2],
+                  [gate_position - 0.35,0.2],
+                  [gate_position - 0.25,0.2]
+                  ]
+        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__":
+    beach_angle = 3.0
+    toe_to_nails_distance = 0.1
+    for friction in [0.0,0.01]:
+        main(friction,beach_angle, toe_to_nails_distance,
+             is_trial_run = True, outputdir_name='testing')
+    
