Index: anuga_work/production/broome_2006/build_broome.py
===================================================================
--- anuga_work/production/broome_2006/build_broome.py	(revision 4153)
+++ anuga_work/production/broome_2006/build_broome.py	(revision 4153)
@@ -0,0 +1,146 @@
+"""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
+#------------------------------------------------------------------------------
+
+copy_code_files(project.output_build_time_dir,__file__, 
+               dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+start_screen_catcher(project.output_build_time_dir)
+
+print 'USER:    ', project.user
+
+#-------------------------------------------------------------------------------
+# 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
+#-------------------------------------------------------------------------------
+print"project.bounding_polygon",project.bounding_polygon
+print"project.combined_dir_name",project.combined_dir_name
+
+# 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
+island_in_dir_name1 = project.island_in_dir_name1
+island_in_dir_name2 = project.island_in_dir_name2
+island_in_dir_name3 = project.island_in_dir_name3
+offshore_in_dir_name = project.offshore_in_dir_name
+offshore1_in_dir_name = project.offshore1_in_dir_name
+
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+island_dir_name = project.island_dir_name
+island_dir_name1 = project.island_dir_name1
+island_dir_name2 = project.island_dir_name2
+island_dir_name3 = project.island_dir_name3
+offshore_dir_name = project.offshore_dir_name
+
+# creates DEM from asc data
+print "creates DEMs from asc 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)
+convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(island_in_dir_name3, basename_out=island_dir_name3, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name,
+        use_cache=True, 
+        verbose=True)
+
+#creates pts file for island DEM
+dem2pts(island_dir_name, use_cache=True, verbose=True)
+dem2pts(island_dir_name1, use_cache=True, verbose=True)
+dem2pts(island_dir_name2, use_cache=True, verbose=True)
+dem2pts(island_dir_name3, use_cache=True, verbose=True)
+
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+print'create Geospatial data2 objects from topographies'
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
+print'create Geospatial data3 objects from topographies'
+G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+print'create Geospatial data4 objects from topographies'
+G4 = Geospatial_data(file_name = island_dir_name1 + '.pts')
+print'create Geospatial data5 objects from topographies'
+G5 = Geospatial_data(file_name = island_dir_name2 + '.pts')
+print'create Geospatial data6 objects from topographies'
+G6 = Geospatial_data(file_name = island_dir_name3 + '.pts')
+print'create Geospatial data7 objects from topographies'
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+print'create Geospatial data8 objects from topographies'
+G_off1 = Geospatial_data(file_name = offshore1_in_dir_name + '.xya')
+
+print'add all geospatial objects'
+G = G1 + G2 + G3 + G4 + G5 + G6 + G_off + G_off1 
+
+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_clipped.export_points_file(project.combined_dir_name + '.txt')
+#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 
+#-------------------------------------------------------------------------
+'''
+
+
+
+
+
+
+
Index: anuga_work/production/broome_2006/project.py
===================================================================
--- anuga_work/production/broome_2006/project.py	(revision 4152)
+++ anuga_work/production/broome_2006/project.py	(revision 4153)
@@ -7,27 +7,40 @@
 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
 
-if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME')
-    user = getenv('USERPROFILE')
+# file and system info
+#---------------------------------
+codename = 'project.py'
 
-else:    
-    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
-    user = getenv('LOGNAME')
-    print 'USER:', user
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+user = get_user_name()
 
 # INUNDATIONHOME is the inundation directory, not the data directory.
 home += sep +'data'
 
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+build_time = time+'_build'
+run_time = time+'_run'
+print 'gtime: ', gtime
+
+tide = 0
+
 #Making assumptions about the location of scenario data
 state = 'western_australia'
-scenario_dir_name = 'broome_tsunami_scenario_2006'
+scenario_name = 'broome'
+scenario = 'broome_tsunami_scenario_2006'
 
 # onshore data provided by WA DLI
 onshore_name = 'dted2_z51' # original
 
-# AHO + DPI data
+#island
+#island_name = 'rott_dli_ext' # original
+
+# offshore
+coast_name = 'coastline'
 offshore_name1 = 'XY100011610'
 offshore_name2 = 'XY100011611'
@@ -53,87 +66,99 @@
 offshore_name22 = 'XYWADPI'
 
-# developed by NM&I
-coast_name = 'coastline'
+#final topo name
+combined_name ='perth_combined_elevation'
+combined_smaller_name = 'perth_combined_elevation_smaller'
 
-# TIN model to fill in data gap
-bathy_interp = 'nearest_neighbour' #'interpolate'
 
-boundary_basename = 'SU-AU' # Mw ?
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
+topographies_time_dir = topographies_dir+build_time+sep
 
-#swollen/ all data output
-basename = 'source'
-codename = 'project.py'
+#input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+island_in_dir_name = topographies_in_dir + island_name
+island_in_dir_name1 = topographies_in_dir + island_name1
+island_in_dir_name2 = topographies_in_dir + island_name2
+island_in_dir_name3 = topographies_in_dir + island_name3
 
-#Derive subdirectories and filenames
-local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
-meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
-datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
-gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
-boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
-outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
-outputtimedir = outputdir + local_time + sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
+offshore1_in_dir_name = topographies_in_dir + offshore1_name
 
-gauge_filename = gaugedir + 'broome_gauges.csv'
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_broome = gaugedir + 'community_broome.csv'
-codedir = getcwd()+sep                            
-codedirname = codedir + 'project.py'
-meshname = outputtimedir + 'mesh_' + basename
+onshore_dir_name = topographies_dir + onshore_name
+island_dir_name = topographies_dir + island_name
+island_dir_name1 = topographies_dir + island_name1
+island_dir_name2 = topographies_dir + island_name2
+island_dir_name3 = topographies_dir + island_name3
 
-# Necessary if using point datasets, rather than grid
-onshore_dem_name = datadir + onshore_name
-offshore_interp_dem_name = datadir + bathy_interp
-offshore_dem_name1 = datadir + offshore_name1
-offshore_dem_name2 = datadir + offshore_name2
-offshore_dem_name3 = datadir + offshore_name3
-offshore_dem_name4 = datadir + offshore_name4
-offshore_dem_name5 = datadir + offshore_name5
-offshore_dem_name6 = datadir + offshore_name6
-offshore_dem_name7 = datadir + offshore_name7
-offshore_dem_name8 = datadir + offshore_name8
-offshore_dem_name9 = datadir + offshore_name9
-offshore_dem_name10 = datadir + offshore_name10
-offshore_dem_name11 = datadir + offshore_name11
-offshore_dem_name12 = datadir + offshore_name12
-offshore_dem_name13 = datadir + offshore_name13
-offshore_dem_name14 = datadir + offshore_name14
-offshore_dem_name15 = datadir + offshore_name15
-offshore_dem_name16 = datadir + offshore_name16
-offshore_dem_name17 = datadir + offshore_name17
-offshore_dem_name18 = datadir + offshore_name18
-offshore_dem_name19 = datadir + offshore_name19
-offshore_dem_name20 = datadir + offshore_name20
-offshore_dem_name21 = datadir + offshore_name21
-offshore_dem_name22 = datadir + offshore_name22
+coast_dir_name = topographies_dir + coast_name
+offshore_dir_name = topographies_dir + offshore_name
 
-coast_dem_name = datadir + coast_name
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+#combined_time_dir_final_name = topographies_time_dir + combined_final_name
 
-combined_dem_name   = datadir + 'broome_combined_elevation'
+meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
+tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
+
+
+boundaries_source = '????'
+#boundaries locations
+boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
+boundaries_in_dir_name = boundaries_in_dir + scenario_name
+boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + scenario_name
+#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+
+#output locations
+output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
+output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
+output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = 'broome_gauges.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
+community_filename = gauges_dir + 'CHINS_v2.csv'
+community_broome = gauges_dir + 'community_broome.csv'
+
+
 
 ###############################
 # Domain definitions
 ###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
-# bounding box for clipping MOST/URS output (much bigger than study area)
-##south = degminsec2decimal_degrees(-19,0,0)
-##north = degminsec2decimal_degrees(-17,15,0)
-##west  = degminsec2decimal_degrees(120,0,0)
-##east  = degminsec2decimal_degrees(122,30,0) 
-##
-##d0 = [south, west]
-##d1 = [south, east]
-##d2 = [north, east]
-##d3 = [north, west]
-##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
-##refzone = zone
+poly_all = read_polygon(polygons_dir+'extent.csv')
+res_poly_all = 100000
 
-# bounding polygon for study area
-polyAll = read_polygon(polygondir+'extent.csv')
+###############################
+# Interior region definitions
+###############################
 
-# plot bounding polygon and make sure BC info surrounds it
-#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
-print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
+poly_1 = read_polygon(polygons_dir+'Broome_Local_Polygon_update.csv')
+res_1 = 20000
+
+poly_2 = read_polygon(polygons_dir+'Broome_Close2_update.csv')
+res_2 = 1000
+
+poly_3 = read_polygon(polygons_dir+'Broome_Coast_update2.csv')
+res_3 = 1000
+#assert zone == refzone
+
+interior_regions = [[poly_1,res_1],[poly_2,res_2]
+                     ,[poly_3,res_3]]
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
 
 ###################################################################
@@ -142,44 +167,9 @@
 
 # exporting asc grid
-eastingmin = 412036.43
-eastingmax = 427184.37
-northingmin = 8007984
-northingmax = 8029830
-
-	
-	
+eastingmin = 406215.87
+eastingmax = 440208.78
+northingmin = 7983427.73
+northingmax = 8032834.52
 
 
-###############################
-# Interior region definitions
-###############################
 
-# broome digitized polygons
-poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv')
-poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv')
-poly_broome3 = read_polygon(polygondir+'Broome_Coast_update2.csv')
-#poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
-
-#plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
-print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
-print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
-print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
-#print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0
-
-for i in poly_broome3:
-    v = is_inside_polygon(i,poly_broome1, verbose=False)
-    if v == False: print v
-
-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/broome_2006/run_broome.py
===================================================================
--- anuga_work/production/broome_2006/run_broome.py	(revision 4152)
+++ anuga_work/production/broome_2006/run_broome.py	(revision 4153)
@@ -1,225 +1,201 @@
-"""Script for running a tsunami inundation scenario for Broome, WA, Australia.
+"""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.outputtimedir 
+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 tsunami wave generated by MOST.
-
-Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
+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
-import os
+from os import sep, mkdir, access, F_OK
+from os.path import dirname, basename
+from shutil import copy
 import time
-from shutil import copy
-from os.path import dirname, basename
-from os import mkdir, access, F_OK, sep
 import sys
 
 # Related major packages
-from anuga.shallow_water import Domain, Reflective_boundary, \
-                            Dirichlet_boundary, Time_boundary, File_boundary
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary, File_boundary, 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
-#-------------------------------------------------------------------------------
-
-# creates copy of code in output dir 
-copy_code_files(project.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' )
-myid = 0
-numprocs = 1
-#start_screen_catcher(project.outputtimedir, myid, numprocs)
-
-print 'USER:    ', project.user
-
-#-------------------------------------------------------------------------------
-# Preparation of topographic data
-#
-# Convert ASC 2 DEM 2 PTS using source data and store result in source data
-#-------------------------------------------------------------------------------
+#------------------------------------------------------------------------------
+
+start_screen_catcher(project.output_run_time_dir, myid, numprocs)
 
 # filenames
-onshore_dem_name = project.onshore_dem_name
-offshore_interp_dem_name = project.offshore_interp_dem_name
-coast_points = project.coast_dem_name
-meshname = project.meshname+'.msh'
-
-# creates DEM from asc data
-convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
-
-#creates pts file for onshore DEM
-dem2pts(onshore_dem_name, use_cache=True, verbose=True)
-
-# creates DEM from asc data
-convert_dem_from_ascii2netcdf(offshore_interp_dem_name, use_cache=True, verbose=True)
-
-#creates pts file for offshore interpolated DEM
-dem2pts(offshore_interp_dem_name, use_cache=True, verbose=True)
-
-print 'create offshore'
-G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\
-     Geospatial_data(file_name = project.offshore_interp_dem_name + '.pts')
-print 'create onshore'
-G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
-print 'create coast'
-G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
-print 'add'
-G = G1 + G2 + G3
-print 'export points'
-G.export_points_file(project.combined_dem_name + '.pts')
-G.export_points_file(project.combined_dem_name + '.xya')
-
-#----------------------------------------------------------------------------
-# Create the triangular mesh based on overall clipping polygon with a tagged
+boundaries_name = project.scenario
+meshes_dir_name = project.meshes_dir_name+'.msh'
+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()
+
+print 'USER: ', project.user
+print 'min triangles', project.trigs_min,
+print 'Note: This is generally about 20% less than the final amount'
+
+#--------------------------------------------------------------------------
+# 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
-#-------------------------------------------------------------------------------
-
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
-remainder_res = 500000
-local_res = 25000
-broome_res = 5000
-coast_res = 500
-interior_regions = [[project.poly_broome1, local_res],
-                    [project.poly_broome2, broome_res],
-                    [project.poly_broome3, coast_res]]
-
-from project import number_mesh_triangles
-print 'estimate of number of triangles', \
-      number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
-
+#--------------------------------------------------------------------------
+
+if myid == 0:     
+
+    print 'start create mesh from regions'
+    create_mesh_from_regions(project.poly_all,
+                         boundary_tags={'back': [1, 2], 'side': [0,3],
+                                        'ocean': [4, 5, 6]},
+                         maximum_triangle_area=project.res_poly_all,
+                         interior_regions=project.interior_regions,
+                         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
-_ = cache(create_mesh_from_regions,
-          project.polyAll,
-           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
-                              'e3': [3], 'e4':[4], 'e5': [5],
-			      'e6': [6]},
-           'maximum_triangle_area': remainder_res,
-           'filename': meshname,
-           'interior_regions': interior_regions},
-          verbose = True, evaluate=False)
-print 'created mesh'
-
-#-------------------------------------------------------------------------------                                 
-# Setup computational domain
-#-------------------------------------------------------------------------------                                 
-domain = Domain(meshname, use_cache = True, verbose = True)
-
-print 'Number of triangles = ', len(domain)
-print 'The extent is ', domain.get_extent()
-print domain.statistics()
-
-domain.set_name(project.basename)
-domain.set_datadir(project.outputtimedir)
-domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-domain.set_minimum_storable_height(0.01)
-
-#-------------------------------------------------------------------------------                                 
+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
-#-------------------------------------------------------------------------------
-
-tide = 0.0
-domain.set_quantity('stage', tide)
-domain.set_quantity('friction', 0.0) 
-domain.set_quantity('elevation', 
-                    filename = project.combined_dem_name + '.pts',
+#-------------------------------------------------------------------------
+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_dir_name+'.txt',
+#                    filename = project.combined_smaller_name_dir+'.xya',
                     use_cache = True,
                     verbose = True,
-                    alpha = 0.1
-                    )
-
-#-------------------------------------------------------------------------------                                 
-# Setup boundary conditions
-#-------------------------------------------------------------------------------
-'''
-print 'start urs2sww'
-print '', project.boundary_basename
-from anuga.shallow_water.data_manager import urs2sww
-
-south = project.south 
-north = project.north
-west  = project.west
-east  = project.east
-
-#note only need to do when an SWW file for the MOST boundary doesn't exist
-cache(urs2sww,
-      (source_dir + project.boundary_basename,
-       source_dir + project.boundary_basename), 
-      {'verbose': True,
-       'minlat': south,
-       'maxlat': north,
-       'minlon': west,
-       'maxlon': east,
-       #'origin': domain.geo_reference.get_origin(),
-       'mean_stage': tide,
-       'zscale': 1,                 #Enhance tsunami
-       'fail_on_NaN': False,
-       'inverted_bathymetry': True},
-      #evaluate = True,
-       verbose = True,
-       dependencies = source_dir + project.boundary_basename + '.sww')
-
-'''
+                    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()
 
-#Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
-#                    domain, verbose = True)
+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])
-
-# 7 min square wave starting at 1 min, 6m high
-Bw = Time_boundary(domain = domain,
-                   f=lambda t: [(60<t<480)*10, 0, 0])
-
-domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
-		      'e5': Bd,  'e6': Bd} )
-
-
-#-------------------------------------------------------------------------------                                 
+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
-#-------------------------------------------------------------------------------
-import time
+#----------------------------------------------------------------------------
+
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 240, finaltime = 480): 
+for t in domain.evolve(yieldstep = 60, finaltime = 10000):
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'e14')     
+    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)
 
-print 'finished'
