Index: development/analytical solutions/Analytical_solution_Sampson_cross_mesh.py
===================================================================
--- development/analytical solutions/Analytical_solution_Sampson_cross_mesh.py	(revision 2648)
+++ development/analytical solutions/Analytical_solution_Sampson_cross_mesh.py	(revision 2648)
@@ -0,0 +1,168 @@
+"""Example of shallow water wave equation analytical solution
+consists of a flat water surface profile in a parabolic basin
+with linear friction. The analytical solution was derived by Sampson in 2002.
+
+   Copyright 2004
+   Christopher Zoppou, Stephen Roberts
+   ANU
+
+Specific methods pertaining to the 2D shallow water equation
+are imported from shallow_water
+for use with the generic finite volume framework
+
+Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
+numerical vector named conserved_quantities.
+"""
+
+######################
+# Module imports
+#
+
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+from pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary, gravity, linear_friction
+from math import sqrt, cos, sin, pi, exp
+from pyvolution.mesh_factory import rectangular_cross
+from pyvolution.quantity import Quantity
+from utilities.polygon import inside_polygon
+from Numeric import asarray
+from least_squares import Interpolation
+
+#-------------------------------
+# Domain
+n = 100
+m = 100
+lenx = 10000.0
+leny = 10000.0
+origin = (-5000.0, -5000.0)
+
+points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
+domain = Domain(points, elements, boundary)
+
+#----------------
+# Order of scheme
+domain.default_order = 1
+domain.smooth = True
+
+domain.quantities['linear_friction'] = Quantity(domain)
+
+#---------------------------------------------------------------------------
+# Reconstruct forcing terms with linear friction instead of manning friction
+domain.forcing_terms = []
+domain.forcing_terms.append(gravity)
+domain.forcing_terms.append(linear_friction)
+
+print domain.forcing_terms
+
+#-------------------------------------
+# Provide file name for storing output
+domain.store = True
+domain.format = 'sww'
+domain.filename = 'sampson_first_order_cross'
+
+print 'Number of triangles = ', len(domain)
+
+#-----------------------------------------
+#Reduction operation for get_vertex_values
+from utilities.numerical_tools import mean
+domain.reduction = mean
+#domain.reduction = min  #Looks better near steep slopes
+
+
+#-----------------
+#Initial condition
+#
+print 'Initial condition'
+t = 0.0
+h0 = 10.
+a = 3000.
+g = 9.81
+tau =0.001
+B = 5
+A = 0
+p = sqrt(8*g*h0/a/a)
+s = sqrt(p*p-tau*tau)/2
+t = 0.
+
+#------------------------------
+#Set bed-elevation and friction
+def x_slope(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a)
+    return z
+
+domain.set_quantity('elevation', x_slope)
+domain.set_quantity('linear_friction', tau)
+
+#-------------------
+#Set the water stage
+def stage(x,y):
+    z = x_slope(x,y)
+    n = x.shape[0]
+    h = 0*x
+    for i in range(n):
+        h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \
+               +tau*B/2*sin(s*t)))*x[i] \
+          -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i]
+    if h[i] < z[i]:
+        h[i] = z[i]
+    return h
+
+domain.set_quantity('stage', stage)
+
+
+#---------
+# Boundary
+print 'Boundary conditions'
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([0.0, 0.0, 0.0])
+domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
+
+#---------------------------------------------
+# Find triangle that contains the point points
+# and print to file
+points = [0.,0.]
+for n in range(len(domain.triangles)):
+    t = domain.triangles[n]
+    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+    if inside_polygon(points,tri):
+        print 'Point is within triangle with vertices '+'%s'%tri
+        n_point = n
+
+print 'n_point = ',n_point
+t = domain.triangles[n_point]
+tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+filename=domain.filename
+file = open(filename,'w')
+
+#----------
+# Evolution
+import time
+t0 = time.time()
+
+Stage     = domain.quantities['stage']
+Xmomentum = domain.quantities['xmomentum']
+Ymomentum = domain.quantities['ymomentum']
+
+for t in domain.evolve(yieldstep = 20.0, finaltime = 10000.0 ):
+    domain.write_time()
+
+    tri_array = asarray(tri)
+    t_array = asarray([[0,1,2]])
+    interp = Interpolation(tri_array,t_array,[points])
+
+
+    stage     = Stage.get_values(location='centroids',indices=[n_point])[0]
+    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0]
+    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0]
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
+
+file.close()
+print 'That took %.2f seconds' %(time.time()-t0)
Index: development/analytical solutions/Analytical_solution_Yoon_cross_mesh.py
===================================================================
--- development/analytical solutions/Analytical_solution_Yoon_cross_mesh.py	(revision 2648)
+++ development/analytical solutions/Analytical_solution_Yoon_cross_mesh.py	(revision 2648)
@@ -0,0 +1,169 @@
+"""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 Yoon and Chou.
+
+   Copyright 2005
+   Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia
+
+"""
+
+#---------------
+# Module imports
+#import sys
+#from os import sep
+#sys.path.append('..'+sep+'pyvolution')
+from pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary
+from math import sqrt, cos, sin, pi
+from pyvolution.mesh_factory import rectangular_cross
+from utilities.polygon import inside_polygon
+from Numeric import asarray
+from pyvolution.least_squares import Interpolation
+
+#-------------------------------
+# Domain
+n = 100
+m = 100
+delta_x = 80.0
+delta_y = 80.0
+lenx = delta_x*n
+leny = delta_y*m
+origin = (-4000.0, -4000.0)
+
+points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin)
+domain = Domain(points, elements, boundary)
+
+#----------------
+# Order of scheme
+domain.default_order = 1
+
+domain.smooth = True
+
+#-------------------------------------
+# Provide file name for storing output
+domain.store = False
+domain.format = 'sww'
+domain.filename = 'yoon_mesh_second_order_cross'
+print 'Number of triangles = ', len(domain)
+
+#----------------------------------------------------------
+# Decide which quantities are to be stored at each timestep
+domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+#------------------------------------------
+# 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)
+T = pi/omega
+
+#------------------
+# Set bed elevation
+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 initial water level
+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('stage', level)
+
+#---------
+# Boundary
+print 'Boundary conditions'
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([0.0, 0.0, 0.0])
+domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
+
+#---------------------------------------------
+# Find triangle that contains the point points
+# and print to file
+points = [0.,0.]
+for n in range(len(domain.triangles)):
+    t = domain.triangles[n]
+    tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+    if inside_polygon(points,tri):
+        print 'Point is within triangle with vertices '+'%s'%tri
+        n_point = n
+
+print 'n_point = ',n_point
+t = domain.triangles[n_point]
+tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]]
+
+filename=domain.filename
+file = open(filename,'w')
+
+#----------
+# Evolution
+import time
+t0 = time.time()
+
+time_array = []
+stage_array = []
+Stage     = domain.quantities['stage']
+Xmomentum = domain.quantities['xmomentum']
+Ymomentum = domain.quantities['ymomentum']
+
+for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ):
+    domain.write_time()
+
+    tri_array = asarray(tri)
+    t_array = asarray([[0,1,2]])
+    interp = Interpolation(tri_array,t_array,[points])
+
+
+    stage     = Stage.get_values(location='centroids',indices=[n_point])[0]
+    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0]
+    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0]
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
+
+    time_array.append(t)
+    stage_array.append(stage)
+
+file.close()
+print 'That took %.2f seconds' %(time.time()-t0)
+
+
+from pylab import *
+ion()
+hold(False)
+plot(time_array, stage_array, 'r.-')
+#title('Gauge %s' %name)
+xlabel('time(s)')
+ylabel('stage (m)')
+legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+#savefig(name, dpi = 300)
+
+#raw_input('Next')
+show()
+
+
