Index: /inundation/analytical solutions/Analytical solution_oblique_shock_01.py
===================================================================
--- /inundation/analytical solutions/Analytical solution_oblique_shock_01.py	(revision 1909)
+++ /inundation/analytical solutions/Analytical solution_oblique_shock_01.py	(revision 1909)
@@ -0,0 +1,116 @@
+"""
+Example of shallow water wave equation
+consisting of an asymetrical converging channel.
+
+Copyright 2004
+Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
+Geoscience Australia, 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
+#
+
+#Were these used?
+#import visualise2_chris as visualise
+#import Image, ImageGrab
+
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+
+from shallow_water import Domain, Constant_height
+from shallow_water import Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary
+
+from math import sqrt, cos, sin, pi
+from mesh_factory import oblique_cross
+
+
+######################
+# Domain
+#
+leny = 30.
+lenx = 40.
+n = 50
+m = 60
+n=n*4
+m=m*4
+
+points, elements, boundary = oblique_cross(m, n, lenx, leny)
+domain = Domain(points, elements, boundary)
+
+# Order of solver
+domain.default_order=2
+
+# Store output
+domain.store=True
+
+# Output format
+domain.format="sww" #NET.CDF binary format
+                    # "dat" for ASCII
+
+# Provide file name for storing output
+domain.filename="oblique"
+
+# Visualization smoothing
+domain.smooth=True
+domain.visualise=True
+
+#######################
+#Bed-slope and friction
+def x_slope(x, y):
+    return 0*x
+
+domain.set_quantity('elevation', x_slope)
+domain.set_quantity('friction', 0.0)
+
+######################
+# Boundary conditions
+#
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary([1.82, 8.49, 0.0])
+
+
+domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R})
+
+######################
+#Initial condition
+h = 0.5
+domain.set_quantity('stage', Constant_height(x_slope, h) )
+
+
+
+#for id, face in domain.boundary:
+#    print '(',id,',',face,') ',domain.boundary[(id,face)]
+
+######################
+#Evolution
+import time
+t0 = time.time()
+for t in domain.evolve(yieldstep = 0.1, finaltime = 5):
+    domain.write_time()
+    print domain.quantities['stage'].get_values(location='centroids',indices=[3399])
+    print domain.quantities['xmomentum'].get_values(location='centroids',indices=[3399])
+    print domain.quantities['ymomentum'].get_values(location='centroids',indices=[3399])
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
+#FIXME: Compute average water depth on either side of shock and compare
+#to expected values. And also Froude numbers.
+
+
+#print "saving file?"
+#im = ImageGrab.grab()
+#im.save("ccube.eps")
+
Index: /inundation/pyvolution/mesh_factory.py
===================================================================
--- /inundation/pyvolution/mesh_factory.py	(revision 1908)
+++ /inundation/pyvolution/mesh_factory.py	(revision 1909)
@@ -676,2 +676,85 @@
 
     return points, elements, boundary
+
+
+
+
+def oblique_cross(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
+    """Setup a oblique grid of triangles
+    with m segments in the x-direction and n segments in the y-direction
+
+    """
+
+    from Numeric import array
+    import math
+
+    from config import epsilon
+
+
+    deltax = lenx/float(m)
+    deltay = leny/float(n)
+    a = 0.75*lenx*math.tan(theta/180.*math.pi)
+    x1 = lenx
+    y1 = 0
+    x2 = lenx
+    y2 = leny
+    x3 = 0.25*lenx
+    y3 = leny
+    x4 = x3
+    y4 = 0
+    a2 = a/(x1-x4)
+    a1 = -a2*x4
+    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
+    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3
+
+    # Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        x = deltax*i
+        for j in range(n+1):
+            y = deltay*j
+            if x > 0.25*lenx:
+                y = a1 + a2*x + a3*y + a4*x*y
+
+            vertices[i,j] = len(points)
+            points.append([x + origin[0], y + origin[1]])
+
+    # Construct 4 triangles per element
+    elements = []
+    boundary = {}
+    for i in range(m):
+        for j in range(n):
+            v1 = vertices[i,j+1]
+            v2 = vertices[i,j]
+            v3 = vertices[i+1,j+1]
+            v4 = vertices[i+1,j]
+            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
+            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
+            v5 = len(points)
+            points.append([x, y])
+
+            #Update boundary dictionary and create elements
+                        #Create left triangle
+            if i == 0:
+                boundary[(len(elements), 1)] = 'left'
+            elements.append([v2,v5,v1])
+
+            #Create bottom triangle
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v5,v2])
+
+            #Create right triangle
+            if i == m-1:
+                boundary[(len(elements), 1)] = 'right'
+            elements.append([v3,v5,v4])
+
+            #Create top triangle
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v5,v3])
+
+
+    return points, elements, boundary
