Index: anuga_work/production/patong/run_patong.py
===================================================================
--- anuga_work/production/patong/run_patong.py	(revision 6175)
+++ anuga_work/production/patong/run_patong.py	(revision 6178)
@@ -23,9 +23,6 @@
 from os import sep
 import os
-from os.path import dirname, basename
-from os import mkdir, access, F_OK
-from shutil import copy
+from os.path import dirname
 import time
-import sys
 
 # Related major packages
@@ -35,7 +32,6 @@
 from anuga.shallow_water import Reflective_boundary
 from anuga.shallow_water import Field_boundary
-from Numeric import allclose
 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water import create_domain_from_regions
 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
 from anuga.caching import myhash
@@ -51,201 +47,140 @@
 myid = 0
 
-def run_model(**kwargs):
-    
-    #------------------------------------------------------------------------------
-    # Copy scripts to time stamped output directory and capture screen
-    # output to file
-    #------------------------------------------------------------------------------
-
-    #copy script must be before screen_catcher
-
-    print 'output_dir',kwargs['output_dir']
-    
-    copy_code_files(kwargs['output_dir'],__file__, 
-             dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-    store_parameters(**kwargs)
-
-    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
 
     
-    #-----------------------------------------------------------------------
-    # Domain definitions
-    #-----------------------------------------------------------------------
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
 
-    # Read in boundary from ordered sts file
-    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
+#copy script must be before screen_catcher
 
-    # Reading the landward defined points, this incorporates the original clipping
-    # polygon minus the 100m contour
-    landward_bounding_polygon = read_polygon(project.landward_dir)
+output_dir = project.output_run_time_dir
+print 'output_dir', output_dir
 
-    # Combine sts polyline with landward points
-    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
-    
-    # counting segments
-    N = len(urs_bounding_polygon)-1
+#copy_code_files(output_dir,__file__, 
+#         dirname(project.__file__)+sep+ project.__name__+'.py' )
 
-    # boundary tags refer to project.landward 4 points equals 5 segments start at N
-    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
-
-    #--------------------------------------------------------------------------
-    # Create the triangular mesh based on overall clipping polygon with a tagged
-    # boundary and interior regions defined in project.py along with
-    # resolutions (maximal area of per triangle) for each polygon
-    #--------------------------------------------------------------------------
-
-    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
-    # causes problems with the ability to cache set quantity which takes alot of times
-        
-    print 'start create mesh from regions'
-
-    create_mesh_from_regions(bounding_polygon,
-                             boundary_tags=boundary_tags,
-                             maximum_triangle_area=project.res_poly_all,
-                             interior_regions=project.interior_regions,
-                             filename=project.meshes_dir_name,
-                             use_cache=True,
-                             verbose=True)
-   
-    #-------------------------------------------------------------------------
-    # Setup computational domain
-    #-------------------------------------------------------------------------
-    print 'Setup computational domain'
-
-    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
-    print 'memory usage before del domain',mem_usage()
-       
-    print domain.statistics()
-    print 'triangles',len(domain)
-    
-    kwargs['act_num_trigs']=len(domain)
+#start_screen_catcher(output_dir, myid, numprocs)
 
 
-    #-------------------------------------------------------------------------
-    # Setup initial conditions
-    #-------------------------------------------------------------------------
-    print 'Setup initial conditions'
+#-----------------------------------------------------------------------
+# Domain definitions
+#-----------------------------------------------------------------------
 
-    # sets the initial stage in the offcoast region only
-    IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
-                             geo_reference = domain.geo_reference)
-    domain.set_quantity('stage', IC)
-    #domain.set_quantity('stage',kwargs['tide'] )
-    domain.set_quantity('friction', kwargs['friction']) 
-    
-    print 'Start Set quantity',kwargs['elevation_file']
+# Read in boundary from ordered sts file
+urs_boundary_name = os.path.join(project.boundaries_dir,
+                                 project.scenario_name)
+urs_bounding_polygon = create_sts_boundary(urs_boundary_name)
 
-    domain.set_quantity('elevation', 
-                        filename = kwargs['elevation_file'],
-                        use_cache = False,
-                        verbose = True,
-                        alpha = kwargs['alpha'])
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour
+landward_bounding_polygon = read_polygon(project.landward_dir)
 
-    # Add buildings from file
-    building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
+# Combine sts polyline with landward points
+bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
 
-    L = []
-    for key in building_polygons:
-        poly = building_polygons[key]
-        elev = building_heights[key]
-        L.append((poly, elev))
-    domain.add_quantity('elevation', Polygon_function(L, default=0.0))
-    
-    print 'Finished Set quantity'
+# Counting segments
+N = len(urs_bounding_polygon)-1
 
-##   #------------------------------------------------------
-##    # Distribute domain to implement parallelism !!!
-##    #------------------------------------------------------
-##
-##    if numprocs > 1:
-##        domain=distribute(domain)
+# boundary tags refer to project.landward 4 points equals 5 segments start at N
+boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
 
-    #------------------------------------------------------
-    # Set domain parameters
-    #------------------------------------------------------
-    print 'domain id', id(domain)
-    domain.set_name(kwargs['scenario_name'])
-    domain.set_datadir(kwargs['output_dir'])
-    domain.set_default_order(2)                 # Apply second order scheme
-    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
-    domain.set_store_vertices_uniquely(False)
-    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.tight_slope_limiters = 1
-    print 'domain id', id(domain)
+#--------------------------------------------------------------------------
+# Create the computational domain based on overall clipping polygon with
+# a tagged boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#--------------------------------------------------------------------------
 
-    #-------------------------------------------------------------------------
-    # Setup boundary conditions 
-    #-------------------------------------------------------------------------
-    print 'Available boundary tags', domain.get_boundary_tags()
-    print 'domain id', id(domain)
-   
-    boundary_urs_out=project.boundaries_dir + sep + project.scenario_name
+domain = create_domain_from_regions(bounding_polygon,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.res_poly_all,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes_dir_name,
+                                    use_cache=True,
+                                    verbose=True)
 
-    Br = Reflective_boundary(domain)
-    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    
-    print 'Available boundary tags', domain.get_boundary_tags()
-    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
-                   domain, mean_stage= project.tide,
-                   time_thinning=1,
-                   time_limit=project.finaltime,
-                   default_boundary=Bd,
-                   use_cache=True,
-                   verbose=True,
-                   boundary_polygon=bounding_polygon)
+print 'memory usage after creation of domain', mem_usage()
+print domain.statistics()
 
-    domain.set_boundary({'back': Br,
-                         'side': Bd,
-                         'ocean': Bf}) 
 
-    kwargs['input_start_time']=domain.starttime
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+print 'Setup initial conditions'
 
-    print'finish set boundary'
+# sets the initial stage in the offcoast region only
+IC = Polygon_function([(project.poly_mainland, 0)],
+                      default=project.tide,
+                      geo_reference=domain.geo_reference)
+domain.set_quantity('stage', IC)
+domain.set_quantity('friction', project.friction) 
+domain.set_quantity('elevation', 
+                    filename=project.combined_dir_name+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project.alpha)
 
-    #----------------------------------------------------------------------------
-    # Evolve system through time
-    #--------------------------------------------------------------------
-    t0 = time.time()
+# Add buildings from file
+print 'Creating building polygons'    
+building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
 
-    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
-                       ,skip_initial_step = False): 
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')
+L = []
+for key in building_polygons:
+    poly = building_polygons[key]
+    elev = building_heights[key]
+    L.append((poly, elev))
 
-    # these outputs should be checked with the resultant inundation map
-    x, y = domain.get_maximum_inundation_location()
-    q = domain.get_maximum_inundation_elevation()
-    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+buildings = Polygon_function(L, default=0.0)
+domain.add_quantity('elevation', buildings)
 
-    print 'Simulation took %.2f seconds' %(time.time()-t0)
 
-    #kwargs 'completed' must be added to write the final parameters to file
-    kwargs['completed']=str(time.time()-t0)
-     
-    store_parameters(**kwargs)
-     
-    
-    
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    
-    kwargs={}
-    kwargs['file_name']=project.dir_comment
-    kwargs['finaltime']=project.finaltime
-    kwargs['output_dir']=project.output_run_time_dir
-    kwargs['elevation_file']=project.combined_dir_name+'.pts'
-    kwargs['scenario_name']=project.scenario_name
-    kwargs['tide']=project.tide
-    kwargs['alpha'] = project.alpha
-    kwargs['friction']=project.friction
-    #kwargs['num_cpu']=numprocs
-    #kwargs['host']=project.host
-    #kwargs['starttime']=project.starttime
-    #kwargs['yieldstep']=project.yieldstep
-    #kwargs['user']=project.user
-    #kwargs['time_thinning'] = project.time_thinning
-      
-    run_model(**kwargs)
+#------------------------------------------------------
+# Distribute domain to implement parallelism !!!
+#------------------------------------------------------
+
+#  if numprocs > 1:
+#      domain=distribute(domain)
+
+#------------------------------------------------------
+# Set domain parameters
+#------------------------------------------------------
+domain.set_name(project.scenario_name)
+domain.set_datadir(output_dir)
+domain.set_default_order(2)                 # Apply second order scheme
+domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
+
+#-------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------
+print 'Set boundary conditions'
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([project.tide,0,0])
+Bf = Field_boundary(urs_boundary_name+'.sts', 
+                    domain,
+                    mean_stage= project.tide,
+                    time_thinning=project.time_thinning,
+                    default_boundary=Bd,
+                    use_cache=True,
+                    verbose=True,
+                    boundary_polygon=bounding_polygon)
+
+domain.set_boundary({'back': Br,
+                     'side': Bd,
+                     'ocean': Bf}) 
+
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#--------------------------------------------------------------------
+t0 = time.time()
+
+for t in domain.evolve(yieldstep=project.yieldstep,
+                       finaltime=project.finaltime)
+    print domain.timestepping_statistics()
+    print domain.boundary_statistics(tags='ocean')
+
+print 'Simulation took %.2f seconds' %(time.time()-t0)
       
    
