Index: /development/momentum_sink/CCS.py
===================================================================
--- /development/momentum_sink/CCS.py	(revision 2417)
+++ /development/momentum_sink/CCS.py	(revision 2417)
@@ -0,0 +1,101 @@
+"""Read in sww file, interpolate at specified locations and
+    Cross section at specified T
+
+"""
+from os import sep
+import Numeric
+import project
+from pyvolution.util import file_function
+from create_buildings import porosity, WidtH
+# from run_friction import fric
+from pylab import *
+    
+
+swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
+swwfile_F = project.outputdir + sep  + 'friction_n=20_3135.sww'
+
+gauge_depth = Numeric.arrayrange(0, 2.5*WidtH, 5)
+gauge_breadth = 100
+gauge_locations = []
+
+for GD in gauge_depth:
+    gauge_location = [GD,gauge_breadth]
+    gauge_locations.append(gauge_location)
+
+
+#Read model output
+quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+
+f_B = file_function(swwfile_B,
+                  quantities = quantities,
+                  interpolation_points = gauge_locations,
+                  verbose = True,
+                  use_cache = True)
+
+f_F = file_function(swwfile_F,
+                  quantities = quantities,
+                  interpolation_points = gauge_locations,
+                  verbose = True,
+                  use_cache = True)
+
+
+T = Numeric.arrayrange(0,500,20)
+#T = [ 20, 50, 100, 150, 200 ]
+
+for t in T:
+    model_time = []
+    stages_B = []
+    stages_F = []
+    elevations = []
+    momenta = []
+    velocity = []
+    for i, GL in enumerate(gauge_depth):
+
+        # mines data from specified time values in sww file
+        wB = f_B(t, point_id = i)[0]
+        zB = f_B(t, point_id = i)[1]
+        uhB = f_B(t, point_id = i)[2]
+        vhB = f_B(t, point_id = i)[3]
+
+        wF = f_F(t, point_id = i)[0]        
+        uhF = f_B(t, point_id = i)[2]
+        vhF = f_B(t, point_id = i)[3]
+        #m = sqrt(uh*uh + vh*vh)   #Absolute momentum
+        velB = sqrt(uhB*uhB + vhB*vhB) / (wB + 1.e-30) #Absolute velocity
+        velF = sqrt(uhF*uhF + vhF*vhF) / (wF + 1.e-30) #Absolute velocity                   
+        #print velB, "buiding velocity"
+        #print velF, "friction velocity"
+
+        stages_B.append(wB)
+        stages_F.append(wF)
+        #elevations.append(zB)  #Should be constant over time
+        #momenta.append(mB)
+        #velocityB.append(velB)
+        #velocityF.append(velF)
+        
+    #print len(stages), "stages"
+    diff = abs(Numeric.array(stages_B) - Numeric.array(stages_F))
+    #print diff, "difference"
+
+
+    ion()
+    hold(True)
+    subplot(211)
+    
+    plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r')
+    title('Time_%d_ vs water depth (n=%f, BP=%f)' %(t, fric, porosity))
+    xlabel('gauge length(m)')
+    ylabel('water depths (m)')
+    #savefig('Comp_Time%d_B_vs_F' %t)
+    hold(False)
+#savefig('Final_comp_n=%s' %fric)
+subplot(212)
+subplots_adjust(hspace = 0.4)
+plot(gauge_depth, diff, 'g')
+title('Final Differences in water levels')
+xlabel('gauge length(m)')
+ylabel('water differences (m)')
+savefig('Final_Water_diff_n=%s' %fric)
+print "finished"
+
+
Index: /development/momentum_sink/PCS.py
===================================================================
--- /development/momentum_sink/PCS.py	(revision 2416)
+++ /development/momentum_sink/PCS.py	(revision 2417)
@@ -7,17 +7,18 @@
 import project
 from pyvolution.util import file_function
-#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
-#from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+from create_buildings import porosity, WidtH
+# from run_friction import fric
 from pylab import *
-#from compare_sww import gauge_locations
 
-    
+fric = 1    
 #swwfile = project.newoutputname + '.sww'
 #swwfile = project.outputname
-swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
-swwfile_F = project.outputdir + sep  + 'Buildings_3662.sww'
-#swwfile_F = project.outputdir + sep  + 'friction_3135.sww'
 
-gauge_depth = Numeric.arrayrange(0, 700, 10)
+#swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
+#swwfile_B = project.outputdir + sep  + 'Buildings_3662.sww'
+swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
+swwfile_F = project.outputdir + sep  + 'friction_n=20_3135.sww'
+
+gauge_depth = Numeric.arrayrange(0, 2.5*WidtH, 5)
 gauge_breadth = 100
 gauge_locations = []
@@ -55,4 +56,6 @@
     velocity = []
     for i, GL in enumerate(gauge_depth):
+
+        # mines data from specified time values in sww file
         wB = f_B(t, point_id = i)[0]
         zB = f_B(t, point_id = i)[1]
@@ -82,13 +85,21 @@
 
     ion()
-    hold(False)
+    hold(True)
+    subplot(211)
+    
     plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r')
-    title('Time_%d_ vs water depth' %t)
+    title('Time_%d_ vs water depth (n=%f, BP=%f)' %(t, fric, porosity))
     xlabel('gauge length(m)')
     ylabel('water depths (m)')
-    savefig('Comp_Time%d_B_vs_F' %t)
-    #plot(gauge_depth, diff, 'g')
-    #title('differences')
-    # raw_input('Next')
+    #savefig('Comp_Time%d_B_vs_F' %t)
+    hold(False)
+#savefig('Final_comp_n=%s' %fric)
+subplot(212)
+subplots_adjust(hspace = 0.4)
+plot(gauge_depth, diff, 'g')
+title('Final Differences in water levels')
+xlabel('gauge length(m)')
+ylabel('water differences (m)')
+savefig('Final_Water_diff_n=%s' %fric)
 print "finished"
 
Index: /development/momentum_sink/create_buildings.py
===================================================================
--- /development/momentum_sink/create_buildings.py	(revision 2416)
+++ /development/momentum_sink/create_buildings.py	(revision 2417)
@@ -7,4 +7,14 @@
 
 
+WidtH = 200 # width of boudary in metres
+depth = 20 # depth of building side to oncoming wave
+breadth = 20 # breadth of building
+print "building footprint"
+print depth * breadth , "m^2"
+block = 625
+BL = block**0.5
+
+porosity = breadth/BL
+print porosity, " Building porosity"
 
 def create_mesh(maximum_triangle_area,
@@ -18,5 +28,5 @@
     m = Mesh()
     # Boundary of problem
-    WidtH = 200 # width of boudary in metres
+    
     #W = WidtH/8
     #L = W
@@ -26,17 +36,11 @@
 
     # inner polygons => building boundaries
-    depth = 4 # depth of building side to oncoming wave
+    
     wh = depth/2 
-    breadth = 4 # breadth of building
+    
     lh = breadth/2 
           
-    print "building footprint"
-    print depth * breadth , "m^2"
-    block = 625
-    BL = block**0.5
     ForDep = (0.2*WidtH) + (BL/2)
     RearDep = 1.2*WidtH 
-    porosity = breadth/BL
-    print porosity, " Building porosity"
         
     Breadths = Numeric.arrayrange( (BL/2), WidtH, (BL))
Index: /development/momentum_sink/loop_friction.py
===================================================================
--- /development/momentum_sink/loop_friction.py	(revision 2417)
+++ /development/momentum_sink/loop_friction.py	(revision 2417)
@@ -0,0 +1,90 @@
+
+#Convention for strings representing files
+# #_file has the extention
+#           #name does not have the extension
+
+
+# set friction value for zone
+
+# fwojhrt
+import time
+
+import project_friction
+
+from pyvolution.pmesh2domain import pmesh_to_domain_instance 
+from caching import cache
+from pyvolution.shallow_water import Domain, Reflective_boundary,\
+     File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary,\
+     Set_region
+from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
+
+
+from friction_block import create_mesh
+#from building_generator import create_mesh
+
+from pmesh.mesh import importMeshFromFile
+
+Friction = [1,2,3,5,7,10,20,30,50]
+for fric in Friction:
+    meshname = project_friction.meshname
+    outputname = project_friction.outputname
+
+    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,
+                                     'triangles_in_name':True}
+                                    ,dependencies = ['friction_block.py']
+                                    #,evaluate = True      
+                                    )
+    #meshname = 'build.tsh'
+    #outputname = outputname[:-4] + '_' + str(triagle_count) + outputname[-4:]
+
+    print 'Initialising the mesh took %.2f seconds' %(time.time()-t0)  
+
+
+    #meshname = importMeshFromFile('build.tsh')
+
+    #Setup domain
+    domain = cache(pmesh_to_domain_instance, (meshname, Domain),
+                   dependencies = [meshname]                    
+                   ,verbose = False
+                   )               
+
+
+    domain.set_name(project_friction.basename + '%d_%d' %(fric, triagle_count))
+    domain.set_datadir(project_friction.outputdir)
+    domain.store = True
+    domain.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
+
+    print 'Number of triangles = ', len(domain)
+    print 'The extent is ', domain.get_extent()
+
+    #Setup Initial Conditions
+    domain.set_quantity('friction', 0.01)
+    domain.set_quantity('stage', 0)
+    domain.set_region(Set_region('mound', 'friction', fric)) #, location='unique vertices'))
+    #Setup Boundary Conditions
+    print domain.get_boundary_tags()
+
+    domain.starttime = 0  #Obtained from MOST
+
+    Br = Reflective_boundary(domain)
+    Bt = Transmissive_boundary(domain)
+    Bdw = Dirichlet_boundary([20,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': Bdw, 'back': Bdb, 'exterior':Bdw} ) 
+
+    #Evolve
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 10, finaltime = 500):
+        domain.write_time()      
+
+    print 'That took %.2f seconds' %(time.time()-t0)
+
Index: /development/momentum_sink/project_friction.py
===================================================================
--- /development/momentum_sink/project_friction.py	(revision 2416)
+++ /development/momentum_sink/project_friction.py	(revision 2417)
@@ -12,5 +12,5 @@
 
 scenario_dir_name = 'momentum_sink_2005'
-basename = 'friction'
+basename = 'friction_n='
 
 if sys.platform == 'win32':
Index: /development/momentum_sink/run_friction.py
===================================================================
--- /development/momentum_sink/run_friction.py	(revision 2416)
+++ /development/momentum_sink/run_friction.py	(revision 2417)
@@ -4,4 +4,9 @@
 #           #name does not have the extension
 
+
+# set friction value for zone
+fric=10
+
+# fwojhrt
 import time
 
@@ -60,5 +65,5 @@
 domain.set_quantity('friction', 0.01)
 domain.set_quantity('stage', 0)
-domain.set_region(Set_region('mound', 'friction', 10)) #, location='unique vertices'))
+domain.set_region(Set_region('mound', 'friction', fric)) #, location='unique vertices'))
 #Setup Boundary Conditions
 print domain.get_boundary_tags()
