Index: /development/momentum_sink/CCS.py
===================================================================
--- /development/momentum_sink/CCS.py	(revision 2445)
+++ /development/momentum_sink/CCS.py	(revision 2446)
@@ -12,8 +12,8 @@
 from numerical_tools import norm, corr, err  
 
-swwfile_B = project.outputdir + sep  + 'Buildings_3260.sww'
 
 
-gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 20)
+
+gauge_depth = Numeric.arrayrange(46, 1.5*WidtH, 25)
 print len(gauge_depth),"gauge_depth length"
 gauge_breadth = 100
@@ -24,4 +24,6 @@
     gauge_locations.append(gauge_location)
 
+
+swwfile_B = project.outputdir + sep  + 'Straight_offset_11.sww'
 
 #Read model output
@@ -38,5 +40,5 @@
 least_norm=[1]
 norm_OK=[]
-Friction = Numeric.arrayrange(1,51,1)
+Friction = Numeric.arrayrange(1,80,1)
 for fric in Friction:
     swwfile_F = project.outputdir + sep  + 'friction_n=' + repr(fric) + '_3135.sww'
@@ -47,5 +49,5 @@
                       use_cache = True)
 
-    t=500
+    t=1000 # specifies times to look at cross sections
     
     stages_B = []
@@ -54,16 +56,17 @@
 
         # mines data from specified time values in sww file
+        # Building 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]
-
+        # Friction file
         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                   
-
         stages_B.append(wB)
         stages_F.append(wF)
@@ -73,10 +76,9 @@
     Ar_F = Numeric.array(stages_F)
     diff = abs(Ar_B - Ar_F)
-    norm_F = err(Ar_B, Ar_F, 2, relative = True)
+    norm_F = err(Ar_B, Ar_F, 2, relative = True)# 2nd norm (RMS) test
      
     print " "
     print norm_F, "norm N_F"
     print " "
-    G = corr(Ar_F, Ar_B)
     norm_OK.append(norm_F)
 
@@ -91,5 +93,5 @@
 print " "
 print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
-print fric_OK, " <- This Manning's n is OK"
+print fric_OK, " <- This Manning's n is best fit from 2nd norm test"
 print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
 print least_norm, "2nd normal of Friction vs Buildings  "
Index: /development/momentum_sink/PCS.py
===================================================================
--- /development/momentum_sink/PCS.py	(revision 2445)
+++ /development/momentum_sink/PCS.py	(revision 2446)
@@ -11,14 +11,15 @@
 from pylab import *
 
-fric = 40    
+fric = 25    
 #swwfile = project.newoutputname + '.sww'
 #swwfile = project.outputname
 
 #swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
-swwfile_B = project.outputdir + sep  + 'Buildings_3662.sww'
+swwfile_B = project.outputdir + sep  + 'Straight_offset.sww'
 #swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
-swwfile_F = project.outputdir + sep  + 'friction_n=40_3135.sww'
+swwfile_F = project.outputdir + sep  + 'friction_n=51_3135.sww'
 
-gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 5)
+gauge_depth = Numeric.arrayrange(46, 1.5*WidtH, 25) # Random special
+#gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 25)
 gauge_breadth = 100
 gauge_locations = []
@@ -45,5 +46,5 @@
 
 
-T = Numeric.arrayrange(0,1000,20)
+T = Numeric.arrayrange(0,1000,10)
 #T = [ 20, 50, 100, 150, 200 ]
 
Index: /development/momentum_sink/create_buildings.py
===================================================================
--- /development/momentum_sink/create_buildings.py	(revision 2445)
+++ /development/momentum_sink/create_buildings.py	(revision 2446)
@@ -2,12 +2,16 @@
 import Numeric
 import math
+import random
 # add inundation dir to your pythonpath
 from pmesh.mesh import Mesh
 from coordinate_transforms.geo_reference import Geo_reference
 
+depth = 11 # depth of building side to oncoming wave
+breadth = 11 # breadth of building, width of building to oncoming wave
+#from loop_buildings import depth
 
 WidtH = 200 # width of boudary in metres
-depth = 15 # depth of building side to oncoming wave
-breadth = 15 # breadth of building
+#breadth = depth
+
 print "building footprint"
 print depth * breadth , "m^2"
@@ -15,6 +19,8 @@
 BL = block**0.5
 
-porosity = breadth/BL
+porosity = breadth/BL				
 print porosity, " Building porosity"
+
+#random.uniform(0,1) #proper implemantation of random generator
 
 def create_mesh(maximum_triangle_area,
@@ -25,12 +31,11 @@
     triangles in the mesh to the mesh file name.
     """
+
+    
     # create a mesh instance of class Mesh
     m = Mesh()
+
     # Boundary of problem
-    
-    #W = WidtH/8
-    #L = W
     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,2], 'wave':[3], 'back':[1]})
 
@@ -41,5 +46,5 @@
     #Th = (45 *(3.142/180)) # sets an initial rotation      
     Th = 0
-    ForDep = (0.2*WidtH) + (BL/2)
+    ForDep = (0.2*WidtH) + (BL-whs)
     RearDep = 1.2*WidtH 
         
@@ -61,9 +66,7 @@
             lh4 = (+lhs) * math.cos(Th) - (-whs) * math.sin(Th)
             polygon = [[D+wh1,B+lh1],[D+wh2,B+lh2],[D+wh3,B+lh3],[D+wh4,B+lh4]]
-            #polygon = [[D-wh,B-lh],[D+wh,B-lh],[D+wh,B+lh],[D-wh,B+lh]]
-            m.add_hole_from_polygon(polygon, tags={'wall':[0,1,2,3]})# Adds holes with reflective boundaries.        
-            Th = Th + (37.3 *(3.142/180)) # keeps rotating individual buildings.
+            m.add_hole_from_polygon(polygon, tags={'wall':[0,1,2,3]})        
+            #Th = Th + (90 *(3.14159/180)) # keeps rotating individual buildings.
 
-            #print polygon
     m.generate_mesh(maximum_triangle_area=maximum_triangle_area)
     triangle_count = m.get_triangle_count()
@@ -73,5 +76,5 @@
     else:
         if triangles_in_name is True:
-            mesh_file = mesh_file[:-4] + '_' + str(triangle_count) \
+            mesh_file = mesh_file[:-4] + '_Ran_' + str(depth) + '_' + str(triangle_count) \
                         + mesh_file[-4:]
         m.export_mesh_file(mesh_file)
Index: /development/momentum_sink/loop_buildings.py
===================================================================
--- /development/momentum_sink/loop_buildings.py	(revision 2446)
+++ /development/momentum_sink/loop_buildings.py	(revision 2446)
@@ -0,0 +1,78 @@
+
+#Convention for strings representing files
+# #_file has the extention
+#           #name does not have the extension
+
+import time
+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
+from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
+
+import Numeric
+
+
+
+#from building_generator import create_mesh
+
+
+
+DR = [3,5,7,9]
+for depth in DR:
+    import project
+    
+    from create_buildings import create_mesh
+    from pmesh.mesh import importMeshFromFile
+    meshname = project.meshname
+    outputname = project.outputname
+    t0 = time.time()
+    meshname, triagle_count = cache(create_mesh,(100),
+                                    {'mesh_file':meshname,
+                                     'triangles_in_name':True}
+                                    ,dependencies = ['create_buildings.py']
+                                    #,evaluate = True      
+                                    )
+    
+    print 'Initialising the mesh took %.2f seconds' %(time.time()-t0)  
+
+    #Setup domain
+    domain = cache(pmesh_to_domain_instance, (meshname, Domain),
+                   dependencies = [meshname]                    
+                   ,verbose = False
+                   )               
+
+
+    domain.set_name(project.basename + '_brd=%dm_%d' %(depth, triagle_count))
+    domain.set_datadir(project.outputdir)
+    domain.store = True
+    domain.quantities_to_be_stored = ['stage']
+
+    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)
+
+    #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<300)*4, 0, 0])
+
+    domain.set_boundary( {'wall': Br,'wave': Bdw, 'back': Bdb, 'exterior':Bdw} ) 
+
+    #Evolve
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 1, finaltime = 10):
+        domain.write_time()      
+
+    print 'That took %.2f seconds' %(time.time()-t0)
Index: /development/momentum_sink/loop_friction.py
===================================================================
--- /development/momentum_sink/loop_friction.py	(revision 2445)
+++ /development/momentum_sink/loop_friction.py	(revision 2446)
@@ -24,5 +24,5 @@
 from pmesh.mesh import importMeshFromFile
 
-Friction = Numeric.arrayrange(1,50,1)
+Friction = Numeric.arrayrange(50,80,1)
 for fric in Friction:
     meshname = project_friction.meshname
Index: /development/momentum_sink/run_buildings.py
===================================================================
--- /development/momentum_sink/run_buildings.py	(revision 2445)
+++ /development/momentum_sink/run_buildings.py	(revision 2446)
@@ -15,5 +15,5 @@
 
 
-from create_buildings import create_mesh
+from create_buildings import create_mesh, depth
 #from building_generator import create_mesh
 
@@ -46,5 +46,5 @@
 
 
-domain.set_name(project.basename + '_%d' %triagle_count)
+domain.set_name(project.basename + '_Br=%d_%d' %(depth, triagle_count))
 domain.set_datadir(project.outputdir)
 domain.store = True
@@ -75,5 +75,5 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 10, finaltime = 500):
+for t in domain.evolve(yieldstep = 10, finaltime = 1000):
     domain.write_time()      
 
