Index: undation/ga/storm_surge/analytical solutions/Analytical_solution_Samnpson_import_mesh.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/Analytical_solution_Samnpson_import_mesh.py	(revision 557)
+++ 	(revision )
@@ -1,109 +1,0 @@
-"""Example of shallow water wave equation analytical solution
-consists of a parabolic profile in a parabolic basin. Analytical
-solutiuon to this problem was derived by Carrier and Greenspan
-and used by Toon and Chou.
-The analytical solution was derived cy Sampson in 2002.
-
-   Copyright 2004
-   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
-   Geoscience Australia
-   
-"""
-
-######################
-# Module imports 
-#
-
-import sys
-from os import sep
-sys.path.append('..'+sep+'pyvolution')
-
-from shallow_water import Domain, Dirichlet_boundary
-from math import sqrt, cos, sin, pi
-from mesh_factory import strang_mesh
-
-
-######################
-# Domain
-# Strang_domain will search through the file and test to see if there are
-# two or three entries. Two entries are for points and three for triangles.
-
-points, elements = strang_mesh('yoon_circle.pt')
-domain = Domain(points, elements) 
-
-domain.default_order = 2
-domain.smooth = True
-
-
-# Provide file name for storing output
-domain.store = True
-domain.format = 'sww'
-domain.filename = 'yoon_strang_second_order'
-
-print 'Number of triangles = ', len(domain)
-
-#Reduction operation for get_vertex_values
-from util import mean
-domain.reduction = mean 
-#domain.reduction = min  #Looks better near steep slopes
-
-
-######################
-#Initial condition
-#
-print 'Initial condition'
-t = 0.0
-D0 = 1.
-L = 2500.
-R0 = 2000.
-g = 9.81
-
-A = (L**4 - R0**4)/(L**4 + R0**4)
-omega = 2./L*sqrt(2.*g*D0)
-print 'omega = ',omega
-T = pi/omega
-print T
-
-
-#Set bed-elevation and friction(None)
-def x_slope(x,y):
-    n = x.shape[0]
-    z = 0*x
-    for i in range(n):
-        r = sqrt(x[i]*x[i] + y[i]*y[i])
-        z[i] = -D0*(1.-r*r/L/L)
-    return z
-
-domain.set_quantity('elevation', x_slope)
-
-#Set the water level (w = z+h)
-def level(x,y):
-    z = x_slope(x,y)
-    n = x.shape[0]
-    h = 0*x
-    for i in range(n):
-        r = sqrt(x[i]*x[i] + y[i]*y[i])
-        h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
-                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
-    if h[i] < z[i]:
-        h[i] = z[i]
-    return h   
-
-domain.set_quantity('level', level)
-
-
-############
-#Boundary
-domain.set_boundary({'exterior': Dirichlet_boundary([0.0, 0.0, 0.0])})
-    
-
-######################
-#Evolution
-import time
-t0 = time.time()
-for t in domain.evolve(yieldstep = 10.0, finaltime = 300):
-    domain.write_time()
-
-print 'That took %.2f seconds' %(time.time()-t0)
-
-    
Index: /inundation/ga/storm_surge/analytical solutions/Analytical_solution_Sampson_import_mesh.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/Analytical_solution_Sampson_import_mesh.py	(revision 557)
+++ /inundation/ga/storm_surge/analytical solutions/Analytical_solution_Sampson_import_mesh.py	(revision 558)
@@ -26,4 +26,5 @@
 from math import sqrt, cos, sin, pi, exp
 from mesh_factory import strang_mesh
+from quantity import Quantity
 
 
@@ -39,4 +40,9 @@
 domain.default_order = 2
 domain.smooth = True
+
+domain.quantities['linear_friction'] = Quantity(domain)
+domain.set_quantity('friction', 0.0)
+#del domain.quantities['friction']
+
 
 # Provide file name for storing output
@@ -79,4 +85,6 @@
 
 domain.set_quantity('elevation', x_slope)
+domain.set_quantity('linear_friction', tau)
+
 
 #Set the water level
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 557)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 558)
@@ -877,5 +877,5 @@
     uh = domain.quantities['xmomentum'].centroid_values
     vh = domain.quantities['ymomentum'].centroid_values
-    eta = domain.quantities['linear_friction'].centroid_values    
+    tau = domain.quantities['linear_friction'].centroid_values    
     
     xmom_update = domain.quantities['xmomentum'].semi_implicit_update
@@ -886,7 +886,7 @@
     
     for k in range(N):
-        if eta[k] >= eps:
+        if tau[k] >= eps:
             if h[k] >= eps:
-            	S = -eta[k]/h[k]
+            	S = -tau[k]/h[k]
 
             	#Update momentum
