Index: /development/momentum_sink/compare_sww.py
===================================================================
--- /development/momentum_sink/compare_sww.py	(revision 2394)
+++ /development/momentum_sink/compare_sww.py	(revision 2395)
@@ -5,11 +5,24 @@
 import project
 from pyvolution.util import file_function
+import Numeric
+#from friction_block #import test
+
+#print "test", test
+#import sys
+#sys.exit() 
 
 
 #swwfile = project.outputname + '.sww'
-swwfile = project.outputdir + sep  + 'friction_12330.sww'
-gauge_locations = [[10,10],[20,20]]
+#swwfile = project.outputdir + sep  + 'FrictionBlock_0.01_3.sww'
+swwfile = project.outputdir + sep  + 'Buildings_3662.sww'
+gauge_depth = Numeric.arrayrange(0, 700, 20)
+gauge_breadth = 100
+gauge_locations = []
 
-quantities = ['stage','elevation']
+for GD in gauge_depth:
+    gauge_location = [GD,gauge_breadth]
+    gauge_locations.append(gauge_location)
+
+quantities = ['stage','elevation','xmomentum']
 gauges = file_function(swwfile,
                   quantities = quantities,
@@ -18,5 +31,5 @@
                   use_cache = True)
 
-times = [0,1,2,3,4]
+times = [150,500]
 
 for gauge, _ in enumerate(gauge_locations):
Index: /development/momentum_sink/create_buildings.py
===================================================================
--- /development/momentum_sink/create_buildings.py	(revision 2394)
+++ /development/momentum_sink/create_buildings.py	(revision 2395)
@@ -18,15 +18,15 @@
 
     # Boundary of problem
-    WidtH = 500 # width of boudary in metres
+    WidtH = 200 # width of boudary in metres
     #W = WidtH/8
     #L = W
-    outer_polygon = [[0,0],[1.5*WidtH,0],[1.5*WidtH,WidtH],[0,WidtH]]
+    outer_polygon = [[0,0],[5*WidtH,0],[5*WidtH,WidtH],[0,WidtH]]
     print outer_polygon
     m.add_region_from_polygon(outer_polygon, tags={'wall':[0,1,2], 'wave':[3]})
 
     # inner polygons => building boundaries
-    depth = 18 # depth of building side to oncoming wave
+    depth = 20 # depth of building side to oncoming wave
     wh = depth/2 
-    breadth = 18 # breadth of building
+    breadth = 20 # breadth of building
     lh = breadth/2 
       
@@ -46,5 +46,5 @@
 
     for D in Depths:
-        Breadths = Breadths + BL/2
+        #Breadths = Breadths + BL/2 #Used to offset buildings
         for B in Breadths:
             polygon = [[D-wh,B-lh],[D+wh,B-lh],[D+wh,B+lh],[D-wh,B+lh]]
Index: /development/momentum_sink/friction_block.py
===================================================================
--- /development/momentum_sink/friction_block.py	(revision 2394)
+++ /development/momentum_sink/friction_block.py	(revision 2395)
@@ -6,5 +6,5 @@
 from coordinate_transforms.geo_reference import Geo_reference
 
-
+#test = 300
 def create_mesh(maximum_triangle_area,
                 mesh_file=None,
@@ -16,5 +16,6 @@
     # create a mesh instance of class Mesh
     m = Mesh()
-
+    
+    #print "test", test  
     # Boundary of problem
     WidtH = 200 # width of boudary in metres
@@ -28,5 +29,5 @@
     #Boundary
     dict = {}
-    dict['points'] = [[0,0],[0.2*WidtH,0],[1.2*WidtH,0],[3.5*WidtH,0],[3.5*WidtH,WidtH] \
+    dict['points'] = [[0,0],[0.2*WidtH,0],[1.2*WidtH,0],[5*WidtH,0],[5*WidtH,WidtH] \
                       ,[1.2*WidtH,WidtH],[0.2*WidtH,WidtH],[0,WidtH]]
     
@@ -37,5 +38,5 @@
                             'wall',                           
                             'wall',
-                            'wall',
+                            'back',
                             'wall',
                             'wall',
@@ -50,7 +51,9 @@
     m.generate_mesh(maximum_triangle_area=maximum_triangle_area)
     triangle_count = m.get_triangle_count()
+  
+
     
     if mesh_file is None:    
-        return m, triangle_count
+        return m, triangle_count, Width
     else:
         if triangles_in_name is True:
Index: /development/momentum_sink/process_gauges.py
===================================================================
--- /development/momentum_sink/process_gauges.py	(revision 2395)
+++ /development/momentum_sink/process_gauges.py	(revision 2395)
@@ -0,0 +1,96 @@
+"""Read in sww file, interpolate at specified locations and plot time series
+Additionally, store augmented building database with max inundation depths
+"""
+
+import project
+from pyvolution.util import file_function
+from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+from pylab import *
+
+
+plot = False
+    
+swwfile = project.outputname + '.sww'
+
+def get_gauges_from_file(filename):
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+
+    gauges = []
+    gaugelocation = []
+    for line in lines[1:]:
+        fields = line.split(',')
+        location = fields[0]
+        easting = float(fields[1])
+        northing = float(fields[2])
+     
+        gauges.append([easting, northing])
+        gaugelocation.append(location)
+
+    #Return gauges and raw data for subsequent storage    
+    return gauges, lines, gaugelocation
+
+gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
+
+
+#Read model output
+quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+f = file_function(swwfile,
+                  quantities = quantities,
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+lines[0] = lines[0].strip() +\
+               ',MAX INUNDATION DEPTH (m)' +\
+               ',MAX MOMENTUM (m^2/s)\n'            
+               
+from math import sqrt
+N = len(gauges)
+for k, g in enumerate(gauges):
+    if k%((N+10)/10)==0:
+        print 'Doing row %d of %d' %(k, N)
+
+    model_time = []
+    stages = []
+    elevations = []
+    momenta = []
+
+    max_depth = 0
+    max_momentum = 0    
+    for i, t in enumerate(f.T):
+        w = f(t, point_id = k)[0]
+        z = f(t, point_id = k)[1]
+        uh = f(t, point_id = k)[2]
+        vh = f(t, point_id = k)[3]        
+
+        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
+        
+        model_time.append(t)        
+        stages.append(w)
+        elevations.append(z)  #Should be constant over time
+        momenta.append(m) 
+
+        if w-z > max_depth:
+            max_depth = w-z
+        if m > max_momentum:            
+            max_momentum = m    
+
+    #Augment building data
+    lines[k+1] = lines[k+1].strip() +\
+                     ',%f' %max_depth +\
+                     ',%f\n' %max_momentum
+    
+    print lines[k]
+    
+#Store new building file with mak depths added
+
+#FN = 'inundation_augmented_' + project.gauge_filename
+FN = project.gauge_outname
+fid = open(FN, 'w')
+for line in lines:
+    fid.write(line)
+fid.close()    
+
Index: /development/momentum_sink/run_buildings.py
===================================================================
--- /development/momentum_sink/run_buildings.py	(revision 2394)
+++ /development/momentum_sink/run_buildings.py	(revision 2395)
@@ -65,14 +65,15 @@
 Br = Reflective_boundary(domain)
 Bt = Transmissive_boundary(domain)
-Bd = Dirichlet_boundary([10,0,0])
+Bdw = Dirichlet_boundary([36,0,0])
+Bdb = Dirichlet_boundary([0,0,0])
 Bw = Time_boundary(domain=domain,
                    f=lambda t: [(60<t<660)*4, 0, 0])
 
-domain.set_boundary( {'wall': Br,'wave': Bd, 'back': Bt, 'exterior':Bd} ) 
+domain.set_boundary( {'wall': Br,'wave': Bdw, 'back': Bdb, 'exterior':Bdw} ) 
 
 #Evolve
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 1, finaltime = 20):
+for t in domain.evolve(yieldstep = 10, finaltime = 1000):
     domain.write_time()      
 
Index: /development/momentum_sink/run_friction.py
===================================================================
--- /development/momentum_sink/run_friction.py	(revision 2394)
+++ /development/momentum_sink/run_friction.py	(revision 2395)
@@ -26,4 +26,6 @@
 t0 = time.time()
 
+#meshname, triagle_count =create_mesh(100,mesh_file=meshname,triangles_in_name=True)
+
 meshname, triagle_count = cache(create_mesh,(100),
                                 {'mesh_file':meshname,
@@ -50,5 +52,5 @@
 domain.set_datadir(project_friction.outputdir)
 domain.store = True
-domain.quantities_to_be_stored = ['stage']
+domain.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
 
 print 'Number of triangles = ', len(domain)
@@ -58,5 +60,5 @@
 domain.set_quantity('friction', 0.01)
 domain.set_quantity('stage', 0)
-domain.set_region(Set_region('mound', 'friction', 0.8)) #, location='unique vertices'))
+domain.set_region(Set_region('mound', 'friction', 3)) #, location='unique vertices'))
 #Setup Boundary Conditions
 print domain.get_boundary_tags()
@@ -66,14 +68,15 @@
 Br = Reflective_boundary(domain)
 Bt = Transmissive_boundary(domain)
-Bd = Dirichlet_boundary([36,0,0])
+Bdw = Dirichlet_boundary([36,0,0])
+Bdb = Dirichlet_boundary([0,0,0])
 Bw = Time_boundary(domain=domain,
                    f=lambda t: [(60<t<660)*4, 0, 0])
 
-domain.set_boundary( {'wall': Br,'wave': Bd, 'back': Bt, 'exterior':Bd} ) 
+domain.set_boundary( {'wall': Br,'wave': Bdw, 'back': Bdb, 'exterior':Bdw} ) 
 
 #Evolve
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 1, finaltime = 80):
+for t in domain.evolve(yieldstep = 4, finaltime = 400):
     domain.write_time()      
 
Index: /development/momentum_sink/run_rotate.py
===================================================================
--- /development/momentum_sink/run_rotate.py	(revision 2394)
+++ /development/momentum_sink/run_rotate.py	(revision 2395)
@@ -74,5 +74,5 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 1, finaltime = 20):
+for t in domain.evolve(yieldstep = 1, finaltime = 60):
     domain.write_time()      
 
