Index: /inundation/ga/storm_surge/analytical solutions/oblique_test_chris.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/oblique_test_chris.py	(revision 540)
+++ /inundation/ga/storm_surge/analytical solutions/oblique_test_chris.py	(revision 541)
@@ -12,11 +12,21 @@
 # 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 domain import oblique_domain, rectangular_domain
 
-from shallow_water import Friction, Constant_height
-import visualise2_chris as visualise
-import Image, ImageGrab
+from math import sqrt, cos, sin, pi
+from mesh_factory import oblique
+
 
 ######################
@@ -28,5 +38,6 @@
 lenx = 40.
 
-domain = oblique_domain(m, n, lenx, leny )
+points, elements, boundary = oblique(m, n, lenx, leny)
+domain = Domain(points, elements, boundary) 
 
 # Order of solver
@@ -51,11 +62,12 @@
     return 0*x
 
-domain.set_field_values(x_slope, Friction(0))
+domain.set_quantity('elevation', x_slope)
+domain.set_quantity('friction', 0.0)
 
 ######################
 # Boundary conditions
 #
-R = Reflective_boundary()
-T = Transmissive_boundary()
+R = Reflective_boundary(domain)
+T = Transmissive_boundary(domain)
 D = Dirichlet_boundary([1.0, 8.57, 0.0])
 
@@ -65,26 +77,19 @@
 #Initial condition
 h = 0.5
-domain.set_conserved_quantities(Constant_height(x_slope, h), None, None)
+domain.set_quantity('level', Constant_height(x_slope, h) )
+
 
 ######################
-# Visualisation using Visual Python, comment out if not required
-#visualise.create_surface(domain)
-  
-######################
 #Evolution
-#raw_input('Press key to start simulation')
 import time
 t0 = time.time()
-for t in domain.evolve(max_timestep = 0.1, yieldstep = 0.5, finaltime = 50):
+for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
     domain.write_time()
-
-# Remove the following if Visual Python not required    
-#    visualise.update(domain, index=0, smooth=True)
     
 print 'That took %.2f seconds' %(time.time()-t0)
 
-print "saving file?"
-im = ImageGrab.grab()
-im.save("ccube.eps")
+#print "saving file?"
+#im = ImageGrab.grab()
+#im.save("ccube.eps")
 
 
Index: undation/ga/storm_surge/analytical solutions/visualise2_chris.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/visualise2_chris.py	(revision 540)
+++ 	(revision )
@@ -1,255 +1,0 @@
-from visual import *
-
-def normal(v1,v2):
-    """Safe computation of normalised cross product
-    """
-    try:
-        n = norm(cross(v1, v2))
-    except ZeroDivisionError:
-        print v1, v2
-        raise Exception
-        #n = vector(0,0,0)
-
-    return n    
-
-
-class Triangle:
-    def __init__(self,v0,v1,v2,color=(1,1,1), frame=None):
-        """Create one flat triangular panel with two sides
-        """
-
-        #Outward normal
-        
-        self.normal=normal(v1-v0, v2-v0)
-
-        self.panel = faces(frame = frame)
-        self.backpanel = faces(frame = frame)        
-        for v in (v0,v1,v2):
-            self.panel.append(pos=v, normal = self.normal, color=color)
-
-        for v in (v2,v1,v0):
-            self.backpanel.append(pos=v, normal = -self.normal, color=color)            
-
-            
-    def move(self,v):
-        """Move panel in direction given by vector v
-        """
-
-        for vertex in self.panel.pos:
-            vertex += v            
-
-    def update_height(self,d,fix_baseline=False):
-        """Change height of triangle by displacing the third vertex by
-        d in the direction perpendicular to the baseline (v1-v0) and the
-        outward normal vector
-        """
-
-        v0 = self.panel.pos[0]
-        v1 = self.panel.pos[1]
-        v2 = self.panel.pos[2]        
-        
-        n = normal(v1-v0, self.normal)
-
-        if fix_baseline:
-            self.panel.pos[2] -= d*n
-        else:
-            self.panel.pos[:2] += d*n        
-
-    def update_color(self, c):
-        """Change color to c
-        """
-
-        self.panel.color = c
-        
-
-    def set_vertexheights(self, heights, floor_heights = None):
-        from Numeric import zeros, Float
-
-        from config import minimum_allowed_height as hmin
-        if floor_heights is None:
-            floor_heights = zeros(heights.shape, Float)
-
-        all_vertices_below_threshold = True    
-        for k in range(3):
-            w = heights[k]
-            z = floor_heights[k]
-
-            if w-z >= hmin:
-                all_vertices_below_threshold = False
-            
-            vertex = self.panel.pos[k]
-            vertex[2] = w
-
-            vertex = self.backpanel.pos[2-k]
-            vertex[2] = w
-
-
-        #Update color to visualise dry areas 
-        if all_vertices_below_threshold: 
-            self.panel.color = self.bottom_color
-            self.backpanel.color = self.bottom_color                
-        else:
-            self.panel.color = self.top_color
-            self.backpanel.color = self.top_color
-
-
-        #update normal    
-        v0 = self.panel.pos[0]
-        v1 = self.panel.pos[1]
-        v2 = self.panel.pos[2]
-
-
-        n = normal(v1-v0, v2-v0)
-        self.panel.normal=n
-        self.backpanel.normal=-n
-
-        
-
-
-def create_surface(domain, index=0, smooth=True):
-    """Link surface of domains to their visual representations
-    """
-
-    fr = frame() #Default frame to contain all objects
-    s=0
-
-    for volume in domain:
-
-        v0, v1, v2 = volume.get_vertex_objects()        
-
-        z0 = volume.get_conserved_quantities('vertex', 0)[index]
-        z1 = volume.get_conserved_quantities('vertex', 1)[index]
-        z2 = volume.get_conserved_quantities('vertex', 2)[index]
-        
-        V0 = vector(v0.get_x(), v0.get_y(), z0)
-        V1 = vector(v1.get_x(), v1.get_y(), z1)
-        V2 = vector(v2.get_x(), v2.get_y(), z2)
-
-
-        #Top surface
-        c0 = 0.1
-        c1 = 0.4
-        c2 = 0.99
-        if s:
-            c2 = 0.7*c2   #To show triangles in slightly different shades
-
-        s = 1-s
-
-        col = (c0,c1,c2)
-            
-        volume.visual_top = Triangle(V0,V1,V2,color=col,frame=fr)
-        volume.visual_top.top_color = col
-
-    
-        #Bottom surface
-        #v0, v1, v2 = volume.vertices
-        #v0 = vector(v0.x, v0.y, v0.field_values[index])
-        #v1 = vector(v1.x, v1.y, v1.field_values[index])
-        #v2 = vector(v2.x, v2.y, v2.field_values[index])
-        z0 = volume.get_field_values('vertex', 0)[0]
-        z1 = volume.get_field_values('vertex', 1)[0]
-        z2 = volume.get_field_values('vertex', 2)[0]
-        
-        V0 = vector(v0.get_x(), v0.get_y(), z0)
-        V1 = vector(v1.get_x(), v1.get_y(), z1)
-        V2 = vector(v2.get_x(), v2.get_y(), z2)        
-
-        c0 = 0.3
-        c1 = 0.3
-        c2 = 0.3
-        if s:
-            c2 = 0.4*c2   #To show triangles in slightly different shades
-
-        col = (c0,c1,c2)   
-        volume.visual_bottom = Triangle(V0, V1, V2, color=col,frame=fr)
-
-        volume.visual_top.bottom_color=col
-        
-    update(domain, index, smooth)
-    #print 'Scale', scene.scale    
-
-
-def update(domain, index = 0, smooth = True):
-    """Update vertex heights.
-    The argument index refers to which conserved quantity to visualise.    
-    If smooth is set True, vertex values will be averaged
-    yielding a smoother surface.
-    """
-
-    from Numeric import array
-
-    #print scene.forward
-    if smooth:
-        #Get all average point values
-        vertex_heights = {}
-        for k in range(len(domain)):
-            vol = domain[k]
-
-            for i in range(len(vol.vertices)):
-                vertex = vol.vertices[i]
-                if vertex_heights.has_key(vertex):
-                    vertex_heights[vertex].append(
-                        vol.get_conserved_quantities('vertex',i)[index])
-                else:
-                    vertex_heights[vertex] = []
-                    vertex_heights[vertex].append(
-                        vol.get_conserved_quantities('vertex',i)[index])
-
-    
-    for volume in domain:
-        if smooth:
-            #The averages
-            x = zeros(3, Float)
-            i=0
-            for vertex in volume.vertices:
-                A = array(vertex_heights[vertex])
-                x[i] = sum(A)/len(A)
-                i+=1
-        else:        
-            #The true values
-            x = [volume.get_conserved_quantities('vertex',0)[index],
-                 volume.get_conserved_quantities('vertex',1)[index],
-                 volume.get_conserved_quantities('vertex',2)[index]]
-
-            
-        #Do it
-        floor_heights = array([pos[2] for pos in volume.visual_bottom.panel.pos])
-        
-        volume.visual_top.set_vertexheights(x, floor_heights)
-
-
-
-
-
-#scene.width = 1000
-#scene.height = 800
-
-#Original
-#scene.center = (0.5,0.5,0)
-#scene.forward = vector(-0.1, 0.5, -0.5)
-
-#Temporary (for bedslope)
-#scene.forward = vector(0.0006, 0.7, -0.03)    
-
-
-#Temporary for hackett - begin
-#scene.autoscale = 0
-#scene.scale = (0.002, 0.002, 0.01) #Scale z so that countours stand out more
-#scene.center = (300.0,500.0,-10)
-#Temporary for hackett - end
-
-    
-scene.ambient = 0.4
-scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4),
-                (-0.2, 0.2, 0.1)]
-
-scene2 = display(fullscreen=1,center=(0,0,0),
-                 background=(1,1,1),userspin = 1,
-                 uniform = 1,autoscale = 1,ambient = 0.5,
-                 lights = [(0., 0.9, 0.3)])
-
-
-
-
-
-
Index: /inundation/ga/storm_surge/pyvolution/mesh_factory.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 540)
+++ /inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 541)
@@ -63,73 +63,73 @@
 
 
-# def oblique_mesh(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0),
-#                   point_class=Point, element_class=Triangle, mesh_class=Mesh):
-#     """Setup a oblique grid of triangles
-#     with m segments in the x-direction and n segments in the y-direction
-
-#     Triangle refers to the actual class or subclass to be instantiated:
-#     e.g. if Volume is a subclass of Triangle,
-#     this function can be invoked with the keywords
-#       oblique_mesh(...,Triangle=Volume, Mesh=Domain)
-#     """
-
-#     from Numeric import array
-#     from visual import rate
-#     import math
-
-#     from config import epsilon
-
-#     E = m*n*2        #Number of triangular elements
-#     P = (m+1)*(n+1)  #Number of initial vertices
-
-#     initialise_consecutive_datastructure(P+E, E,
-#                                          point_class,
-#                                          element_class,
-#                                          mesh_class)
-
-#     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 = {}
-
-#     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
+def oblique(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] = Point(x + origin[0],y + origin[1])
-
-#     # Construct 2 elements per element        
-#     elements = []
-#     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]
- 
-#             elements.append(element_class(v4,v3,v2)) #Lower                
-#             elements.append(element_class(v1,v2,v3)) #Upper    
-
-#     M = mesh_class(elements)
-
-#     #Set a default tagging
-
+            vertices[i,j] = len(points)
+            points.append([x + origin[0], y + origin[1]])
+
+    # Construct 2 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]
+            
+            #Update boundary dictionary and create elements
+            if i == m-1:
+                boundary[(len(elements), 2)] = 'right'
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'                
+            elements.append([v4,v3,v2]) #Lower
+
+            if i == 0:
+                boundary[(len(elements), 2)] = 'left'
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top' 
+             
+            elements.append([v1,v2,v3]) #Upper
+
+    return points, elements, boundary       
+
+#OLD FORMULA FOR SETTING BOUNDARY TAGS - OBSOLETE NOW             
 #     for id, face in M.boundary:
 
@@ -159,6 +159,8 @@
 #             print face, id, id%m, m, n
 #             raise 'Unknown boundary'          
-    
+#   
 #     return M
+
+
 
 def from_polyfile(name):
