Index: /anuga_work/production/busselton_2006/build_busselton.py
===================================================================
--- /anuga_work/production/busselton_2006/build_busselton.py	(revision 4132)
+++ /anuga_work/production/busselton_2006/build_busselton.py	(revision 4132)
@@ -0,0 +1,163 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+
+# Application specific imports
+import project   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+start_screen_catcher(project.output_build_time_dir)
+
+print 'time stamp: ',project.gtime
+print 'USER: ', project.user
+
+
+copy_code_files(project.output_build_time_dir,__file__, 
+               dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+# Do for coarse and fine data
+# Fine pts file to be clipped to area of interest
+#-------------------------------------------------------------------------------
+
+# topography directory filenames
+onshore_in_dir_name = project.onshore_in_dir_name
+coast_in_dir_name = project.coast_in_dir_name
+#island_in_dir_name = project.island_in_dir_name
+offshore_in_dir_name = project.offshore_in_dir_name
+
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+#island_dir_name = project.island_dir_name
+offshore_dir_name = project.offshore_dir_name
+
+# creates DEM from asc data
+print "creates DEMs from ascii data"
+convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
+#convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name,
+#        easting_min=project.eastingmin,
+#        easting_max=project.eastingmax,
+#        northing_min=project.northingmin,
+#        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
+
+#creates pts file for island DEM
+#dem2pts(island_dir_name, use_cache=True, verbose=True)
+
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
+#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+
+print'add all geospatial objects'
+G = G1 + G2 + G_off 
+
+print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project.bounding_polygon)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+G.export_points_file(project.combined_dir_name + '.xya')
+#G_clipped.export_points_file(project.combined_dir_name + '.xya')
+
+'''
+print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww
+
+print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
+print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
+
+#import sys; sys.exit()
+
+#if access(project.boundaries_dir,F_OK) == 0:
+#    mkdir (project.boundaries_dir)
+
+from caching import cache
+cache(urs2sww,
+      (boundaries_in_dir_name,
+       project.boundaries_dir_name1), 
+      {'verbose': True,
+       'minlat': project.south_boundary,
+       'maxlat': project.north_boundary,
+       'minlon': project.west_boundary,
+       'maxlon': project.east_boundary,
+#       'minlat': project.south,
+#       'maxlat': project.north,
+#       'minlon': project.west,
+#       'maxlon': project.east,
+       'mint': 0, 'maxt': 40000,
+#       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': project.tide,
+#       'zscale': 1,                 #Enhance tsunami
+       'fail_on_NaN': False},
+       verbose = True,
+       )
+#       dependencies = source_dir + project.boundary_basename + '.sww')
+
+'''
+
+
+
+
+
+
+
Index: /anuga_work/production/busselton_2006/project.py
===================================================================
--- /anuga_work/production/busselton_2006/project.py	(revision 4131)
+++ /anuga_work/production/busselton_2006/project.py	(revision 4132)
@@ -7,5 +7,5 @@
 import sys
 from time import localtime, strftime, gmtime
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
 from anuga.utilities.system_tools import get_user_name
@@ -26,5 +26,5 @@
 build_time = time+'_build'
 run_time = time+'_run'
-print 'gtime: ', gtime
+#print 'gtime: ', gtime
 
 tide = 0.6
@@ -105,8 +105,6 @@
 # bounding polygon for study area
 poly_all = read_polygon(polygons_dir+'poly_all.csv')
-res_poly_all = 100000
+res_poly_all = 1000000
 
-
-'''
 ###################################################################
 # Clipping regions for export to asc and regions for clipping data
@@ -123,27 +121,13 @@
 ###############################
 
-# broome digitized polygons
-poly_busselton1 = read_polygon(polygons_dir+'buss_Local_Polygon_update.csv')
-res_busselton1 = 500000
-poly_busselton2 = read_polygon(polygons_dir+'buss_Close2_update.csv')
-res_busselton2 = 500000
-poly_busselton3 = read_polygon(polygons_dir+'buss_Coast_update.csv')
-res_busselton3 = 500000
+#digitized polygons
+poly_busselton1 = read_polygon(polygons_dir+'neg20_pos10_polygon.csv')
+res_busselton1 = 100000
+poly_busselton2 = read_polygon(polygons_dir+'neg5_pos5_poly_.csv')
+res_busselton2 = 10000
 
 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
 
-interior_regions = [[poly_busselton1,res_busselton1],[poly_busselton2,res_busselton2],
-                    [poly_busselton3,res_busselton3]]
-#trigs_all = polygon_area(poly_all)/res_poly_all
-#trigs_bus1 = polygon_area(poly_busselton1)/res_busselton1
-#trigs_bus2 = polygon_area(poly_busselton2)/res_busselton2
-#trigs_bus3 = polygon_area(poly_busselton3)/res_busselton3
-#trigs_min = trigs_bound + trigs_pos + trigs_cbd + trigs_penguin
-
-#print 'Area of bounding poly', trigs_all
-#print 'Area of busselton1', trigs_bus1
-#print 'Area of busselton2', trigs_bus2
-#print 'Area of busselton3', trigs_bus3
-#print 'min number triangles', trigs_min
+interior_regions = [[poly_busselton1,res_busselton1],[poly_busselton2,res_busselton2]]
 
 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
@@ -151,17 +135,3 @@
 print 'min number triangles', trigs_min
 
-def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
-    from anuga.utilities.polygon import polygon_area
-    
-    # TO DO check if any of the regions fall inside one another
-    no_triangles = 0.0
-    area = polygon_area(bounding_poly)
-    for i,j in interior_regions:
-        this_area = polygon_area(i)
-        no_triangles += this_area/j
-        area -= this_area
-        print j, this_area/1000000., area/1000000.
-    no_triangles += area/remainder_res
-    return int(no_triangles/0.7)
-'''
 
Index: /anuga_work/production/busselton_2006/run_busselton.py
===================================================================
--- /anuga_work/production/busselton_2006/run_busselton.py	(revision 4132)
+++ /anuga_work/production/busselton_2006/run_busselton.py	(revision 4132)
@@ -0,0 +1,220 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+from Numeric import allclose
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+from anuga.utilities.polygon import plot_polygons, polygon_area
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+start_screen_catcher(project.output_run_time_dir, myid, numprocs)
+
+print 'time stamp: ',project.gtime
+print 'USER: ', project.user
+print 'min triangles', project.trigs_min
+
+# filenames
+boundaries_name = project.scenario
+meshes_dir_name = project.meshes_dir_name+'.msh'
+#boundaries_time_dir_name = project.boundaries_time_dir_name
+boundaries_dir_name = project.boundaries_dir_name
+
+tide = project.tide
+
+# creates copy of code in output dir
+if myid == 0:
+    copy_code_files(project.output_run_time_dir,__file__, 
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
+barrier()
+
+#--------------------------------------------------------------------------
+# 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
+#--------------------------------------------------------------------------
+
+if myid == 0:
+
+    print 'start create mesh from regions'
+#    interior_regions = [[project.poly_busselton1, project.res_bus1],
+#                    [project.poly_busselton2, project.res_bus2],
+#                    [project.poly_busselton3, project.res_bus3]]    
+
+#    from anuga.utilities.polygon import plot_polygons
+#    figname = project.output_run_time_dir + 'poly_pic'
+#    plot_polygons([project.poly_coast,project.poly_pipeline,project.poly_facility,
+#               project.bounding_polygon],
+#               figname,
+#               verbose = True)
+#    import sys; sys.exit()
+
+    print 'start create mesh from regions'
+    create_mesh_from_regions(project.poly_all,
+                         boundary_tags={'back': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 14],
+                                        'side': [10,13], 'ocean': [11, 12]},
+                         maximum_triangle_area=project.res_poly_all,
+                         interior_regions=project.interior_regions,
+#                         filename=meshes_time_dir_name,
+                         filename=meshes_dir_name,
+                         use_cache=True,
+                         verbose=True)
+
+# to sync all processors are ready 
+barrier()
+
+#-------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------
+print 'Setup computational domain'
+#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
+domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
+print domain.statistics()
+'''
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww, ferret2sww
+
+# put above distribute
+print 'boundary file is: ',boundaries_dir_name
+from caching import cache
+if myid == 0:
+    cache(ferret2sww,
+          (boundaries_in_dir_name,
+    #       boundaries_time_dir_name), 
+           boundaries_dir_name), 
+          {'verbose': True,
+           'minlat': project.south_boundary,
+           'maxlat': project.north_boundary,
+           'minlon': project.west_boundary,
+           'maxlon': project.east_boundary,
+#           'minlat': project.south,
+#           'maxlat': project.north,
+#           'minlon': project.west,
+#           'maxlon': project.east,
+           'mint': 0, 'maxt': 35100,
+           'origin': domain.geo_reference.get_origin(),
+           'mean_stage': project.tide,
+#           'zscale': 1,                 #Enhance tsunami
+           'fail_on_NaN': False},
+           verbose = True,
+           )
+barrier()
+'''
+
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+if myid == 0:
+
+    print 'Setup initial conditions'
+
+    domain.set_quantity('stage', tide)
+    domain.set_quantity('friction', 0.01) 
+    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
+    print 'Start Set quantity'
+
+    domain.set_quantity('elevation', 
+#                    filename = project.combined_smaller_dir_name+'.xya',
+                    filename = project.combined_dir_name+'.pts',
+                    use_cache = True,
+                    verbose = True,
+                    alpha = 0.1)
+    print 'Finished Set quantity'
+barrier()
+
+#------------------------------------------------------
+# Distribute domain to implement parallelism !!!
+#------------------------------------------------------
+
+if numprocs > 1:
+    domain=distribute(domain)
+
+#------------------------------------------------------
+# Set domain parameters
+#------------------------------------------------------
+
+domain.set_name(project.scenario_name)
+domain.set_datadir(project.output_run_time_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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+
+#-------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------
+print 'Available boundary tags', domain.get_boundary_tags()
+
+print 'Reading Boundary file'
+print 'domain id', id(domain)
+#Bf = File_boundary(boundaries_dir_name + '.sww',
+#                  domain, time_thinning=5, use_cache=True, verbose=True)
+
+print 'finished reading boundary file'
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+Bo = Dirichlet_boundary([tide+5.0,0,0])
+
+print'set_boundary'
+domain.set_boundary({'back': Br,
+                     'side': Bd,
+                     'ocean': Bd}) 
+print'finish set boundary'
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#----------------------------------------------------------------------------
+
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 60, finaltime = 10000):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'ocean')
+    if allclose(t, 120):
+        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
+
+    if allclose(t, 1020):
+        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
+
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
