Index: /anuga_core/source/anuga/pmesh/mesh.py
===================================================================
--- /anuga_core/source/anuga/pmesh/mesh.py	(revision 4905)
+++ /anuga_core/source/anuga/pmesh/mesh.py	(revision 4906)
@@ -53,6 +53,4 @@
     kinds = _kinds()
     
-SET_COLOUR='red'
-
 #FIXME: this is not tested.
 from anuga.utilities.numerical_tools import gradient
@@ -299,141 +297,4 @@
                                    self.getTag())
         
-class Triangle(MeshObject):
-    """
-    A triangle element, defined by 3 vertices.
-    Attributes based on the Triangle program.
-    """
-
-    def __init__(self, vertex1, vertex2, vertex3, attribute = None,
-                 neighbors = None ):
-        """
-        Vertices, the initial arguments, are listed in counterclockwise order.
-        """
-        self.vertices= [vertex1,vertex2, vertex3 ]
-        
-        if attribute is None:
-            self.attribute =""
-        else:
-            self.attribute = attribute #this is a string
-            
-        if neighbors is None:
-            self.neighbors=[]
-        else:
-            self.neighbors=neighbors
-
-    def replace(self,new_triangle):
-        self = new_triangle
-
-    def longestSideID(self):
-        ax = self.vertices[0].x
-        ay = self.vertices[0].y
-        
-        bx = self.vertices[1].x
-        by = self.vertices[1].y
-        
-        cx = self.vertices[2].x
-        cy = self.vertices[2].y
-
-        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
-        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
-        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
- 
-        len = [lenA,lenB,lenC]
-        return len.index(max(len))
-
-    def rotate(self,offset):
-        """
-        permute the order of the sides of the triangle
-        offset must be 0,1 or 2
-        """
-
-        if offset == 0:
-            pass
-        else:
-            if offset == 1:
-                self.vertices = [self.vertices[1],self.vertices[2],
-                                 self.vertices[0]]
-                self.neighbors = [self.neighbors[1],self.neighbors[2],
-                                  self.neighbors[0]]
-            if offset == 2:
-                self.vertices = [self.vertices[2],self.vertices[0],
-                                 self.vertices[1]]
-                self.neighbors = [self.neighbors[2],self.neighbors[0],
-                                  self.neighbors[1]]
-
-    def rotate_longest_side(self):
-        self.rotate(self.longestSideID())
-
-    def getVertices(self):
-        return self.vertices
-    
-    def get_vertices(self):
-        """
-        Return a list of the vertices.  The x and y values will be relative
-        Easting and Northings for the zone of the current geo_ref.
-        """
-        return self.vertices
-    
-    def calcArea(self):
-        ax = self.vertices[0].x
-        ay = self.vertices[0].y
-        
-        bx = self.vertices[1].x
-        by = self.vertices[1].y
-        
-        cx = self.vertices[2].x
-        cy = self.vertices[2].y
-        
-        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
-    
-    def calcP(self):
-        #calculate the perimeter
-        ax = self.vertices[0].x
-        ay = self.vertices[0].y
-        
-        bx = self.vertices[1].x
-        by = self.vertices[1].y
-        
-        cx = self.vertices[2].x
-        cy = self.vertices[2].y
-
-        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
-        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
-        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
-
-        return a+b+c
-            
-    def setNeighbors(self,neighbor1=None, neighbor2=None, neighbor3=None):
-        """
-        neighbor1 is the triangle opposite vertex1 and so on.
-        Null represents no neighbor
-        """
-        self.neighbors = [neighbor1, neighbor2, neighbor3]
-
-    def setAttribute(self,attribute):
-        """
-        neighbor1 is the triangle opposite vertex1 and so on.
-        Null represents no neighbor
-        """
-        self.attribute = attribute #this is a string
-        
-    def __repr__(self):
-        return "[%s,%s]" % (self.vertices,self.attribute)
-        
-
-    def draw(self, canvas, tags, scale=1, xoffset=0, yoffset=0,
-             colour="green"):
-        """
-        Draw a triangle, returning the objectID
-        """
-        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
-                                     scale*-1*(self.vertices[1].y + yoffset),
-                                     scale*(self.vertices[0].x + xoffset),
-                                     scale*-1*(self.vertices[0].y + yoffset),
-                                     scale*(self.vertices[2].x + xoffset),
-                                     scale*-1*(self.vertices[2].y + yoffset),
-                                     tags = tags,
-                                     outline = colour,fill = '')
-
 class Segment(MeshObject):
     """
@@ -631,17 +492,5 @@
         self.tri_mesh=None
         
-        self.setID={}
-        #a dictionary of names.
-        #multiple sets are allowed, but the gui does not yet
-        #support this
-        
-        self.setID['None']=0
-        #contains the names of the sets pointing to the indexes
-        #in the list. 
-        
-	self.sets=[[]]
-        #Contains the lists of triangles (triangle sets)
-
-       
+        
         self.visualise_graph = True
 
@@ -2338,554 +2187,4 @@
             self.regions.append(Object)
   
-############################################
-
-        
-    def refineSet(self,setName):
-        Triangles = self.sets[self.setID[setName]]
-        Refine(self,Triangles)
-
-    def selectAllTriangles(self):
-        A=[]
-        A.extend(self.meshTriangles)
-        if not('All' in self.setID.keys()):
-            self.setID['All']=len(self.sets)
-            self.sets.append(A)
-        else:
-            self.sets[self.setID['All']]=A
-        return 'All'
-        # and objectIDs
-
-
-    def clearSelection(self):
-        A = []
-        if not('None' in self.setID.keys()):
-            self.setID['None']=len(self.sets)
-            self.sets.append(A)
-        return 'None'
-
-    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
-    #FIXME Draws over previous triangles - may bloat canvas
-        Triangles = self.sets[self.setID[setName]]
-        for triangle in Triangles:
-            triangle.draw(canvas,1,
-                          scale = SCALE,
-                          colour = colour)
-	    
-    def undrawSet(self,canvas,setName,SCALE,colour='green'):
-    #FIXME Draws over previous lines - may bloat canvas
-        Triangles = self.sets[self.setID[setName]]
-        for triangle in Triangles:
-            triangle.draw(canvas,1,
-                          scale = SCALE,
-                          colour = colour)
-
-    def weed(self,Vertices,Segments):
-        #Depreciated
-        #weed out existing duplicates
-        print 'len(self.getUserSegments())'
-        print len(self.getUserSegments())
-        print 'len(self.getUserVertices())'
-        print len(self.getUserVertices())
-
-        point_keys = {}
-        for vertex in Vertices:
-            point = (vertex.x,vertex.y)
-            point_keys[point]=vertex
-        #inlined would looks very ugly
-
-        line_keys = {}
-        for segment in Segments:
-            vertex1 = segment.vertices[0]
-            vertex2 = segment.vertices[1]
-            point1 = (vertex1.x,vertex1.y)
-            point2 = (vertex2.x,vertex2.y)
-            segment.vertices[0]=point_keys[point1]
-            segment.vertices[1]=point_keys[point2]
-            vertex1 = segment.vertices[0]
-            vertex2 = segment.vertices[1]
-            point1 = (vertex1.x,vertex1.y)
-            point2 = (vertex2.x,vertex2.y)
-            line1 = (point1,point2)
-            line2 = (point2,point1)
-            if not (line_keys.has_key(line1) \
-                 or line_keys.has_key(line2)):
-                 line_keys[line1]=segment
-        Vertices=point_keys.values()
-        Segments=line_keys.values()
-        return Vertices,Segments
-
-    def segs_to_dict(self,segments):
-        dict={}
-        for segment in segments:
-            vertex1 = segment.vertices[0]
-            vertex2 = segment.vertices[1]
-            point1 = (vertex1.x,vertex1.y)
-            point2 = (vertex2.x,vertex2.y)
-            line = (point1,point2)
-            dict[line]=segment
-        return dict
-
-    def seg2line(self,s):
-        return ((s.vertices[0].x,s.vertices[0].y,)\
-                (s.vertices[1].x,s.vertices[1].y))
-
-    def line2seg(self,line,tag=None):
-        point0 = self.point2ver(line[0])
-        point1 = self.point2ver(line[1])
-        return Segment(point0,point1,tag=tag)
-
-    def ver2point(self,vertex):
-        return (vertex.x,vertex.y)
-
-    def point2ver(self,point):
-        return Vertex(point[0],point[1])
-
-    def smooth_polySet(self,min_radius=0.05):
-        #for all pairs of connecting segments:
-        #    propose a new segment that replaces the 2
-
-        #    If the difference between the new segment
-        #    and the old lines is small: replace the
-        #    old lines.
-
-        seg2line = self.seg2line
-        ver2point= self.ver2point
-        line2seg = self.line2seg
-        point2ver= self.point2ver
-
-        #create dictionaries of lines -> segments
-        userSegments = self.segs_to_dict(self.userSegments)
-        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
-
-        #lump user and alpha segments
-        for key in alphaSegments.keys():
-            userSegments[key]=alphaSegments[key]
-
-        #point_keys = tuple -> vertex
-        #userVertices = vertex -> [line,line] - lines from that node
-        point_keys = {}
-        userVertices={}
-        for vertex in self.getUserVertices():
-            point = ver2point(vertex)
-            if not point_keys.has_key(point):
-                point_keys[point]=vertex
-                userVertices[vertex]=[]
-        for key in userSegments.keys():
-            line = key
-            point_0 = key[0]
-            point_1 = key[1]
-            userVertices[point_keys[point_0]].append(line)
-            userVertices[point_keys[point_1]].append(line)
-
-        for point in point_keys.keys():
-            try:
-            #removed keys can cause keyerrors
-                vertex = point_keys[point]
-		lines = userVertices[vertex]
-    
-                #if there are 2 lines on the node
-		if len(lines)==2:
-		    line_0 = lines[0]
-		    line_1 = lines[1]
-    
-                    #if the tags are the the same on the 2 lines
-		    if userSegments[line_0].tag == userSegments[line_1].tag:
-			tag = userSegments[line_0].tag 
-    
-                        #point_a is one of the next nodes, point_b is the other
-			if point==line_0[0]:
-			    point_a = line_0[1]
-			if point==line_0[1]:
-			    point_a = line_0[0]
-			if point==line_1[0]:
-			    point_b = line_1[1]
-			if point==line_1[1]:
-			    point_b = line_1[0]
-    
-    
-                        #line_2 is proposed
-			line_2 = (point_a,point_b)
-
-			#calculate the area of the triangle between
-			#the two existing segments and the proposed
-			#new segment
-			ax = point_a[0]
-			ay = point_a[1]
-			bx = point_b[0]
-			by = point_b[1]
-			cx = point[0]
-			cy = point[1]
-			area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
-
-			#calculate the perimeter
-			len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
-			len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
-			len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
-			perimeter = len_a+len_b+len_c
-
-			#calculate the radius
-			r = area/(2*perimeter)
-
-			#if the radius is small: then replace the existing
-			#segments with the new one
-			if r < min_radius:
-			    if len_c < min_radius: append = False
-			    else: append = True
-			    #if the new seg is also time, don't add it
-			    if append:
-				segment = self.line2seg(line_2,tag=tag)
-
-			    list_a=userVertices[point_keys[point_a]]
-			    list_b=userVertices[point_keys[point_b]]
-
-			    if line_0 in list_a:
-				list_a.remove(line_0)
-			    else:
-				list_a.remove(line_1)
-
-			    if line_0 in list_b:
-				list_b.remove(line_0)
-			    else:
-				list_b.remove(line_1)
-
-			    if append:
-				list_a.append(line_2)
-				list_b.append(line_2)
-			    else:
-				if len(list_a)==0:
-				    userVertices.pop(point_keys[point_a])
-				    point_keys.pop(point_a)
-				if len(list_b)==0:
-				    userVertices.pop(point_keys[point_b])
-				    point_keys.pop(point_b)
-
-			    userVertices.pop(point_keys[point])
-			    point_keys.pop(point)
-			    userSegments.pop(line_0)
-			    userSegments.pop(line_1)
-
-			    if append:
-				userSegments[line_2]=segment
-            except:
-                pass
-
-        #self.userVerticies = userVertices.keys()
-        #self.userSegments = []
-        #for key in userSegments.keys():
-        #    self.userSegments.append(userSegments[key])
-        #self.alphaUserSegments = []
-
-        self.userVerticies = []
-        self.userSegments = []
-        self.alphaUserSegments = []
-
-        return userVertices,userSegments,alphaSegments
-
-    def triangles_to_polySet(self,setName):
-        #self.smooth_polySet()
-
-        seg2line = self.seg2line
-        ver2point= self.ver2point
-        line2seg = self.line2seg
-        point2ver= self.point2ver
-
-        from Numeric import array,allclose
-        #turn the triangles into a set
-        Triangles = self.sets[self.setID[setName]]
-        Triangles_dict = {}
-        for triangle in Triangles:
-            Triangles_dict[triangle]=None
- 
-
-        #create a dict of points to vertexes (tuple -> object)
-        #also create a set of vertexes (object -> True)
-        point_keys = {}
-        userVertices={}
-        for vertex in self.getUserVertices():
-            point = ver2point(vertex)
-            if not point_keys.has_key(point):
-                point_keys[point]=vertex
-                userVertices[vertex]=True
-
-        #create a dict of lines to segments (tuple -> object)
-        userSegments = self.segs_to_dict(self.userSegments)
-        #append the userlines in an affine linespace
-        affine_lines = Affine_Linespace()
-        for line in userSegments.keys():
-            affine_lines.append(line)
-        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
-        for line in alphaSegments.keys():
-            affine_lines.append(line)
-        
-        for triangle in Triangles:
-            for i in (0,1,2):
-                #for every triangles neighbour:
-                if not Triangles_dict.has_key(triangle.neighbors[i]):
-                #if the neighbour is not in the set:
-                    a = triangle.vertices[i-1]
-                    b = triangle.vertices[i-2]
-                    #Get possible matches:
-                    point_a = ver2point(a)
-                    point_b = ver2point(b)
-                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
-                    line = (point_a,point_b)
-                    tag = None
-
-
-                    #this bit checks for matching lines
-                    possible_lines = affine_lines[line] 
-                    possible_lines = unique(possible_lines)
-                    found = 0                            
-                    for user_line in possible_lines:
-                        if self.point_on_line(midpoint,user_line):
-                            found+=1
-                            assert found<2
-                            if userSegments.has_key(user_line):
-                                parent_segment = userSegments.pop(user_line)
-                            if alphaSegments.has_key(user_line):
-                                parent_segment = alphaSegments.pop(user_line)
-                            tag = parent_segment.tag
-                            offspring = [line]
-                            offspring.extend(self.subtract_line(user_line,
-                                                                line))
-                            affine_lines.remove(user_line)
-                            for newline in offspring:
-                                line_vertices = []
-                                for point in newline:
-                                    if point_keys.has_key(point):
-                                        vert = point_keys[point]
-                                    else:
-                                        vert = Vertex(point[0],point[1])
-                                        userVertices[vert]=True
-                                        point_keys[point]=vert
-                                    line_vertices.append(vert)
-                                segment = Segment(line_vertices[0],
-                                                  line_vertices[1],tag)
-                                userSegments[newline]=segment
-                                affine_lines.append(newline)
-                            #break
-                    assert found<2
-
-
-
-                    #if no matching lines
-                    if not found:
-                        line_vertices = []
-		        for point in line:
-			    if point_keys.has_key(point):
-			        vert = point_keys[point]
-			    else:
-				vert = Vertex(point[0],point[1])
-				userVertices[vert]=True
-				point_keys[point]=vert
-			    line_vertices.append(vert)
-			segment = Segment(line_vertices[0],
-                                          line_vertices[1],tag)
-                        userSegments[line]=segment
-                        affine_lines.append(line)
-        
-        self.userVerticies = []
-        self.userSegments = []
-        self.alphaUserSegments = []
-
-        return userVertices,userSegments,alphaSegments
-
-    def subtract_line(self,parent,child):
-        #Subtracts child from parent
-        #Requires that the child is a 
-        #subline of parent to work.
-
-        from Numeric import allclose,dot,array
-        A= parent[0]
-        B= parent[1]
-        a = child[0]
-        b = child[1]
-
-        A_array = array(parent[0])
-        B_array = array(parent[1])
-        a_array   = array(child[0])
-        b_array   = array(child[1])
-
-        assert not A == B
-        assert not a == b
-
-        answer = []
-
-        #if the new line does not share a 
-        #vertex with the old one
-        if not (allclose(A_array,a_array)\
-             or allclose(B_array,b_array)\
-             or allclose(A_array,b_array)\
-             or allclose(a_array,B_array)):
-            if dot(A_array-a_array,A_array-a_array) \
-            < dot(A_array-b_array,A_array-b_array):
-                sibling1 = (A,a)
-                sibling2 = (B,b)
-                return [sibling1,sibling2]
-            else:
-                sibling1 = (A,b)
-                sibling2 = (B,a)
-                return [sibling1,sibling2]
-
-        elif allclose(A_array,a_array):
-            if allclose(B_array,b_array):
-                return []
-	    else:
-                sibling = (b,B)
-                return [sibling]
-    	elif allclose(B_array,b_array):
-            sibling = (a,A)
-            return [sibling]
-
-        elif allclose(A_array,b_array):
-            if allclose(B,a):
-                return []
-	    else:
-                sibling = (a,B)
-                return [sibling]
-        elif allclose(a_array,B_array):
-            sibling = (b,A)
-            return [sibling]
-
-    def point_on_line(self,point,line):
-        #returns true within a tolerance of 3 degrees
-        x=point[0]
-        y=point[1]
-        x0=line[0][0]
-        x1=line[1][0]
-        y0=line[0][1]
-        y1=line[1][1]
-        from Numeric import array, dot, allclose
-        from math import sqrt
-        tol = 3. #DEGREES
-        tol = tol*3.1415/180
-
-        a = array([x - x0, y - y0]) 
-        a_normal = array([a[1], -a[0]])      
-        len_a_normal = sqrt(sum(a_normal**2)) 
-
-        b = array([x1 - x0, y1 - y0])  	       
-        len_b = sqrt(sum(b**2)) 
-    
-        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
-            #Point is somewhere on the infinite extension of the line
-
-            len_a = sqrt(sum(a**2))     	       	                
-            if dot(a, b) >= 0 and len_a <= len_b:
-               return True
-            else:   
-               return False
-        else:
-          return False
-
-    def line_length(self,line):
-        x0=line[0][0]
-        x1=line[1][0]
-        y0=line[0][1]
-        y1=line[1][1]
-        return ((x1-x0)**2-(y1-y0)**2)**0.5     
-
-    def threshold(self,setName,min=None,max=None,attribute_name='elevation'):
-        """
-        threshold using  d
-        """
-        triangles = self.sets[self.setID[setName]]
-        A = []
-
-        if attribute_name in self.attributeTitles:
-            i = self.attributeTitles.index(attribute_name)
-        else: i = -1#no attribute
-        if not max == None:
-            for t in triangles:
-                if (min<self.av_att(t,i)<max):
-                    A.append(t)
-        else:
-            for t in triangles:
-                if (min<self.av_att(t,i)):
-                    A.append(t)
-        self.sets[self.setID[setName]] = A
-
-    def general_threshold(self,setName,min=None,max=None\
-              ,attribute_name = 'elevation',function=None):
-        """
-        Thresholds the triangles 
-        """
-        from visual.graph import arange,ghistogram,color as colour
-        triangles = self.sets[self.setID[setName]]
-        A = []
-        data=[]
-        #data is for the graph
-
-        if attribute_name in self.attributeTitles:
-            i = self.attributeTitles.index(attribute_name)
-        else: i = -1
-        if not max == None:
-            for t in triangles:
-                value=function(t,i)
-                if (min<value<max):
-                    A.append(t)
-                data.append(value)
-        else:
-            for t in triangles:
-                value=function(t,i)
-                if (min<value):
-                    A.append(t)
-                data.append(value)
-        self.sets[self.setID[setName]] = A
-
-        if self.visualise_graph:
-            if len(data)>0:
-                max=data[0]
-                min=data[0]
-                for value in data:
-                    if value > max:
-                        max = value
-                    if value < min:
-                        min = value
-
-                inc = (max-min)/100
-
-                histogram = ghistogram(bins=arange(min,max,inc),\
-                             color = colour.red)
-                histogram.plot(data=data)
-        
-    def av_att(self,triangle,i):
-        if i==-1: return 1
-        else:
-            #evaluates the average attribute of the vertices of a triangle.
-            V = triangle.getVertices()
-            a0 = (V[0].attributes[i])
-            a1 = (V[1].attributes[i])
-            a2 = (V[2].attributes[i])
-            return (a0+a1+a2)/3
-
-    def Courant_ratio(self,triangle,index):
-        """
-        Uses the courant threshold
-        """
-        e = self.av_att(triangle,index)
-        A = triangle.calcArea()
-        P = triangle.calcP()
-        r = A/(2*P)
-        e = max(0.1,abs(e))
-        return r/e**0.5
-
-    def Gradient(self,triangle,index):
-        V = triangle.vertices
-        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
-        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
-        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
-            print ((grad_x**2)+(grad_y**2))**(0.5)
-        return ((grad_x**2)+(grad_y**2))**(0.5)
-    
-
-    def append_triangle(self,triangle):
-        self.meshTriangles.append(triangle)
-
-    def replace_triangle(self,triangle,replacement):
-        i = self.meshTriangles.index(triangle)
-        self.meshTriangles[i]=replacement
-        assert replacement in self.meshTriangles
-
 def importUngenerateFile(ofile):
     """
@@ -3210,628 +2509,4 @@
     return u
 
-"""Refines triangles
-
-   Implements the #triangular bisection?# algorithm.
- 
-   
-"""
-
-def Refine(mesh, triangles):
-    """
-    Given a general_mesh, and a triangle number, split
-    that triangle in the mesh in half. Then to prevent
-    vertices and edges from meeting, keep refining 
-    neighbouring triangles until the mesh is clean.
-    """
-    state = BisectionState(mesh)
-    for triangle in triangles:
-        if not state.refined_triangles.has_key(triangle):
-            triangle.rotate_longest_side()
-            state.start(triangle)
-            Refine_mesh(mesh, state)
-
-def Refine_mesh(mesh, state):
-    """
-    """
-    state.getState(mesh)
-    refine_triangle(mesh,state)
-    state.evolve()
-    if not state.end:
-        Refine_mesh(mesh,state)
-
-def refine_triangle(mesh,state):
-    split(mesh,state.current_triangle,state.new_point)
-    if state.case == 'one':
-        state.r[3]=state.current_triangle#triangle 2
-
-        new_triangle_id = len(mesh.meshTriangles)-1
-        new_triangle = mesh.meshTriangles[new_triangle_id]
-
-        split(mesh,new_triangle,state.old_point)
-        state.r[2]=new_triangle#triangle 1.2
-        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
-        r = state.r
-        state.repairCaseOne()
-
-    if state.case == 'two':
-        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
-
-        new_triangle = state.current_triangle
-
-        split(mesh,new_triangle,state.old_point)
-
-        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
-        state.r[4]=new_triangle#triangle 2.2
-        r = state.r
-
-        state.repairCaseTwo()
-
-    if state.case == 'vertex':
-        state.r[2]=state.current_triangle#triangle 2
-        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
-        r = state.r
-        state.repairCaseVertex()
-        
-    if state.case == 'start':
-        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
-        state.r[3]=state.current_triangle#triangle 2
-
-    if state.next_case == 'boundary':
-        state.repairCaseBoundary()
-
-
-def split(mesh, triangle, new_point):
-	"""
-	Given a mesh, triangle_id and a new point,
-	split the corrosponding triangle into two
-	new triangles and update the mesh.
-	"""
-
-	new_triangle1 = Triangle(new_point,triangle.vertices[0],
-                                 triangle.vertices[1],
-                                 attribute = triangle.attribute,
-                                 neighbors = None)
-	new_triangle2 = Triangle(new_point,triangle.vertices[2],
-                                 triangle.vertices[0],
-                                 attribute = triangle.attribute,
-                                 neighbors = None)
-
-        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
-        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
-
-        mesh.meshTriangles.append(new_triangle1)
-
-        triangle.vertices = new_triangle2.vertices
-        triangle.neighbors = new_triangle2.neighbors
-
-
-class State:
-
-    def __init__(self):
-        pass
-
-class BisectionState(State):
-
-
-    def __init__(self,mesh):
-        self.len = len(mesh.meshTriangles)
-        self.refined_triangles = {}
-        self.mesh = mesh
-        self.current_triangle = None
-        self.case = 'start'
-        self.end = False
-        self.r = [None,None,None,None,None]
-
-    def start(self, triangle):
-        self.current_triangle = triangle
-        self.case = 'start'
-        self.end = False
-        self.r = [None,None,None,None,None]
-
-    def getState(self,mesh):
-        if not self.case == 'vertex':
-            self.new_point=self.getNewVertex(mesh, self.current_triangle)
-            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
-            self.neighbour = self.current_triangle.neighbors[0]
-            if not self.neighbour is None:
-                self.neighbour.rotate_longest_side()
-            self.next_case = self.get_next_case(mesh,self.neighbour,
-                                                self.current_triangle)
-        if self.case == 'vertex':
-            self.new_point=self.old_point
-
-
-    def evolve(self):
-        if self.case == 'vertex':
-            self.end = True
-
-        self.last_case = self.case
-        self.case = self.next_case
-
-        self.old_point = self.new_point
-        self.current_triangle = self.neighbour
-
-        if self.case == 'boundary':
-            self.end = True
-        self.refined_triangles[self.r[2]]=1
-        self.refined_triangles[self.r[3]]=1
-        if not self.r[4] is None:
-            self.refined_triangles[self.r[4]]=1
-        self.r[0]=self.r[2]
-        self.r[1]=self.r[3]
-
-
-    def getNewVertex(self,mesh,triangle):
-	coordinate1 = triangle.vertices[1]
-	coordinate2 = triangle.vertices[2]
-	a = ([coordinate1.x*1.,coordinate1.y*1.])
-	b = ([coordinate2.x*1.,coordinate2.y*1.])
-        attributes = []
-        for i in range(len(coordinate1.attributes)):
-            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
-            attributes.append(att)
-	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
-	newVertex = Vertex(new_coordinate[0],new_coordinate[1],
-                           attributes = attributes)
-        mesh.maxVertexIndex+=1
-        newVertex.index = mesh.maxVertexIndex
-        mesh.meshVertices.append(newVertex)
-	return newVertex
-
-    def get_next_case(self, mesh,neighbour,triangle):
-	"""
-	Given the locations of two neighbouring triangles, 
-	examine the interior indices of their vertices (i.e. 
-	0,1 or 2) to determine what how the neighbour needs
-	to be refined.
-	"""
-	if (neighbour is None):
-		next_case = 'boundary'
-	else:
-		if triangle.vertices[1].x==neighbour.vertices[2].x:
-		    if triangle.vertices[1].y==neighbour.vertices[2].y:
-			next_case = 'vertex'
-                if triangle.vertices[1].x==neighbour.vertices[0].x:
-		    if triangle.vertices[1].y==neighbour.vertices[0].y:
-			next_case = 'two'
-		if triangle.vertices[1].x==neighbour.vertices[1].x:
-		    if triangle.vertices[1].y==neighbour.vertices[1].y:
-			next_case = 'one'
-	return next_case
-
-
-
-    def repairCaseVertex(self):
-
-        r = self.r
-
-
-        self.link(r[0],r[2])
-        self.repair(r[0])
-
-        self.link(r[1],r[3])
-        self.repair(r[1])
-
-        self.repair(r[2])
-
-        self.repair(r[3])
-
-
-    def repairCaseOne(self):
-        r = self.rkey
-
-
-        self.link(r[0],r[2])
-        self.repair(r[0])
-
-        self.link(r[1],r[4])
-        self.repair(r[1])
-
-        self.repair(r[4])
-
-    def repairCaseTwo(self):
-        r = self.r
-
-        self.link(r[0],r[4])
-        self.repair(r[0])
-
-        self.link(r[1],r[3])
-        self.repair(r[1])
-
-        self.repair(r[4])
-
-    def repairCaseBoundary(self):
-        r = self.r
-        self.repair(r[2])
-        self.repair(r[3])
-
-
-
-    def repair(self,triangle):
-        """
-        Given a triangle that knows its neighbours, this will
-        force the neighbours to comply.
-
-        However, it needs to compare the vertices of triangles
-        for this implementation 
-
-        But it doesn't work for invalid neighbour structures
-        """
-        n=triangle.neighbors
-        for i in (0,1,2):
-            if not n[i] is None:
-                for j in (0,1,2):#to find which side of the list is broken
-                    if not (n[i].vertices[j] in triangle.vertices):
-                    #ie if j is the side of n that needs fixing
-                        k = j
-                n[i].neighbors[k]=triangle
-
-
-
-    def link(self,triangle1,triangle2):
-        """
-        make triangle1 neighbors point to t
-                #count = 0riangle2
-        """
-        count = 0
-        for i in (0,1,2):#to find which side of the list is broken
-            if not (triangle1.vertices[i] in triangle2.vertices):
-                j = i
-                count+=1
-        assert count == 1
-        triangle1.neighbors[j]=triangle2
-
-class Discretised_Tuple_Set:
-    """
-    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
-    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
-    a[(10000)]=[] #NOT KEYERROR
-
-    a.append[(0.01)]
-    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
-
-    #NOT IMPLIMENTED
-    a.remove[(0.01)]
-    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
-
-    a.remove[(0.17)]
-    => {(0.0):[(0.02),(0.01)],0.2:[]}
-    #NOT IMPLIMENTED
-    at a.precision = 2:
-        a.round_up_rel[0.0]=
-        a.round_flat[0.0]=
-        a.round_down_rel[0.0]=
-
-        a.up((0.1,2.04))=
-
-    If t_rel = 0, nothing gets rounded into
-    two bins. If t_rel = 0.5, everything does.
-
-    Ideally, precision can be set high, so that multiple
-    entries are rarely in the same bin. And t_rel should
-    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
-    so that it is rare to put items in mutiple bins.
-
-    Ex bins per entry = product(a,b,c...,n)
-    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
-    b = 1 or 2 ... 
-
-    BUT!!! to avoid missing the right bin:
-    (-10)**(precision+1)*t_rel must be greater than the 
-    greatest possible variation that an identical element
-    can display.
-
-
-    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
-    but not .6 - this looks wrong, but note that *everything* will round,
-    so .6 wont be missed as everything close to it will check in .7 and .5.
-    """
-    def __init__(self,p_rel = 6,t_rel = 0.01):
-        self.__p_rel__ = p_rel
-        self.__t_rel__ = t_rel
-
-        self.__p_abs__ = p_rel+1
-        self.__t_abs__ = t_rel
-
-        assert t_rel <= 0.5
-        self.__items__ = {}
-        from math import frexp
-        self.frexp = frexp
-        roundings = [self.round_up_rel,\
-        self.round_down_rel,self.round_flat_rel,\
-        self.round_down_abs,self.round_up_abs,\
-        self.round_flat_abs]#
-
-        self.roundings = roundings
-
-    def __repr__(self):
-        return '%s'%self.__items__
-
-    def rounded_keys(self,key):
-        key = tuple(key)
-        keys = [key]
-        keys = self.__rounded_keys__(key)
-        return (keys)
-
-    def __rounded_keys__(self,key):
-        keys = []
-        rounded_key=list(key)
-        rounded_values=list(key)
-
-        roundings = list(self.roundings)
-
-        #initialise rounded_values
-        round = roundings.pop(0)
-        for i in range(len(rounded_values)):
-            rounded_key[i]=round(key[i])
-            rounded_values[i]={}
-            rounded_values[i][rounded_key[i]]=None
-        keys.append(tuple(rounded_key))
-        
-        for round in roundings:
-            for i in range(len(rounded_key)):
-                rounded_value=round(key[i])
-                if not rounded_values[i].has_key(rounded_value):
-                    #ie unless round_up_rel = round_down_rel
-                    #so the keys stay unique
-                    for j in range(len(keys)):
-                        rounded_key = list(keys[j])
-                        rounded_key[i]=rounded_value
-                        keys.append(tuple(rounded_key))
-        return keys
-
-    def append(self,item):
-        keys = self.rounded_keys(item)
-        for key in keys:
-            if self.__items__.has_key(key):
-                self.__items__[key].append(item)
-            else:
-                self.__items__[key]=[item]
-
-    def __getitem__(self,key):
-        answer = []
-        keys = self.rounded_keys(key)
-        for key in keys:
-            if self.__items__.has_key(key):
-                answer.extend(self.__items__[key])
-        #if len(answer)==0:
-        #    raise KeyError#FIXME or return KeyError
-        #                  #FIXME or just return []?
-        else:
-            return answer #FIXME or unique(answer)?
-
-    def __delete__(self,item):
-        keys = self.rounded_keys(item)
-        answer = False
-        #if any of the possible keys contains
-        #a list, return true
-        for key in keys:        
-            if self.__items__.has_key(key):
-                if item in self.__items__[key]:
-                    self.__items__[key].remove(item)
-
-    def remove(self,item):
-        self.__delete__(item)
-
-    def __contains__(self,item):
-
-        keys = self.rounded_keys(item)
-        answer = False
-        #if any of the possible keys contains
-        #a list, return true
-        for key in keys:        
-            if self.__items__.has_key(key):
-                if item in self.__items__[key]:
-                    answer = True
-        return answer
-
-
-    def has_item(self,item):
-        return self.__contains__(item)
-
-    def round_up_rel2(self,value):
-         t_rel=self.__t_rel__
-         #Rounding up the value
-         m,e = self.frexp(value)
-         m = m/2
-         e = e + 1
-         #m is the mantissa, e the exponent
-         # 0.5 < |m| < 1.0
-         m = m+t_rel*(10**-(self.__p_rel__))
-         #bump m up
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_down_rel2(self,value):
-         t_rel=self.__t_rel__
-         #Rounding down the value
-         m,e = self.frexp(value)
-         m = m/2
-         e = e + 1
-         #m is the mantissa, e the exponent
-         # 0.5 < m < 1.0
-         m = m-t_rel*(10**-(self.__p_rel__))
-         #bump the |m| down, by 5% or whatever
-         #self.p_rel dictates
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_flat_rel2(self,value):
-    #redundant
-         m,e = self.frexp(value)
-         m = m/2
-         e = e + 1
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_up_rel(self,value):
-         t_rel=self.__t_rel__
-         #Rounding up the value
-         m,e = self.frexp(value)
-         #m is the mantissa, e the exponent
-         # 0.5 < |m| < 1.0
-         m = m+t_rel*(10**-(self.__p_rel__))
-         #bump m up
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_down_rel(self,value):
-         t_rel=self.__t_rel__
-         #Rounding down the value
-         m,e = self.frexp(value)
-         #m is the mantissa, e the exponent
-         # 0.5 < m < 1.0
-         m = m-t_rel*(10**-(self.__p_rel__))
-         #bump the |m| down, by 5% or whatever
-         #self.p_rel dictates
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_flat_rel(self,value):
-    #redundant
-         m,e = self.frexp(value)
-         m = round(m,self.__p_rel__)
-         return m*(2.**e)
-
-    def round_up_abs(self,value):
-         t_abs=self.__t_abs__
-         #Rounding up the value
-         m = value+t_abs*(10**-(self.__p_abs__))
-         #bump m up
-         m = round(m,self.__p_abs__)
-         return m
-
-    def round_down_abs(self,value):
-         t_abs=self.__t_abs__
-         #Rounding down the value
-         m = value-t_abs*(10**-(self.__p_abs__))
-         #bump the |m| down, by 5% or whatever
-         #self.p_rel dictates
-         m = round(m,self.__p_abs__)
-         return m
-
-    def round_flat_abs(self,value):
-    #redundant?
-         m = round(value,self.__p_abs__)
-         return m
-
-    def keys(self):
-        return self.__items__.keys()
-
-
-class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
-    """
-    This is a discretised tuple set, but 
-    based on a mapping. The mapping MUST
-    return a sequence.
-
-    example: 
-    def weight(animal):
-        return [animal.weight]
-    
-    a = Mapped_Discretised_Tuple_Set(weight)
-    a.append[cow]
-    a.append[fox]
-    a.append[horse]
-
-    a[horse]    -> [cow,horse]
-    a[dog]      -> [fox]
-    a[elephant] -> []
-    """
-    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
-        Discretised_Tuple_Set.__init__\
-         (self,p_rel,t_rel = t_rel)
-        self.mapping = mapping
-
-    def rounded_keys(self,key):
-        mapped_key = tuple(self.mapping(key))
-        keys = self.__rounded_keys__(mapped_key)
-        return keys
-
-class Affine_Linespace(Mapped_Discretised_Tuple_Set):
-    """
-    The affine linespace creates a way to record and compare lines.
-    Precision is a bit of a hack, but it creates a way to avoid 
-    misses caused by round offs (between two lines of different
-    lenghts, the short one gets rounded off more). 
-    I am starting to think that a quadratic search would be faster.
-    Nearly.
-    """
-    def __init__(self,p_rel=4,t_rel=0.2):
-        Mapped_Discretised_Tuple_Set.__init__\
-            (self,self.affine_line,\
-            p_rel=p_rel,t_rel=t_rel)
-
-        roundings = \
-        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
-        self.roundings = roundings
-        #roundings = \
-        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
-        #self.roundings = roundings
-
-    def affine_line(self,line):
-        point_1 = line[0]
-        point_2 = line[1]
-        #returns the equation of a line
-        #between two points, in the from
-        #(a,b,-c), as in ax+by-c=0
-        #or line *dot* (x,y,1) = (0,0,0)
-
-        #Note that it normalises the line
-        #(a,b,-c) so Line*Line = 1.
-        #This does not change the mathematical
-        #properties, but it makes comparism
-        #easier.
-
-        #There are probably better algorithms.
-        x1 = point_1[0]
-        x2 = point_2[0]
-        y1 = point_1[1]
-        y2 = point_2[1]
-        dif_x = x1-x2
-        dif_y = y1-y2
-
-        if dif_x == dif_y == 0:
-            msg = 'points are the same'
-            raise msg
-        elif abs(dif_x)>=abs(dif_y):
-            alpha = (-dif_y)/dif_x
-            #a = alpha * b
-            b = -1.
-            c = (x1*alpha+x2*alpha+y1+y2)/2.
-            a = alpha*b
-        else:
-            beta = dif_x/(-dif_y)
-            #b = beta * a
-            a = 1.
-            c = (x1+x2+y1*beta+y2*beta)/2.
-            b = beta*a
-        mag = abs(a)+abs(b)
-        #This does not change the mathematical
-        #properties, but it makes comparism possible.
-
-        #note that the gradient is b/a, or (a/b)**-1.
-        #so 
-
-        #if a == 0:
-        #    sign_a = 1.
-        #else:
-        #    sign_a = a/((a**2)**0.5)
-        #if b == 0:
-        #    sign_b = 1.
-        #else:
-        #    sign_b = b/((b**2)**0.5)
-        #if c == 0:
-        #    sign_c = 1.
-        #else:
-        #    sign_c = c/((c**2)**0.5)
-        #a = a/mag*sign_a
-        #b = b/mag*sign_b
-        #c = c/mag*sign_c
-        a = a/mag
-        b = b/mag
-        c = c/mag
-        return a,b,c
-
 
 # FIXME (DSG-DSG)
Index: /anuga_core/source/anuga/pmesh/test_mesh.py
===================================================================
--- /anuga_core/source/anuga/pmesh/test_mesh.py	(revision 4905)
+++ /anuga_core/source/anuga/pmesh/test_mesh.py	(revision 4906)
@@ -36,20 +36,4 @@
                         'Point DistanceToPoint is wrong!')
         
-    def testTriangle(self):
-        a = Vertex (0.0, 0.0)
-        b = Vertex (0.0, 2.0)
-        c = Vertex (2.0,0.0)
-        d = Vertex (0.0, 4.0)
-        e = Vertex (2.0, 2.0)
-        f = Vertex (4.0,0.0)
-        
-        t1 = Triangle(b,a,c)       
-        t2 = Triangle(b,c,e)      
-        t3 = Triangle(e,c,f)      
-        t4 = Triangle(d,b,e)
-        t2.setNeighbors(t3,t4,t1)
-        
-        self.failUnless( t2.neighbors[2].vertices[0] == b, 'Triangle initialisation is wrong!')
-    
         
     def testSegment(self):
@@ -274,21 +258,4 @@
                         'Vertex deleted, instead of segment.')
 
-    def testTriangleArea(self):
-        a = Vertex (10.0, 10.0)
-        b = Vertex (10.0, 20.0)
-        c = Vertex (20.0,10.0)
-        
-        d = Vertex (-20.0, 0.0)
-        e = Vertex (-20.0, -20.0)
-        f = Vertex (20.0,-20.0)
-        
-        t1 = Triangle(b,a,c)      
-        t2 = Triangle(e,d,f)
-        
-#         print "t1", t1
-#         print "t1 area ", t1.calcArea()
-#         print "t2", t2
-#         print "t2 area ", t2.calcArea()
-        self.failUnless( t1.calcArea() == 50 and t2.calcArea() == 400, 'Triangle area is wrong!')
     def testisUserSegmentNew (self):
         mesh = Mesh()
@@ -1753,156 +1720,4 @@
 
 
-#___________beginning of Peters tests
-
-    def test_set_stuff(self):
-        """
-        Documentation
-        """
-        #making a test mesh
-        p0=[0.,2.]
-        p1=[1.,2.]
-        p2=[0.,1.]
-        p3=[1.,1.]
-        p4=[0.,0.]
-        p5=[2.,0.]
-        p6=[-1.,1.]
-        point_list = [p0,p1,p2,p3,p4,p5,p6]
-
-        a0=[0]
-        a1=[0]
-        a2=[100]
-        a3=[0]
-        a4=[0]
-        a5=[0]
-        a6=[0]
-        attribute=[a0,a1,a2,a3,a4,a5,a6]
-        
-        t0=[0,3,1]
-        t1=[0,2,3]
-        t2=[2,4,3]
-        t3=[4,5,3]
-        t4=[1,3,5]
-        t5=[2,6,4]
-
-        n0=[4,-1,2]
-        n1=[2,0,-1]
-        n2=[3,1,5]
-        n3=[4,2,-1]
-        n4=[3,-1,0]
-        n5=[-1,2,-1]
-
-        tri_list = [t0,t1,t2,t3,t4,t5]
-        n_list = [n0,n1,n2,n3,n4,n5]
-        for i in range(6):
-            for j in (0,1,2):
-                a=attribute[tri_list[i][j]]
-                tri_list[i][j]=point_list[tri_list[i][j]]
-                tri_list[i][j]=Vertex(tri_list[i][j][0]\
-                                      ,tri_list[i][j][1],a)
-            neighbours=n_list[i]
-            tri_list[i]=Triangle(tri_list[i][0],\
-                                 tri_list[i][1],tri_list[i][2]\
-                                ,neighbors=neighbours)
-
-        #testing selectAll
-        mesh = Mesh()
-        mesh.attributeTitles=['attrib']
-
-        mesh.meshTriangles=tri_list
-
-        mesh.selectAllTriangles()
-        A=mesh.sets[mesh.setID['All']]
-        assert list_comp(tri_list,A)
-
-       #testing threshold
-        mesh = Mesh()
-        mesh.attributeTitles=['attrib']
-
-        mesh.meshTriangles=tri_list
-        mesh.selectAllTriangles()
-        mesh.threshold('All',min=30,max=35,attribute_name = 'attrib')
-        A = [tri_list[1],tri_list[2],tri_list[5]]
-        B = mesh.sets[mesh.setID['All']]
-        assert list_comp(A,B)
-       
-
-        A = [tri_list[3],tri_list[2],tri_list[5]]
-        assert not list_comp(A,B)
-
-        #testing 
-
-    def test_Discretised_Tuple_Set_rounding(self):
-        #This is the hardest bit of DST
-
-        tol = 0.1
-        a=Discretised_Tuple_Set(p_rel=1,t_rel= tol)
-        m = 0.541
-        m_up = 0.6
-        m_down = 0.5
-        assert m_up == a.round_up_rel(m)
-        assert m_down == a.round_down_rel(m)
-
-        tol = 0.1
-        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
-        m = 0.539
-        m_up = 0.5
-        m_down = 0.5
-        assert m_up == a.round_up_rel(m)
-        assert m_down == a.round_down_rel(m)
-
-        tol = 0.5
-        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
-
-
-        m = 0.6
-        m_up = 0.7
-        m_down = 0.5
-        assert m_up == a.round_up_rel(m)
-        assert m_down == a.round_down_rel(m)
-
-        m = 0.599
-        m_up = 0.6
-        m_down = 0.5
-        assert m_up == a.round_up_rel(m)
-        assert m_down == a.round_down_rel(m)
-
-    def test_Discretised_Tuple_Set_get(self):
-        
-        tol = 0.25
-        a=Discretised_Tuple_Set(p_rel=1,t_rel = tol)
-        b = (1.1,1.1)
-        a.append(b)
-        list = [(1.2,1.),(1.,1.),(1.,1.2),(1.2,1.2)]
-        for key in list:
-            assert a[key][0]==b
-            assert len(a[key])==1
-        
-        c = (2.1,1.)
-        a.append(c)
-        assert a[(2.,1.)][0]==c
-        assert a[(2.2,1.)][0]==c
-
-    def test_mapped_Discretised_Tuple_Set(self):
-
-        def map(sequence):
-            return [len(sequence)]
-
-        tol = 0.5
-        a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
-        b = range(20)
-        a.append(b)
-        assert b in a[range(17)] 
-        assert b in a[range(22)]
-
-        tol = 0.01
-        a=Mapped_Discretised_Tuple_Set(map,p_rel=1,t_rel = tol)
-        b = range(20)
-        a.append(b)
-        assert b in a[range(20)] 
-        assert b in a[range(19)] 
-        assert not range(17) in a
-
-#___________end of Peters tests
-
     def test_add_region_from_polygon(self):
         m=Mesh()
Index: /anuga_core/source/obsolete_code/obs_mesh.py
===================================================================
--- /anuga_core/source/obsolete_code/obs_mesh.py	(revision 4906)
+++ /anuga_core/source/obsolete_code/obs_mesh.py	(revision 4906)
@@ -0,0 +1,1321 @@
+
+SET_COLOUR='red'
+
+class Triangle(MeshObject):
+    """
+    A triangle element, defined by 3 vertices.
+    Attributes based on the Triangle program.
+    """
+
+    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
+        """
+        Vertices, the initial arguments, are listed in counterclockwise order.
+        """
+        self.vertices= [vertex1,vertex2, vertex3 ]
+        
+        if attribute is None:
+            self.attribute =""
+        else:
+            self.attribute = attribute #this is a string
+            
+        if neighbors is None:
+            self.neighbors=[]
+        else:
+            self.neighbors=neighbors
+
+    def replace(self,new_triangle):
+        self = new_triangle
+
+    def longestSideID(self):
+        ax = self.vertices[0].x
+        ay = self.vertices[0].y
+        
+        bx = self.vertices[1].x
+        by = self.vertices[1].y
+        
+        cx = self.vertices[2].x
+        cy = self.vertices[2].y
+
+        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
+        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
+        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
+ 
+        len = [lenA,lenB,lenC]
+        return len.index(max(len))
+
+    def rotate(self,offset):
+        """
+        permute the order of the sides of the triangle
+        offset must be 0,1 or 2
+        """
+
+        if offset == 0:
+            pass
+        else:
+            if offset == 1:
+                self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
+                self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
+            if offset == 2:
+                self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
+                self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
+
+    def rotate_longest_side(self):
+        self.rotate(self.longestSideID())
+
+    def getVertices(self):
+        return self.vertices
+    
+    def get_vertices(self):
+        """
+        Return a list of the vertices.  The x and y values will be relative
+        Easting and Northings for the zone of the current geo_ref.
+        """
+        return self.vertices
+    
+    def calcArea(self):
+        ax = self.vertices[0].x
+        ay = self.vertices[0].y
+        
+        bx = self.vertices[1].x
+        by = self.vertices[1].y
+        
+        cx = self.vertices[2].x
+        cy = self.vertices[2].y
+        
+        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
+    
+    def calcP(self):
+        #calculate the perimeter
+        ax = self.vertices[0].x
+        ay = self.vertices[0].y
+        
+        bx = self.vertices[1].x
+        by = self.vertices[1].y
+        
+        cx = self.vertices[2].x
+        cy = self.vertices[2].y
+
+        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
+        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
+        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
+
+        return a+b+c
+            
+    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
+        """
+        neighbor1 is the triangle opposite vertex1 and so on.
+        Null represents no neighbor
+        """
+        self.neighbors = [neighbor1, neighbor2, neighbor3]
+
+    def setAttribute(self,attribute):
+        """
+        neighbor1 is the triangle opposite vertex1 and so on.
+        Null represents no neighbor
+        """
+        self.attribute = attribute #this is a string
+        
+    def __repr__(self):
+        return "[%s,%s]" % (self.vertices,self.attribute)
+        
+
+    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
+        """
+        Draw a triangle, returning the objectID
+        """
+        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
+                                     scale*-1*(self.vertices[1].y + yoffset),
+                                     scale*(self.vertices[0].x + xoffset),
+                                     scale*-1*(self.vertices[0].y + yoffset),
+                                     scale*(self.vertices[2].x + xoffset),
+                                     scale*-1*(self.vertices[2].y + yoffset),
+                                     tags = tags,
+                                     outline = colour,fill = '')
+
+
+        self.setID={}
+        #a dictionary of names.
+        #multiple sets are allowed, but the gui does not yet
+        #support this
+        
+        self.setID['None']=0
+        #contains the names of the sets pointing to the indexes
+        #in the list. 
+        
+	self.sets=[[]]
+        #Contains the lists of triangles (triangle sets)
+
+      ##################
+
+        
+    def refineSet(self,setName):
+        Triangles = self.sets[self.setID[setName]]
+        Refine(self,Triangles)
+
+    def selectAllTriangles(self):
+        A=[]
+        A.extend(self.meshTriangles)
+        if not('All' in self.setID.keys()):
+            self.setID['All']=len(self.sets)
+            self.sets.append(A)
+        else:
+            self.sets[self.setID['All']]=A
+        return 'All'
+        # and objectIDs
+
+
+    def clearSelection(self):
+        A = []
+        if not('None' in self.setID.keys()):
+            self.setID['None']=len(self.sets)
+            self.sets.append(A)
+        return 'None'
+
+    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
+    #FIXME Draws over previous triangles - may bloat canvas
+        Triangles = self.sets[self.setID[setName]]
+        for triangle in Triangles:
+            triangle.draw(canvas,1,
+                          scale = SCALE,
+                          colour = colour)
+	    
+    def undrawSet(self,canvas,setName,SCALE,colour='green'):
+    #FIXME Draws over previous lines - may bloat canvas
+        Triangles = self.sets[self.setID[setName]]
+        for triangle in Triangles:
+            triangle.draw(canvas,1,
+                          scale = SCALE,
+                          colour = colour)
+
+    def weed(self,Vertices,Segments):
+        #Depreciated
+        #weed out existing duplicates
+        print 'len(self.getUserSegments())'
+        print len(self.getUserSegments())
+        print 'len(self.getUserVertices())'
+        print len(self.getUserVertices())
+
+        point_keys = {}
+        for vertex in Vertices:
+            point = (vertex.x,vertex.y)
+            point_keys[point]=vertex
+        #inlined would looks very ugly
+
+        line_keys = {}
+        for segment in Segments:
+            vertex1 = segment.vertices[0]
+            vertex2 = segment.vertices[1]
+            point1 = (vertex1.x,vertex1.y)
+            point2 = (vertex2.x,vertex2.y)
+            segment.vertices[0]=point_keys[point1]
+            segment.vertices[1]=point_keys[point2]
+            vertex1 = segment.vertices[0]
+            vertex2 = segment.vertices[1]
+            point1 = (vertex1.x,vertex1.y)
+            point2 = (vertex2.x,vertex2.y)
+            line1 = (point1,point2)
+            line2 = (point2,point1)
+            if not (line_keys.has_key(line1) \
+                 or line_keys.has_key(line2)):
+                 line_keys[line1]=segment
+        Vertices=point_keys.values()
+        Segments=line_keys.values()
+        return Vertices,Segments
+
+    def segs_to_dict(self,segments):
+        dict={}
+        for segment in segments:
+            vertex1 = segment.vertices[0]
+            vertex2 = segment.vertices[1]
+            point1 = (vertex1.x,vertex1.y)
+            point2 = (vertex2.x,vertex2.y)
+            line = (point1,point2)
+            dict[line]=segment
+        return dict
+
+    def seg2line(self,s):
+        return ((s.vertices[0].x,s.vertices[0].y,)\
+                (s.vertices[1].x,s.vertices[1].y))
+
+    def line2seg(self,line,tag=None):
+        point0 = self.point2ver(line[0])
+        point1 = self.point2ver(line[1])
+        return Segment(point0,point1,tag=tag)
+
+    def ver2point(self,vertex):
+        return (vertex.x,vertex.y)
+
+    def point2ver(self,point):
+        return Vertex(point[0],point[1])
+
+    def smooth_polySet(self,min_radius=0.05):
+        #for all pairs of connecting segments:
+        #    propose a new segment that replaces the 2
+
+        #    If the difference between the new segment
+        #    and the old lines is small: replace the
+        #    old lines.
+
+        seg2line = self.seg2line
+        ver2point= self.ver2point
+        line2seg = self.line2seg
+        point2ver= self.point2ver
+
+        #create dictionaries of lines -> segments
+        userSegments = self.segs_to_dict(self.userSegments)
+        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
+
+        #lump user and alpha segments
+        for key in alphaSegments.keys():
+            userSegments[key]=alphaSegments[key]
+
+        #point_keys = tuple -> vertex
+        #userVertices = vertex -> [line,line] - lines from that node
+        point_keys = {}
+        userVertices={}
+        for vertex in self.getUserVertices():
+            point = ver2point(vertex)
+            if not point_keys.has_key(point):
+                point_keys[point]=vertex
+                userVertices[vertex]=[]
+        for key in userSegments.keys():
+            line = key
+            point_0 = key[0]
+            point_1 = key[1]
+            userVertices[point_keys[point_0]].append(line)
+            userVertices[point_keys[point_1]].append(line)
+
+        for point in point_keys.keys():
+            try:
+            #removed keys can cause keyerrors
+                vertex = point_keys[point]
+		lines = userVertices[vertex]
+    
+                #if there are 2 lines on the node
+		if len(lines)==2:
+		    line_0 = lines[0]
+		    line_1 = lines[1]
+    
+                    #if the tags are the the same on the 2 lines
+		    if userSegments[line_0].tag == userSegments[line_1].tag:
+			tag = userSegments[line_0].tag 
+    
+                        #point_a is one of the next nodes, point_b is the other
+			if point==line_0[0]:
+			    point_a = line_0[1]
+			if point==line_0[1]:
+			    point_a = line_0[0]
+			if point==line_1[0]:
+			    point_b = line_1[1]
+			if point==line_1[1]:
+			    point_b = line_1[0]
+    
+    
+                        #line_2 is proposed
+			line_2 = (point_a,point_b)
+
+			#calculate the area of the triangle between
+			#the two existing segments and the proposed
+			#new segment
+			ax = point_a[0]
+			ay = point_a[1]
+			bx = point_b[0]
+			by = point_b[1]
+			cx = point[0]
+			cy = point[1]
+			area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
+
+			#calculate the perimeter
+			len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
+			len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
+			len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
+			perimeter = len_a+len_b+len_c
+
+			#calculate the radius
+			r = area/(2*perimeter)
+
+			#if the radius is small: then replace the existing
+			#segments with the new one
+			if r < min_radius:
+			    if len_c < min_radius: append = False
+			    else: append = True
+			    #if the new seg is also time, don't add it
+			    if append:
+				segment = self.line2seg(line_2,tag=tag)
+
+			    list_a=userVertices[point_keys[point_a]]
+			    list_b=userVertices[point_keys[point_b]]
+
+			    if line_0 in list_a:
+				list_a.remove(line_0)
+			    else:
+				list_a.remove(line_1)
+
+			    if line_0 in list_b:
+				list_b.remove(line_0)
+			    else:
+				list_b.remove(line_1)
+
+			    if append:
+				list_a.append(line_2)
+				list_b.append(line_2)
+			    else:
+				if len(list_a)==0:
+				    userVertices.pop(point_keys[point_a])
+				    point_keys.pop(point_a)
+				if len(list_b)==0:
+				    userVertices.pop(point_keys[point_b])
+				    point_keys.pop(point_b)
+
+			    userVertices.pop(point_keys[point])
+			    point_keys.pop(point)
+			    userSegments.pop(line_0)
+			    userSegments.pop(line_1)
+
+			    if append:
+				userSegments[line_2]=segment
+            except:
+                pass
+
+        #self.userVerticies = userVertices.keys()
+        #self.userSegments = []
+        #for key in userSegments.keys():
+        #    self.userSegments.append(userSegments[key])
+        #self.alphaUserSegments = []
+
+        self.userVerticies = []
+        self.userSegments = []
+        self.alphaUserSegments = []
+
+        return userVertices,userSegments,alphaSegments
+
+    def triangles_to_polySet(self,setName):
+        #self.smooth_polySet()
+
+        seg2line = self.seg2line
+        ver2point= self.ver2point
+        line2seg = self.line2seg
+        point2ver= self.point2ver
+
+        from Numeric import array,allclose
+        #turn the triangles into a set
+        Triangles = self.sets[self.setID[setName]]
+        Triangles_dict = {}
+        for triangle in Triangles:
+            Triangles_dict[triangle]=None
+ 
+
+        #create a dict of points to vertexes (tuple -> object)
+        #also create a set of vertexes (object -> True)
+        point_keys = {}
+        userVertices={}
+        for vertex in self.getUserVertices():
+            point = ver2point(vertex)
+            if not point_keys.has_key(point):
+                point_keys[point]=vertex
+                userVertices[vertex]=True
+
+        #create a dict of lines to segments (tuple -> object)
+        userSegments = self.segs_to_dict(self.userSegments)
+        #append the userlines in an affine linespace
+        affine_lines = Affine_Linespace()
+        for line in userSegments.keys():
+            affine_lines.append(line)
+        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
+        for line in alphaSegments.keys():
+            affine_lines.append(line)
+        
+        for triangle in Triangles:
+            for i in (0,1,2):
+                #for every triangles neighbour:
+                if not Triangles_dict.has_key(triangle.neighbors[i]):
+                #if the neighbour is not in the set:
+                    a = triangle.vertices[i-1]
+                    b = triangle.vertices[i-2]
+                    #Get possible matches:
+                    point_a = ver2point(a)
+                    point_b = ver2point(b)
+                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
+                    line = (point_a,point_b)
+                    tag = None
+
+
+                    #this bit checks for matching lines
+                    possible_lines = affine_lines[line] 
+                    possible_lines = unique(possible_lines)
+                    found = 0                            
+                    for user_line in possible_lines:
+                        if self.point_on_line(midpoint,user_line):
+                            found+=1
+                            assert found<2
+                            if userSegments.has_key(user_line):
+                                parent_segment = userSegments.pop(user_line)
+                            if alphaSegments.has_key(user_line):
+                                parent_segment = alphaSegments.pop(user_line)
+                            tag = parent_segment.tag
+                            offspring = [line]
+                            offspring.extend(self.subtract_line(user_line,
+                                                                line))
+                            affine_lines.remove(user_line)
+                            for newline in offspring:
+                                line_vertices = []
+                                for point in newline:
+                                    if point_keys.has_key(point):
+                                        vert = point_keys[point]
+                                    else:
+                                        vert = Vertex(point[0],point[1])
+                                        userVertices[vert]=True
+                                        point_keys[point]=vert
+                                    line_vertices.append(vert)
+                                segment = Segment(line_vertices[0],
+                                                  line_vertices[1],tag)
+                                userSegments[newline]=segment
+                                affine_lines.append(newline)
+                            #break
+                    assert found<2
+
+
+
+                    #if no matching lines
+                    if not found:
+                        line_vertices = []
+		        for point in line:
+			    if point_keys.has_key(point):
+			        vert = point_keys[point]
+			    else:
+				vert = Vertex(point[0],point[1])
+				userVertices[vert]=True
+				point_keys[point]=vert
+			    line_vertices.append(vert)
+			segment = Segment(line_vertices[0],
+                                          line_vertices[1],tag)
+                        userSegments[line]=segment
+                        affine_lines.append(line)
+        
+        self.userVerticies = []
+        self.userSegments = []
+        self.alphaUserSegments = []
+
+        return userVertices,userSegments,alphaSegments
+
+    def subtract_line(self,parent,child):
+        #Subtracts child from parent
+        #Requires that the child is a 
+        #subline of parent to work.
+
+        from Numeric import allclose,dot,array
+        A= parent[0]
+        B= parent[1]
+        a = child[0]
+        b = child[1]
+
+        A_array = array(parent[0])
+        B_array = array(parent[1])
+        a_array   = array(child[0])
+        b_array   = array(child[1])
+
+        assert not A == B
+        assert not a == b
+
+        answer = []
+
+        #if the new line does not share a 
+        #vertex with the old one
+        if not (allclose(A_array,a_array)\
+             or allclose(B_array,b_array)\
+             or allclose(A_array,b_array)\
+             or allclose(a_array,B_array)):
+            if dot(A_array-a_array,A_array-a_array) \
+            < dot(A_array-b_array,A_array-b_array):
+                sibling1 = (A,a)
+                sibling2 = (B,b)
+                return [sibling1,sibling2]
+            else:
+                sibling1 = (A,b)
+                sibling2 = (B,a)
+                return [sibling1,sibling2]
+
+        elif allclose(A_array,a_array):
+            if allclose(B_array,b_array):
+                return []
+	    else:
+                sibling = (b,B)
+                return [sibling]
+    	elif allclose(B_array,b_array):
+            sibling = (a,A)
+            return [sibling]
+
+        elif allclose(A_array,b_array):
+            if allclose(B,a):
+                return []
+	    else:
+                sibling = (a,B)
+                return [sibling]
+        elif allclose(a_array,B_array):
+            sibling = (b,A)
+            return [sibling]
+
+    def point_on_line(self,point,line):
+        #returns true within a tolerance of 3 degrees
+        x=point[0]
+        y=point[1]
+        x0=line[0][0]
+        x1=line[1][0]
+        y0=line[0][1]
+        y1=line[1][1]
+        from Numeric import array, dot, allclose
+        from math import sqrt
+        tol = 3. #DEGREES
+        tol = tol*3.1415/180
+
+        a = array([x - x0, y - y0]) 
+        a_normal = array([a[1], -a[0]])      
+        len_a_normal = sqrt(sum(a_normal**2)) 
+
+        b = array([x1 - x0, y1 - y0])  	       
+        len_b = sqrt(sum(b**2)) 
+    
+        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
+            #Point is somewhere on the infinite extension of the line
+
+            len_a = sqrt(sum(a**2))     	       	                
+            if dot(a, b) >= 0 and len_a <= len_b:
+               return True
+            else:   
+               return False
+        else:
+          return False
+
+    def line_length(self,line):
+        x0=line[0][0]
+        x1=line[1][0]
+        y0=line[0][1]
+        y1=line[1][1]
+        return ((x1-x0)**2-(y1-y0)**2)**0.5     
+
+    def threshold(self,setName,min=None,max=None,attribute_name='elevation'):
+        """
+        threshold using  d
+        """
+        triangles = self.sets[self.setID[setName]]
+        A = []
+
+        if attribute_name in self.attributeTitles:
+            i = self.attributeTitles.index(attribute_name)
+        else: i = -1#no attribute
+        if not max == None:
+            for t in triangles:
+                if (min<self.av_att(t,i)<max):
+                    A.append(t)
+        else:
+            for t in triangles:
+                if (min<self.av_att(t,i)):
+                    A.append(t)
+        self.sets[self.setID[setName]] = A
+
+    def general_threshold(self,setName,min=None,max=None\
+              ,attribute_name = 'elevation',function=None):
+        """
+        Thresholds the triangles 
+        """
+        from visual.graph import arange,ghistogram,color as colour
+        triangles = self.sets[self.setID[setName]]
+        A = []
+        data=[]
+        #data is for the graph
+
+        if attribute_name in self.attributeTitles:
+            i = self.attributeTitles.index(attribute_name)
+        else: i = -1
+        if not max == None:
+            for t in triangles:
+                value=function(t,i)
+                if (min<value<max):
+                    A.append(t)
+                data.append(value)
+        else:
+            for t in triangles:
+                value=function(t,i)
+                if (min<value):
+                    A.append(t)
+                data.append(value)
+        self.sets[self.setID[setName]] = A
+
+        if self.visualise_graph:
+            if len(data)>0:
+                max=data[0]
+                min=data[0]
+                for value in data:
+                    if value > max:
+                        max = value
+                    if value < min:
+                        min = value
+
+                inc = (max-min)/100
+
+                histogram = ghistogram(bins=arange(min,max,inc),\
+                             color = colour.red)
+                histogram.plot(data=data)
+        
+    def av_att(self,triangle,i):
+        if i==-1: return 1
+        else:
+            #evaluates the average attribute of the vertices of a triangle.
+            V = triangle.getVertices()
+            a0 = (V[0].attributes[i])
+            a1 = (V[1].attributes[i])
+            a2 = (V[2].attributes[i])
+            return (a0+a1+a2)/3
+
+    def Courant_ratio(self,triangle,index):
+        """
+        Uses the courant threshold
+        """
+        e = self.av_att(triangle,index)
+        A = triangle.calcArea()
+        P = triangle.calcP()
+        r = A/(2*P)
+        e = max(0.1,abs(e))
+        return r/e**0.5
+
+    def Gradient(self,triangle,index):
+        V = triangle.vertices
+        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
+        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
+        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
+            print ((grad_x**2)+(grad_y**2))**(0.5)
+        return ((grad_x**2)+(grad_y**2))**(0.5)
+    
+
+    def append_triangle(self,triangle):
+        self.meshTriangles.append(triangle)
+
+    def replace_triangle(self,triangle,replacement):
+        i = self.meshTriangles.index(triangle)
+        self.meshTriangles[i]=replacement
+        assert replacement in self.meshTriangles
+
+ 
+"""Refines triangles
+
+   Implements the #triangular bisection?# algorithm.
+ 
+   
+"""
+
+def Refine(mesh, triangles):
+    """
+    Given a general_mesh, and a triangle number, split
+    that triangle in the mesh in half. Then to prevent
+    vertices and edges from meeting, keep refining 
+    neighbouring triangles until the mesh is clean.
+    """
+    state = BisectionState(mesh)
+    for triangle in triangles:
+        if not state.refined_triangles.has_key(triangle):
+            triangle.rotate_longest_side()
+            state.start(triangle)
+            Refine_mesh(mesh, state)
+
+def Refine_mesh(mesh, state):
+    """
+    """
+    state.getState(mesh)
+    refine_triangle(mesh,state)
+    state.evolve()
+    if not state.end:
+        Refine_mesh(mesh,state)
+
+def refine_triangle(mesh,state):
+    split(mesh,state.current_triangle,state.new_point)
+    if state.case == 'one':
+        state.r[3]=state.current_triangle#triangle 2
+
+        new_triangle_id = len(mesh.meshTriangles)-1
+        new_triangle = mesh.meshTriangles[new_triangle_id]
+
+        split(mesh,new_triangle,state.old_point)
+        state.r[2]=new_triangle#triangle 1.2
+        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
+        r = state.r
+        state.repairCaseOne()
+
+    if state.case == 'two':
+        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
+
+        new_triangle = state.current_triangle
+
+        split(mesh,new_triangle,state.old_point)
+
+        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
+        state.r[4]=new_triangle#triangle 2.2
+        r = state.r
+
+        state.repairCaseTwo()
+
+    if state.case == 'vertex':
+        state.r[2]=state.current_triangle#triangle 2
+        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
+        r = state.r
+        state.repairCaseVertex()
+        
+    if state.case == 'start':
+        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
+        state.r[3]=state.current_triangle#triangle 2
+
+    if state.next_case == 'boundary':
+        state.repairCaseBoundary()
+
+
+def split(mesh, triangle, new_point):
+	"""
+	Given a mesh, triangle_id and a new point,
+	split the corrosponding triangle into two
+	new triangles and update the mesh.
+	"""
+
+	new_triangle1 = Triangle(new_point,triangle.vertices[0],
+                                 triangle.vertices[1],
+                                 attribute = triangle.attribute,
+                                 neighbors = None)
+	new_triangle2 = Triangle(new_point,triangle.vertices[2],
+                                 triangle.vertices[0],
+                                 attribute = triangle.attribute,
+                                 neighbors = None)
+
+        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
+        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
+
+        mesh.meshTriangles.append(new_triangle1)
+
+        triangle.vertices = new_triangle2.vertices
+        triangle.neighbors = new_triangle2.neighbors
+
+
+class State:
+
+    def __init__(self):
+        pass
+
+class BisectionState(State):
+
+
+    def __init__(self,mesh):
+        self.len = len(mesh.meshTriangles)
+        self.refined_triangles = {}
+        self.mesh = mesh
+        self.current_triangle = None
+        self.case = 'start'
+        self.end = False
+        self.r = [None,None,None,None,None]
+
+    def start(self, triangle):
+        self.current_triangle = triangle
+        self.case = 'start'
+        self.end = False
+        self.r = [None,None,None,None,None]
+
+    def getState(self,mesh):
+        if not self.case == 'vertex':
+            self.new_point=self.getNewVertex(mesh, self.current_triangle)
+            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
+            self.neighbour = self.current_triangle.neighbors[0]
+            if not self.neighbour is None:
+                self.neighbour.rotate_longest_side()
+            self.next_case = self.get_next_case(mesh,self.neighbour,
+                                                self.current_triangle)
+        if self.case == 'vertex':
+            self.new_point=self.old_point
+
+
+    def evolve(self):
+        if self.case == 'vertex':
+            self.end = True
+
+        self.last_case = self.case
+        self.case = self.next_case
+
+        self.old_point = self.new_point
+        self.current_triangle = self.neighbour
+
+        if self.case == 'boundary':
+            self.end = True
+        self.refined_triangles[self.r[2]]=1
+        self.refined_triangles[self.r[3]]=1
+        if not self.r[4] is None:
+            self.refined_triangles[self.r[4]]=1
+        self.r[0]=self.r[2]
+        self.r[1]=self.r[3]
+
+
+    def getNewVertex(self,mesh,triangle):
+	coordinate1 = triangle.vertices[1]
+	coordinate2 = triangle.vertices[2]
+	a = ([coordinate1.x*1.,coordinate1.y*1.])
+	b = ([coordinate2.x*1.,coordinate2.y*1.])
+        attributes = []
+        for i in range(len(coordinate1.attributes)):
+            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
+            attributes.append(att)
+	new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
+	newVertex = Vertex(new_coordinate[0],new_coordinate[1],
+                           attributes = attributes)
+        mesh.maxVertexIndex+=1
+        newVertex.index = mesh.maxVertexIndex
+        mesh.meshVertices.append(newVertex)
+	return newVertex
+
+    def get_next_case(self, mesh,neighbour,triangle):
+	"""
+	Given the locations of two neighbouring triangles, 
+	examine the interior indices of their vertices (i.e. 
+	0,1 or 2) to determine what how the neighbour needs
+	to be refined.
+	"""
+	if (neighbour is None):
+		next_case = 'boundary'
+	else:
+		if triangle.vertices[1].x==neighbour.vertices[2].x:
+		    if triangle.vertices[1].y==neighbour.vertices[2].y:
+			next_case = 'vertex'
+                if triangle.vertices[1].x==neighbour.vertices[0].x:
+		    if triangle.vertices[1].y==neighbour.vertices[0].y:
+			next_case = 'two'
+		if triangle.vertices[1].x==neighbour.vertices[1].x:
+		    if triangle.vertices[1].y==neighbour.vertices[1].y:
+			next_case = 'one'
+	return next_case
+
+
+
+    def repairCaseVertex(self):
+
+        r = self.r
+
+
+        self.link(r[0],r[2])
+        self.repair(r[0])
+
+        self.link(r[1],r[3])
+        self.repair(r[1])
+
+        self.repair(r[2])
+
+        self.repair(r[3])
+
+
+    def repairCaseOne(self):
+        r = self.rkey
+
+
+        self.link(r[0],r[2])
+        self.repair(r[0])
+
+        self.link(r[1],r[4])
+        self.repair(r[1])
+
+        self.repair(r[4])
+
+    def repairCaseTwo(self):
+        r = self.r
+
+        self.link(r[0],r[4])
+        self.repair(r[0])
+
+        self.link(r[1],r[3])
+        self.repair(r[1])
+
+        self.repair(r[4])
+
+    def repairCaseBoundary(self):
+        r = self.r
+        self.repair(r[2])
+        self.repair(r[3])
+
+
+
+    def repair(self,triangle):
+        """
+        Given a triangle that knows its neighbours, this will
+        force the neighbours to comply.
+
+        However, it needs to compare the vertices of triangles
+        for this implementation 
+
+        But it doesn't work for invalid neighbour structures
+        """
+        n=triangle.neighbors
+        for i in (0,1,2):
+            if not n[i] is None:
+                for j in (0,1,2):#to find which side of the list is broken
+                    if not (n[i].vertices[j] in triangle.vertices):
+                    #ie if j is the side of n that needs fixing
+                        k = j
+                n[i].neighbors[k]=triangle
+
+
+
+    def link(self,triangle1,triangle2):
+        """
+        make triangle1 neighbors point to t
+                #count = 0riangle2
+        """
+        count = 0
+        for i in (0,1,2):#to find which side of the list is broken
+            if not (triangle1.vertices[i] in triangle2.vertices):
+                j = i
+                count+=1
+        assert count == 1
+        triangle1.neighbors[j]=triangle2
+
+class Discretised_Tuple_Set:
+    """
+    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
+    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
+    a[(10000)]=[] #NOT KEYERROR
+
+    a.append[(0.01)]
+    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
+
+    #NOT IMPLIMENTED
+    a.remove[(0.01)]
+    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
+
+    a.remove[(0.17)]
+    => {(0.0):[(0.02),(0.01)],0.2:[]}
+    #NOT IMPLIMENTED
+    at a.precision = 2:
+        a.round_up_rel[0.0]=
+        a.round_flat[0.0]=
+        a.round_down_rel[0.0]=
+
+        a.up((0.1,2.04))=
+
+    If t_rel = 0, nothing gets rounded into
+    two bins. If t_rel = 0.5, everything does.
+
+    Ideally, precision can be set high, so that multiple
+    entries are rarely in the same bin. And t_rel should
+    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
+    so that it is rare to put items in mutiple bins.
+
+    Ex bins per entry = product(a,b,c...,n)
+    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
+    b = 1 or 2 ... 
+
+    BUT!!! to avoid missing the right bin:
+    (-10)**(precision+1)*t_rel must be greater than the 
+    greatest possible variation that an identical element
+    can display.
+
+
+    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
+    but not .6 - this looks wrong, but note that *everything* will round,
+    so .6 wont be missed as everything close to it will check in .7 and .5.
+    """
+    def __init__(self,p_rel = 6,t_rel = 0.01):
+        self.__p_rel__ = p_rel
+        self.__t_rel__ = t_rel
+
+        self.__p_abs__ = p_rel+1
+        self.__t_abs__ = t_rel
+
+        assert t_rel <= 0.5
+        self.__items__ = {}
+        from math import frexp
+        self.frexp = frexp
+        roundings = [self.round_up_rel,\
+        self.round_down_rel,self.round_flat_rel,\
+        self.round_down_abs,self.round_up_abs,\
+        self.round_flat_abs]#
+
+        self.roundings = roundings
+
+    def __repr__(self):
+        return '%s'%self.__items__
+
+    def rounded_keys(self,key):
+        key = tuple(key)
+        keys = [key]
+        keys = self.__rounded_keys__(key)
+        return (keys)
+
+    def __rounded_keys__(self,key):
+        keys = []
+        rounded_key=list(key)
+        rounded_values=list(key)
+
+        roundings = list(self.roundings)
+
+        #initialise rounded_values
+        round = roundings.pop(0)
+        for i in range(len(rounded_values)):
+            rounded_key[i]=round(key[i])
+            rounded_values[i]={}
+            rounded_values[i][rounded_key[i]]=None
+        keys.append(tuple(rounded_key))
+        
+        for round in roundings:
+            for i in range(len(rounded_key)):
+                rounded_value=round(key[i])
+                if not rounded_values[i].has_key(rounded_value):
+                    #ie unless round_up_rel = round_down_rel
+                    #so the keys stay unique
+                    for j in range(len(keys)):
+                        rounded_key = list(keys[j])
+                        rounded_key[i]=rounded_value
+                        keys.append(tuple(rounded_key))
+        return keys
+
+    def append(self,item):
+        keys = self.rounded_keys(item)
+        for key in keys:
+            if self.__items__.has_key(key):
+                self.__items__[key].append(item)
+            else:
+                self.__items__[key]=[item]
+
+    def __getitem__(self,key):
+        answer = []
+        keys = self.rounded_keys(key)
+        for key in keys:
+            if self.__items__.has_key(key):
+                answer.extend(self.__items__[key])
+        #if len(answer)==0:
+        #    raise KeyError#FIXME or return KeyError
+        #                  #FIXME or just return []?
+        else:
+            return answer #FIXME or unique(answer)?
+
+    def __delete__(self,item):
+        keys = self.rounded_keys(item)
+        answer = False
+        #if any of the possible keys contains
+        #a list, return true
+        for key in keys:        
+            if self.__items__.has_key(key):
+                if item in self.__items__[key]:
+                    self.__items__[key].remove(item)
+
+    def remove(self,item):
+        self.__delete__(item)
+
+    def __contains__(self,item):
+
+        keys = self.rounded_keys(item)
+        answer = False
+        #if any of the possible keys contains
+        #a list, return true
+        for key in keys:        
+            if self.__items__.has_key(key):
+                if item in self.__items__[key]:
+                    answer = True
+        return answer
+
+
+    def has_item(self,item):
+        return self.__contains__(item)
+
+    def round_up_rel2(self,value):
+         t_rel=self.__t_rel__
+         #Rounding up the value
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         #m is the mantissa, e the exponent
+         # 0.5 < |m| < 1.0
+         m = m+t_rel*(10**-(self.__p_rel__))
+         #bump m up
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_down_rel2(self,value):
+         t_rel=self.__t_rel__
+         #Rounding down the value
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         #m is the mantissa, e the exponent
+         # 0.5 < m < 1.0
+         m = m-t_rel*(10**-(self.__p_rel__))
+         #bump the |m| down, by 5% or whatever
+         #self.p_rel dictates
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_flat_rel2(self,value):
+    #redundant
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_up_rel(self,value):
+         t_rel=self.__t_rel__
+         #Rounding up the value
+         m,e = self.frexp(value)
+         #m is the mantissa, e the exponent
+         # 0.5 < |m| < 1.0
+         m = m+t_rel*(10**-(self.__p_rel__))
+         #bump m up
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_down_rel(self,value):
+         t_rel=self.__t_rel__
+         #Rounding down the value
+         m,e = self.frexp(value)
+         #m is the mantissa, e the exponent
+         # 0.5 < m < 1.0
+         m = m-t_rel*(10**-(self.__p_rel__))
+         #bump the |m| down, by 5% or whatever
+         #self.p_rel dictates
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_flat_rel(self,value):
+    #redundant
+         m,e = self.frexp(value)
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_up_abs(self,value):
+         t_abs=self.__t_abs__
+         #Rounding up the value
+         m = value+t_abs*(10**-(self.__p_abs__))
+         #bump m up
+         m = round(m,self.__p_abs__)
+         return m
+
+    def round_down_abs(self,value):
+         t_abs=self.__t_abs__
+         #Rounding down the value
+         m = value-t_abs*(10**-(self.__p_abs__))
+         #bump the |m| down, by 5% or whatever
+         #self.p_rel dictates
+         m = round(m,self.__p_abs__)
+         return m
+
+    def round_flat_abs(self,value):
+    #redundant?
+         m = round(value,self.__p_abs__)
+         return m
+
+    def keys(self):
+        return self.__items__.keys()
+
+
+class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
+    """
+    This is a discretised tuple set, but 
+    based on a mapping. The mapping MUST
+    return a sequence.
+
+    example: 
+    def weight(animal):
+        return [animal.weight]
+    
+    a = Mapped_Discretised_Tuple_Set(weight)
+    a.append[cow]
+    a.append[fox]
+    a.append[horse]
+
+    a[horse]    -> [cow,horse]
+    a[dog]      -> [fox]
+    a[elephant] -> []
+    """
+    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
+        Discretised_Tuple_Set.__init__\
+         (self,p_rel,t_rel = t_rel)
+        self.mapping = mapping
+
+    def rounded_keys(self,key):
+        mapped_key = tuple(self.mapping(key))
+        keys = self.__rounded_keys__(mapped_key)
+        return keys
+
+class Affine_Linespace(Mapped_Discretised_Tuple_Set):
+    """
+    The affine linespace creates a way to record and compare lines.
+    Precision is a bit of a hack, but it creates a way to avoid 
+    misses caused by round offs (between two lines of different
+    lenghts, the short one gets rounded off more). 
+    I am starting to think that a quadratic search would be faster.
+    Nearly.
+    """
+    def __init__(self,p_rel=4,t_rel=0.2):
+        Mapped_Discretised_Tuple_Set.__init__\
+            (self,self.affine_line,\
+            p_rel=p_rel,t_rel=t_rel)
+
+        roundings = \
+        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
+        self.roundings = roundings
+        #roundings = \
+        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
+        #self.roundings = roundings
+
+    def affine_line(self,line):
+        point_1 = line[0]
+        point_2 = line[1]
+        #returns the equation of a line
+        #between two points, in the from
+        #(a,b,-c), as in ax+by-c=0
+        #or line *dot* (x,y,1) = (0,0,0)
+
+        #Note that it normalises the line
+        #(a,b,-c) so Line*Line = 1.
+        #This does not change the mathematical
+        #properties, but it makes comparism
+        #easier.
+
+        #There are probably better algorithms.
+        x1 = point_1[0]
+        x2 = point_2[0]
+        y1 = point_1[1]
+        y2 = point_2[1]
+        dif_x = x1-x2
+        dif_y = y1-y2
+
+        if dif_x == dif_y == 0:
+            msg = 'points are the same'
+            raise msg
+        elif abs(dif_x)>=abs(dif_y):
+            alpha = (-dif_y)/dif_x
+            #a = alpha * b
+            b = -1.
+            c = (x1*alpha+x2*alpha+y1+y2)/2.
+            a = alpha*b
+        else:
+            beta = dif_x/(-dif_y)
+            #b = beta * a
+            a = 1.
+            c = (x1+x2+y1*beta+y2*beta)/2.
+            b = beta*a
+        mag = abs(a)+abs(b)
+        #This does not change the mathematical
+        #properties, but it makes comparism possible.
+
+        #note that the gradient is b/a, or (a/b)**-1.
+        #so 
+
+        #if a == 0:
+        #    sign_a = 1.
+        #else:
+        #    sign_a = a/((a**2)**0.5)
+        #if b == 0:
+        #    sign_b = 1.
+        #else:
+        #    sign_b = b/((b**2)**0.5)
+        #if c == 0:
+        #    sign_c = 1.
+        #else:
+        #    sign_c = c/((c**2)**0.5)
+        #a = a/mag*sign_a
+        #b = b/mag*sign_b
+        #c = c/mag*sign_c
+        a = a/mag
+        b = b/mag
+        c = c/mag
+        return a,b,c
