Index: inundation/ga/storm_surge/examples/beach.py
===================================================================
--- inundation/ga/storm_surge/examples/beach.py	(revision 668)
+++ inundation/ga/storm_surge/examples/beach.py	(revision 672)
@@ -23,6 +23,6 @@
 from util import read_polygon, Polygon_function
 from math import pi
-from Numeric import choose, greater, ones, sin
-
+from Numeric import choose, greater, ones, sin, exp
+import time
 
 ######################
@@ -33,4 +33,5 @@
 
 #HACK to remove internal boundary (until we get a solution)
+#Is this still necessary?
 for b in domain.boundary.keys():
     if domain.boundary[b] == '':
@@ -46,18 +47,29 @@
 
 def bathymetry(x, y):
-    #cut = 55
-    #return -(x - cut)/5
-
     cut = 75
-    return choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) )) 
-    ####return choose( greater(x, cut), (-20*ones(x.shape), ones(x.shape)))   
+    
+    res = choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) )) 
+    res += (100-y)/50 - 0.75
+    return res
     
 
 def topography(x, y):
-    return 4*sin(pi*x/50) + 1
+
+    import RandomArray
+    #z = 4*sin(pi*x/50) + (100-y)/50 + 1 #Ridge
+    
+    z = -4.0*x/25 + 8  + (100-y)/50 #Beach
+    
+    #z += sin(pi*x/5)/4 
+    #z += RandomArray.normal(0, 1, x.shape)/4 
+    
+    z += 6*exp( -((x-30)/10)**2 ) * exp( -((y-50)/8)**2 )     #Mound
+    z += 4*exp( -((x-30)/4)**2 ) * exp( -((y-26)/10)**2 )     #extra ridge    
+    z += 4*exp( -((x-30)/5)**2 ) * exp( -((y-10)/8)**2 )     #extra ridge        
+    
+    return z
 
 def riverbed(x, y):
-    #return (y-100)/10 + 5*x/100
-    return (y-100)/50 + 4.4
+    return (y-100)/70 - 4.0*x/25 + 8
     
 
@@ -66,15 +78,13 @@
 
 land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
-        [0,0], [0,100]]        
+        [0,0], [0,100]] 
+		       
 	     
-building1 = [[45,80], [50,78], [52,83], [45,83]]	     
+building1 = [[45,80], [49,78], [52,83], [46,83]]	     
 building2 = [[35,75], [40,75], [40,80], [35,80]]	     
-building3 = [[42,63], [47,63], [50,68], [37,68]]	     
-building4 = [[28,70], [28,60], [33,60], [33,70]]	     
-building5 = [[10,70], [10,60], [15,60], [15,70]]	     
-building6 = [[10,60], [10,50], [15,50], [15,60]]	     
-
-
-
+building3 = [[42,63], [46,61], [48,65], [44,67]]	     
+building4 = [[28,56], [28,51], [33,51], [33,56]]	     
+building5 = [[10,70], [10,65], [15,65], [15,70]]	     
+building6 = [[10,50], [10,45], [15,45], [15,50]]	     
 
 river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
@@ -82,10 +92,13 @@
 	      
 print 'Set elevation'    
+t0 = time.time()
 domain.set_quantity('elevation',
     Polygon_function( [(shoreline, bathymetry), (land, topography), 
-                      (building1, 3), (building2, 9), 
-		      (building3, 8), (building4, 7),
-		      (building5, 8), (building6, 6),
+                      (building1, 6), (building2, 8), 
+		      (building3, 6), (building4, 13),
+		      (building5, 10), (building6, 11),
 		      (river, riverbed)]))
+		      
+print 'That took %.2f seconds' %(time.time()-t0)
 		      
 print 'Set level'    		      
@@ -95,5 +108,5 @@
 		    
 print 'Set friction'    		      		    
-domain.set_quantity('friction', 0.03)
+domain.set_quantity('friction', 0.075)
 
 print domain.get_extent()
@@ -111,5 +124,5 @@
 
     if tt > 0.7:
-        return 10000*tt*ones(N, Float)
+        return 30000*tt*ones(N, Float)
     else:
         return zeros(N, Float)
@@ -125,7 +138,7 @@
 Bo = Transmissive_boundary(domain)
 
-#Constant inflow
-Bd = Dirichlet_boundary([0.0, -1.0, 0.0])
-Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])
+#Constant outflow
+Bd = Dirichlet_boundary([-10, 0.0, 0.0])
+Bt = Time_boundary(domain, lambda t: [ 2.8*(1+sin(2*pi*t/100)), 0.0, 0.0])
 
 print 'Available boundary tags are', domain.get_boundary_tags()
@@ -135,7 +148,12 @@
 tags['ocean'] = Bt
 tags['wall'] = Br 
-tags['outflow'] = Bo 
+tags['wall1'] = Br 
+tags['outflow'] = Bd 
 tags['exterior'] = Br
-tags['external'] = Br 
+tags['external'] = Br
+tags['land'] = Bo  
+tags['westbank'] = None    #Riverbank
+tags['eastbank'] = None    #Riverbank
+tags['eastbankN'] = None    #Riverbank
 domain.set_boundary(tags)
 
@@ -145,5 +163,5 @@
 ######################
 #Evolution
-for t in domain.evolve(yieldstep = 0.2, finaltime = 500):
+for t in domain.evolve(yieldstep = 0.2, finaltime = 300):
     domain.write_time()
     
