Index: undation/ga/storm_surge/analytical solutions/domain.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/domain.py	(revision 538)
+++ 	(revision )
@@ -1,1311 +1,0 @@
-"""Classes implementing general 2D triangular domains for
-finite-volume computations.
-
-   Copyright 2004
-   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-   Geoscience Australia
-"""
-
-from mesh import Point, Vector, Mesh, Triangle,\
-     rectangular_mesh, circular_mesh, oblique_mesh,\
-     contracting_channel_mesh,\
-     initialise_consecutive_datastructure as mesh_init
-
-from consecutive import Consecutive
-
-from types import ListType 
-
-class Volume(Triangle):
-    """A finite-volume element
-    """
-
-    from config import number_of_conserved_quantities as M
-    arrays = Triangle.arrays
-    arrays['conserved_quantities_centroid'] = ('Float', M)
-    arrays['conserved_quantities_vertex0'] = ('Float', M)
-    arrays['conserved_quantities_vertex1'] = ('Float', M) 
-    arrays['conserved_quantities_vertex2'] = ('Float', M)   
-    arrays['conserved_quantities_face0'] = ('Float', M)
-    arrays['conserved_quantities_face1'] = ('Float', M)
-    arrays['conserved_quantities_face2'] = ('Float', M)
-
-    #Temporary storage for forcing values (incl flux)
-    arrays['explicit_update'] = ('Float', M)
-    arrays['semi_implicit_update'] = ('Float', M)            
-    
-
-    from config import number_of_field_values as M
-    arrays['field_values_centroid'] = ('Float', M)
-    arrays['field_values_vertex0'] = ('Float', M)
-    arrays['field_values_vertex1'] = ('Float', M) 
-    arrays['field_values_vertex2'] = ('Float', M)   
-    arrays['field_values_face0'] = ('Float', M)
-    arrays['field_values_face1'] = ('Float', M)
-    arrays['field_values_face2'] = ('Float', M)
-
-
-    #Indices of volumes to use for gradient computations            
-    arrays['gradient_info'] = ('Int', 3) 
-
-
-    def __init__(self,
-                 vertex0, vertex1, vertex2,
-                 conserved_quantities = None):
-        
-        Triangle.__init__(self, vertex0, vertex1, vertex2)
-
-        if conserved_quantities is not None:
-            from config import number_of_conserved_quantities as M
-            msg = 'Number of conserved quantities is set to %d in config.py' %M
-            assert len(conserved_quantities) == M, msg
-            self.set_conserved_quantities('centroid', conserved_quantities)
-
-
-
-    def __repr__(self):
-        v = []
-        for id in self.vertices:
-            v.append(Point.instances[id])
-        return 'Volume(%s, %s, %s)' %tuple(v)
-
-
-    def get_conserved_quantities(self, where, i=None):
-        if where in ['face', 'vertex']:
-            if i is None:
-                #Return all of them
-                result = []
-                for i in range(3):
-                    cmd = 'q = self.conserved_quantities_%s%d' %(where, i)
-                    exec(cmd)                
-                    result.append(q)
-            else:
-                cmd = 'result = self.conserved_quantities_%s%d' %(where, i)
-                exec(cmd)
-        elif where == 'centroid':
-            result = self.conserved_quantities_centroid
-        else:
-            msg = 'Unknown argument: %s' %where
-            raise ValueError, msg
-
-        return result
-    
-    def get_field_values(self, where, i=None):
-        if where in ['face', 'vertex']:
-            if i is None:
-                #Return all of them
-                result = []
-                for i in range(3):
-                    cmd = 'q = self.field_values_%s%d' %(where, i)
-                    exec(cmd)                
-                    result.append(q)
-            else:
-                cmd = 'result = self.field_values_%s%d' %(where, i)
-                exec(cmd)
-        elif where == 'centroid':
-            result = self.field_values_centroid
-        else:
-            msg = 'Unknown argument: %s' %where
-            raise ValueError, msg
-
-        return result
-
-
-    def set_conserved_quantities(self, where, arg2, arg3=None):
-        if where in ['face', 'vertex']:
-            i = arg2
-            value = arg3            
-            cmd = 'self.conserved_quantities_%s%d[:] = value' %(where, i)
-            exec(cmd)
-        elif where == 'centroid':
-            value = arg2
-            self.conserved_quantities_centroid[:] = value
-        else:
-            msg = 'Unknown argument: %s' %where
-            raise ValueError, msg
-
-
-    def set_field_values(self, where, arg2, arg3=None):
-        if where in ['face', 'vertex']:
-            i = arg2
-            value = arg3                        
-            cmd = 'self.field_values_%s%d[:] = value' %(where, i)
-            exec(cmd)
-        elif where == 'centroid':
-            value = arg2            
-            self.field_values_centroid[:] = value
-        else:
-            msg = 'Unknown argument: %s' %where
-            raise ValueError, msg
-        
-
-    def precompute(self):
-        """Initialise, e.g. precomputed field_values on edges and at centroid
-        Pre-condition: field_values at vertices have been set
-        """
-
-        f0, f1, f2 = self.get_field_values('vertex')
-
-        # FIXME: interpolate field values as a seperate function
-        # as a class function
-        #Update field values at faces
-        self.set_field_values('face', 0, 0.5*(f1 + f2))
-        self.set_field_values('face', 1, 0.5*(f0 + f2))
-        self.set_field_values('face', 2, 0.5*(f0 + f1))
-
-        #Update field values at centroid 
-        self.set_field_values('centroid', (f0 + f1 + f2)/3)
-
-
-        ######################################################
-        # Precompute information needed by computations of gradients
-        # self.gradient_info
-
-        #Find all neighbouring volumes that are not boundaries
-        boundary_index = 0
-        number_of_boundaries = 0
-        for k in range(3):
-            if self.neighbours[k] < 0:
-                boundary_index = k    #Record position of last boundary
-                number_of_boundaries += 1                
-            else:
-                self.gradient_info[k] = self.neighbours[k]
-
-        if self.get_number_of_boundaries() == 1:
-            #We have one boundary - use own centroid
-            #as a surrogate for the one boundary neighbour
-            self.gradient_info[boundary_index] = self.id
-
-        msg = 'Actual number of boundaries does not match recorded number'
-        assert number_of_boundaries == self.get_number_of_boundaries(), msg
-        
-            
-
-
-    def set_boundary(self, face, boundary_object):
-        """Associate boundary object with volume id and face id.
-	Neighbour index is set to point to boundary data  
-	"""
-	
-        # Allow None so internal boundaries in pmesh can be skipped
-        if boundary_object != None:
-            #Create new boundary value entry
-            B = Boundary_value(self.id, face, boundary_object)
-        
-            #Register index where values are to be stored        
-            self.neighbours[face] = -B.id-1  
-        
-
-    def interpolate_conserved_quantities(cls):
-        """Initialise conserved quantities at centroid from vertex quantities.
-        This is usefull if an initial condition is specified at vertices
-        instead of at the centroid
-        """    
-        q0 = cls.conserved_quantities_vertex0
-        q1 = cls.conserved_quantities_vertex1
-        q2 = cls.conserved_quantities_vertex2        
-
-        #Update field values at centroid.
-        cls.conserved_quantities_centroid[:,:] = (q0 + q1 + q2)/3
-
-	#Update edges. Normally, they are taken care of by the first call 
-	#to distribute. However, if we start from a checkpoint, this is needed.
-        cls.conserved_quantities_face0[:,:] = (q1 + q2)/2	
-        cls.conserved_quantities_face1[:,:] = (q0 + q2)/2	
-        cls.conserved_quantities_face2[:,:] = (q0 + q1)/2			
-    
-    interpolate_conserved_quantities = classmethod(interpolate_conserved_quantities) 
-
-class Domain(Mesh):
-    """Representation of a 2D triangulated domain for use
-    with finite_volume computations.
-    Descendant of class Mesh
-    """
-
-    def __repr__(self):
-        return 'Domain(%s)' %self.elements
-    
-    def __init__(self, volumes = None,
-                 flux_function = None,
-                 limiter = None,
-                 gravity = None,
-                 friction = None,                                  
-                 explicit_forcing_terms = None,
-                 semi_implicit_forcing_terms = None,                 
-                 marker_dict = None,
-                 default_order = 2):
-        
-        Mesh.__init__(self, volumes, marker_dict = marker_dict)
-
-        Boundary_value.initialise(self.number_of_boundaries)
-
-        #Assign user specified flux function
-        self.flux_function = flux_function
-
-        #Assign user specified limiter function
-        self.limiter = limiter
-
-        #Assign user specified gravity function
-        self.gravity = gravity
-
-        #Assign user specified friction function
-        self.friction = friction                
-
-        #Spatial order (1 or 2)
-        self.default_order = default_order        
-        self.order = self.default_order
-        self.smallsteps = 0
-        self.number_of_steps = 0
-        self.number_of_first_order_steps = 0        
-
-        #Assign user specified explicit forcing terms if any
-        # 
-        # Explicit must have the form
-        #
-        #    G(q, t)
-        #
-        # and explicit scheme is
-        #
-        #   q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
-	#
-	#
-	# To call it, though, Volume, t must be passed as the arguments
-
-        
-        if explicit_forcing_terms is None:
-            self.explicit_forcing_terms = []
-        else:
-            self.explicit_forcing_terms = explicit_forcing_terms
-
-        #Assign user specified semi implicit forcing terms if any
-        # 
-        # Semi implicit forcing terms are assumed to have the form
-        #
-        #    G(q, t) = H(q, t) q
-        #
-        # and the semi implicit scheme will then be
-        #
-        #   q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
-        
-        if semi_implicit_forcing_terms is None:
-            self.semi_implicit_forcing_terms = []
-        else:
-            self.semi_implicit_forcing_terms = semi_implicit_forcing_terms                        
-
-        #Model time    
-        self.time = 0.0
-        self.finaltime = None
-        self.min_timestep = self.max_timestep = 0.0        
-	
-	
-        #Checkpoint		
-	self.filename = 'domain'
-	self.checkpoint = False
-
-        #Stored output
-        self.store=False
-	self.format = 'dat'    
-        #self.smooth = False
-        self.smooth = True
-
-        #Reduction operation for get_vertex_values              
-        #from pytools.stats import mean
-        #self.reduction = mean 
-        self.reduction = min  #Looks better near steep slopes
-
-        #Class that holds all volume data
-
-        self.volume_class = Volume
-        
-
-    def get_format(self):
-        return self.format
-
-    def get_name(self):
-        return self.filename    
-    
-    def smooth_vertex_values(self, value_array='field_values',
-                             precision = None):
-        """ Smooths field_values or conserved_quantities data.
-        TODO: be able to smooth individual fields
-        NOTE:  This function does not have a test.
-        """
-
-        from Numeric import concatenate, zeros, Float, Int, array, reshape
-
-        
-        A,V = self.get_vertex_values(xy=False,
-                                     value_array=value_array,
-                                     smooth = True,
-                                     precision = precision)
-
-
-        #Set some field values
-        for volume in self:
-            for i,v in enumerate(volume.vertices):
-                if value_array == 'field_values':
-                    volume.set_field_values('vertex', i, A[v,:])
-                elif value_array == 'conserved_quantities':
-                    volume.set_conserved_quantities('vertex', i, A[v,:])
-
-        if value_array == 'field_values':
-            self.precompute()
-        elif value_array == 'conserved_quantities':
-            Volume.interpolate_conserved_quantities()
-        
-    
-    #Method for outputting model results
-    def get_vertex_values(self, xy=True, value_array='field_values',
-                          indices = None,
-                          smooth = None,
-
-                          precision = None,
-                          reduction = None):
-        """Return either field values or conserved quantities
-        at vertices
-
-        indices is a tuple that specifies which of the quantities
-        should be returned.
-
-        if domain.smooth is True data will be smoothed according
-        to domain.smoother (currently hardwired into this method).
-        In this case vertex coordinates will be de-duplicated
-        If no smoothings is required vertex coordinates and values will
-        be aggregated as a concatenation of of values at
-        vertices 0, vertices 1 and vertices 2
-
-        The connectivity is represented as an integer array, V, of dimension
-        M x 3, where M is the number of volumes. Each row has three indices
-        into the X, Y, A arrays defining the triangle.
-
-        
-        Calling convention
-        if xy is True:
-           X,Y,A,V = get_vertex_values
-
-        else:
-           A,V = get_vertex_values         
-           
-        """
-
-        from Numeric import concatenate, zeros, Float, Int, array, reshape
-
-
-        if smooth is None:
-            smooth = self.smooth
-
-        if precision is None:
-            precision = Float
-
-        if reduction is None:
-            reduction = self.reduction
-
-
-        if indices is None:
-
-            #Get number of columns in array
-            cmd = 'indices = tuple(range(Volume.%s_vertex0.shape[1]))'\
-                  %(value_array)
-            exec(cmd)    
-        
-        N = len(indices)
-        if smooth == True:
-            M = len(self.vertexdict) #Number of unique vertices
-            
-            A = zeros((M,N), precision)
-            V = zeros((len(self),3), Int)            
-
-            #Rebuild triangle structure
-            for k, v in enumerate(self.vertexdict.keys()):            
-                for volume_id, vertex_id in self.vertexdict[v]:                
-                    V[volume_id, vertex_id] = k
-
-
-            #Smoothing loop
-
-            for j, i in enumerate(indices):
-                #For each quantity
-                #
-                #j is an index into returned quantities (0 <= j < N)
-                #i the given index into internal quantities
-                #(e.g. 0: bed elevation or 1: friction)
-
-                for k, v in enumerate(self.vertexdict.keys()):
-                    #Go through all contributions to vertex v
-                    #k is an index into returned vertices (0 <= k < M)
-                    #v is an index into internal vertices
-                    
-                    L = len(self.vertexdict[v])
-
-                    ulist = []
-                    for volume_id, vertex_id in self.vertexdict[v]:
-                        #For all contributing triange, vertex pairs
-                        
-                        cmd = 'u = Volume.%s_vertex%d[%d,%d]'\
-                              %(value_array, vertex_id, volume_id, i)
-                        exec(cmd)
-                        ulist.append(u)
-
-                    A[k,j] = reduction(ulist)    
-
-
-            if xy is True:
-                X = zeros(M, precision)
-                Y = zeros(M, precision)
-
-                for k, v in enumerate(self.vertexdict.keys()):            
-                    X[k], Y[k] = Point.coordinates[v].astype(precision)
-
-                return X, Y, A, V
-            else:
-                return A, V
-        else:
-            #Don't smooth
-
-            m = len(self)  #Number of volumes
-            M = 3*m        #Total number of unique vertices
-
-            A = zeros((M,N), precision)
-
-            #Create definition of triangles
-            tmp = reshape(array(range(m)).astype(Int), (m,1))
-            V = concatenate( (tmp, tmp+m, tmp+2*m), axis = 1 )
-
-            #Get values
-            cmd = 'B = concatenate('
-            cmd += '(Volume.%s_vertex0, Volume.%s_vertex1, Volume.%s_vertex2)'\
-                   %( (value_array,)*3)
-            cmd += ', axis=0)'
-            exec(cmd)
-
-            #Extract desired indices
-            for j, i in enumerate(indices):            
-                A[:,j] = B[:,i].astype(precision) 
-
-            #Do vertex coordinates    
-            if xy is True:                
-                X0, Y0, X1, Y1, X2, Y2 = self.get_vertex_coordinates()
-                X = concatenate((X0, X1, X2), axis=0).astype(precision)
-                Y = concatenate((Y0, Y1, Y2), axis=0).astype(precision)
-                
-                return X, Y, A , V           
-            else:
-                return A, V
-
-
-    def store_xya(self):
-        """Save x,y and bed elevation
-	"""
-
-        import data_manager
-        #Initialise write
-
-
-        self.writer = data_manager.get_dataobject(self.filename,
-                                                  'xya',
-                                                  len(self),
-                                                  mode = 'w')
-
-        X0, Y0, X1, Y1, X2, Y2 = self.get_vertex_coordinates() 
-        self.writer.store_all(X0, X1, X2, Y0, Y1, Y2,
-                              Volume.field_values_vertex0,
-                              Volume.field_values_vertex1,
-                              Volume.field_values_vertex2,
-                              Volume.conserved_quantities_vertex0,
-                              Volume.conserved_quantities_vertex1,
-                              Volume.conserved_quantities_vertex2)
-    	
-    def store_bathymetry(self):
-        """Save x,y and bed elevation
-	"""
-
-        import data_manager
-        
-        #Initialise write
-        self.writer = data_manager.get_dataobject(self, mode = 'w')
-
-        #Store
-        self.writer.store_all()
-
-
-    def store_conserved_quantities(self):
-        """Save conserved quantities e.g. for visualisation
-        using specified format
-	"""
-
-        self.writer.store_timestep()
-        
-	
-	
-    def store_checkpoint(self):
-        """Checkpoint: 
-	   save conserved quantities in given format under given filename  
-	"""
-        #FIXME: Maybe use generic data format as in netcdf
-
-        format = 'cpt'
-        import data_manager
-        writer = data_manager.get_variable_writer(format)
-        
-        writer(self.filename, self.time,
-               Volume.conserved_quantities_vertex0,
-               Volume.conserved_quantities_vertex1,
-               Volume.conserved_quantities_vertex2)
-        
-	
-
-    def get_stored_checkpoints(self):
-        """get stored checkpoints
-	"""
-
-        format = 'cpt'
-        
-        from config import data_dir
-        
-        import glob, os
-	
-	files = glob.glob(data_dir+os.sep+'*.'+format)
-    
-	times = []
-        for file in files:
-
-            #Skip checkpoint of size is wrong
-            i0 = file.find('_size')
-            if i0 >= 0:
-                i0 += 5  #Position where size starts
-                i1 = file.find('_time') #Find where size ends
-                if i1 > i0:
-                    size = int(file[i0:i1])
-                    if size != len(self):
-                        continue
-                    
-                
-                
-            
-	    i0 = file.find('_time')
-            if i0 == -1:
-                #Skip bathymetry file
-                continue
-            
-	    i0 += 5  #Position where time starts           
-	    i1 = file.find('.'+format)	    
-
-	    if i1 > i0:
-	        t = float(file[i0:i1])
-                times.append(t)
-                
-        times.sort()
-        return times
-
-	    
-    def goto_latest_checkpoint(self):
-        """Checkpoint: 
-	if a checkpoint exists read conserved quantities from it.
-	If latest checkpoint cannot be read, try second latest.
-        In fact - this routine never takes the latest checkpoint in case
-        data wasn't written properly.
-	"""
-        import data_manager
-        reader = data_manager.get_variable_reader('cpt')
-
-        times = self.get_stored_checkpoints()
-
-        i = -2 #Always start at second latest point
-        # Try checkpoints starting with the latest until
-        # one is found (the latest may be corrupted due to an interruption)
-	found = False
-        while i >= -len(times):
-            self.time = times[i]
-
-            try:
-                reader(self.filename, self.time,
-                       Volume.conserved_quantities_vertex0,
-                       Volume.conserved_quantities_vertex1,
-                       Volume.conserved_quantities_vertex2)   
-            except:
-                print 'Could not read checkpoint t=%.4f, trying earlier'\
-                      %self.time
-            else:
-                print 'Proceeding from checkpoint t=%.4f' %self.time 
-                Volume.interpolate_conserved_quantities()                    
-		found = True
-
-                #Get handle to existing netcdf file
-                if self.store:
-                    self.writer = data_manager.get_dataobject(self.filename,
-                                                              self.format,
-                                                              len(self), 'a')
-                
-                break
-
-            i -= 1
-
-	if not found:
-	    print 'Could not find a checkpoint, starting from time = 0.0'
-	    self.time = 0.0
-        
-            
-    def delete_old_checkpoints(self):
-        """Delete all old checkpoints except the three latest
-        FIXME: Only consider checkpoints identified by self.filename
-        FIXME: Later on, perhaps retain certain points as specified
-        in a list of times
-        """
-
-        import os
-        from data_manager import create_filename
-        
-        times = self.get_stored_checkpoints()
-        keep = 4 #How many checkpoints to keep
-	
-        M = Volume.number_of_instances
-        for t in times[:-keep]:
-            FN = create_filename(self.filename, 'cpt', M, t)            
-            try:
-                os.remove(FN)
-            except:
-                print 'WARNING: Checkpoint %s could not be removed' %FN
-                
-        
-    
-    def precompute(self):
-        """Initialise domain, e.g. interpolated field_values
-        and check neighbour structure
-        Must be called after field values have been set
-        """
-
-        for volume in self:
-            volume.precompute()
-
-        self.check_integrity()
-
-
-    def set_field_values(self, *args):
-
-        from Numeric import array
-        M, N = Volume.field_values_vertex0.shape
-
-        if len(args) > N:
-            msg = 'Number of functions exceeds number of field values'
-            raise msg
-        
-
-        for v in self:
-            x0,y0,x1,y1,x2,y2 = v.get_instance_vertex_coordinates()
-
-            for i, func in enumerate(args):
-                if func is not None:
-                    z = func(array([x0,x1,x2]), array([y0,y1,y2])) 
-                
-                    v.field_values_vertex0[i] = z[0]
-                    v.field_values_vertex1[i] = z[1]
-                    v.field_values_vertex2[i] = z[2]                
-
-        #Interpolate vertex values to faces and centroids        
-        self.precompute()
-
-
-
-    def set_conserved_quantities(self, *args):
-
-        from Numeric import array
-        M, N = Volume.conserved_quantities_vertex0.shape
-
-        if len(args) > N:
-            msg = 'Number of function exceeds number of conserved_quantities'
-            raise msg
-        
-
-        for v in self:
-            x0,y0,x1,y1,x2,y2 = v.get_instance_vertex_coordinates()
-            
-            for i, func in enumerate(args):
-                if func is not None:
-                    q = func(array([x0,x1,x2]), array([y0,y1,y2]))
-                    
-                    v.conserved_quantities_vertex0[i] = q[0]
-                    v.conserved_quantities_vertex1[i] = q[1]
-                    v.conserved_quantities_vertex2[i] = q[2]
-
-
-        Volume.interpolate_conserved_quantities()
-        
-    def get_vertex_coordinates(self):
-        """Get all vertex coordinates arranged in the same order as volumes
-        """
-
-        return Volume.get_vertex_coordinates()
-        #FIXME: perhaps this thing should be pre-computed    
-
-
-
-    def distribute_to_vertices_and_edges(self):
-        distribute_to_vertices_and_edges(self)        
-        
-
-    def evolve(self, max_timestep = 0.1, yieldstep = None,
-               finaltime = None):
-        """Evolve model from time=0.0 to finaltime using given
-        max_timestep and yielding results every yieldstep.
-
-        The model is determined by the supplied flux_function
-        (self.flux_function) which must be defined as follows:
-
-          def flux_function(ql, qr, normal):
-
-        where ql and qr are the conserved quantities inside and outside
-        each volume (for each of its neighbours) and normal is the
-        corresponding outward normal vector for each neighbour. 
-
-        Evolve is implemented as a generator and is to be called as such, e.g.
-        
-        for t in domain.evolve(timestep, yieldstep, finaltime):
-            <Do something with domain and t>
-        
-        """
-
-        import data_manager
-        from config import min_timestep, epsilon
-
-        if yieldstep is None:
-            yieldstep = max_timestep
-
-
-        self.order = self.default_order
-        
-            
-        self.yieldtime = 0.0 #Time between 'yields'
-
-        #Initialise interval of timestep sizes (for reporting only)
-        self.min_timestep = max_timestep
-        self.max_timestep = min_timestep
-        self.finaltime = finaltime
-        self.number_of_steps = 0
-        self.number_of_first_order_steps = 0                
-
-        #Initial update of vertex and edge values
-        distribute_to_vertices_and_edges(self)
-	
-	
-        #Or maybe restore from latest checkpoint
-	if self.checkpoint is True:
-            self.goto_latest_checkpoint()
-            
-
-        #Store model data, e.g. for visualisation    
-        if self.store is True and self.time == 0.0:
-            self.store_bathymetry()        
-            self.store_conserved_quantities()	
-
-	
-        yield(self.time)  #Yield initial values
-	
-        while True:
-            #Update boundary values
-            update_boundary_values()            
-
-            #Compute all fluxes and timestep suitable for all volumes
-            compute_flux(self, max_timestep)
-
-            #Update timestep to fit yieldstep and finaltime
-            update_timestep(self, yieldstep, finaltime)
-
-            #Update conserved quantities
-            update_conserved_quantities(self)            
-
-            #Update vertex and edge values
-            distribute_to_vertices_and_edges(self)
-                
-            #Update time    
-            self.time += self.timestep
-            self.yieldtime += self.timestep
-            self.number_of_steps += 1
-            if self.order == 1:
-                self.number_of_first_order_steps += 1        
-
-            #Yield results
-            if finaltime is not None and abs(self.time - finaltime) < epsilon:
-                # Yield final time and allow inspection of domain
-                yield(self.time)
-                break
-                
-            if abs(self.yieldtime - yieldstep) < epsilon: 
-                # Yield (intermediate) time and allow inspection of domain
-
-                if self.checkpoint is True:
-                    self.store_checkpoint()
-                    self.delete_old_checkpoints()
-                    
-                    #Store model data, e.g. for visualisation    
-                if self.store is True:
-                    self.store_conserved_quantities()
-
-                    
-                yield(self.time)  
-
-                # Reinitialise
-                self.yieldtime = 0.0                
-                self.min_timestep = max_timestep
-                self.max_timestep = min_timestep
-                self.number_of_steps = 0
-                self.number_of_first_order_steps = 0                        
-            
-    def evolve_to_end(self, max_timestep = 0.1, finaltime = 1.0):
-        """Iterate evolve generator all the way to the end
-        """
-
-        for _ in self.evolve(max_timestep, None, finaltime):            
-            pass
-        
-
-
-    def set_boundary(self, tags):
-        """Given a dictionary of boundaries, keyed by markers,
-        add boundaries to all faces with key markers.
-        Pre-condition: self.boundary has been built.
-
-        NOTE: If a marker is not used to the domain, no error is thrown.
-        FIXME: This woul lead to implementation of a
-        default boundary condition
-        """
-
-        for id, face in self.boundary:
-            volume = Volume.instances[id]
-            
-            marker = self.boundary[ (id, face) ]
-
-            if tags.has_key(marker):
-                volume.set_boundary(face, tags[marker])                
-            else:
-                msg = 'Unknown boundary tag: %s' %marker
-                raise msg
-
-
-            
-    def check_integrity(self):
-        """Check integrity of domain
-        """
-        
-        Mesh.check_integrity(self, Element=Volume)
-
-        assert len(Boundary_value.conserved_quantities) ==\
-               self.number_of_boundaries
-
-
-
-    def write_time(self):
-        if self.min_timestep == self.max_timestep:
-            print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-                  %(self.time, self.min_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-        else:    
-            print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-                  %(self.time, self.min_timestep, self.max_timestep, self.number_of_steps,
-                    self.number_of_first_order_steps)
-
-            
-
-
-    def pmesh_to_domain(cls, fileName, tags = None, setting_function = None):
-        import sys
-        from config import pmesh_filename
-        sys.path.append(pmesh_filename)
-        from load_mesh.loadASCII import import_trianglulation
-
-        try:
-            meshdic = import_trianglulation(fileName)
-        except IOError, e:        
-            msg = 'Could not open file %s ' %fileName
-            raise IOError, msg 
-
-        domain  = Domain.pmesh_dictionary_to_domain(meshdic,
-                                              setting_function = setting_function)
-        if tags:
-            domain.set_boundary(tags)
-        return domain
-
-    def pmesh_dictionary_to_domain(cls, meshdic, setting_function = None):
-        """
-        convert a pmesh dictionary to a list of Volumes.
-        Also, return a list of triangles which have boundary tags
-        meshdic structure;
-        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)  
-        generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
-        generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
-        generated segment marker list: [S1Marker, S2Marker, ...] (list of ints)
-        triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
-        triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
-        triangle attribute list: [T1att, T2att, ...] (list of floats)
-        """
-
-        E = len(meshdic['generatedtrianglelist'])       #Number of triangular elements
-        P = len(meshdic['generatedpointlist'])  #Number of initial vertices
-        #print "E", E
-        #print "P",P
-
-        mesh_init(P+E,E,Point, Volume, Domain)
-
-        mesh_vertices = []
-
-        for index, point in enumerate(meshdic['generatedpointlist']):
-            v=Point(point[0], point[1])
-            mesh_vertices.append(v)
-            v.pmesh_id = index
-
-        
-        volumes = []
-        for triangle in meshdic['generatedtrianglelist']:
-            tObject = Volume( mesh_vertices[triangle[0]],
-                          mesh_vertices[triangle[1]],
-                          mesh_vertices[triangle[2]] )
-            volumes.append(tObject)
-        if setting_function:
-            if not(type(setting_function) is ListType):
-                setting_function = [setting_function]
-            for funct in setting_function:
-                meshdic = funct(meshdic, vertices = mesh_vertices,
-                                volumes = volumes)
-        
-        marker_dict = pmesh_dict_to_marker_dict(volumes,meshdic,mesh_vertices)
-    
-        domain = Domain(volumes = volumes,
-                        marker_dict = marker_dict)
-        #DSG!!! TODO smoothing here is a hack and should be removed.
-        # Can the setting functions be applied after the domain has
-        # been instanciated?
-        #domain.smooth_vertex_values(value_array="field_values")
-        #domain.smooth_vertex_values(value_array="conserved_quantities")
-        domain.precompute()
-        return domain
-
-    #specify pmesh_to_domain as a class function
-    pmesh_to_domain = classmethod(pmesh_to_domain)
-    pmesh_dictionary_to_domain = classmethod(pmesh_dictionary_to_domain)
-
-
-class Boundary_value(Consecutive):
-
-    from config import number_of_conserved_quantities as M
-    arrays = {'conserved_quantities': ('Float', M)}
-
-    from config import number_of_field_values as M    
-    arrays['field_values'] = ('Float', M)
-
-
-    def __init__(self, vol, face, boundary_object):
-        Consecutive.__init__(self)
-        
-        self.volume = vol
-        self.face = face
-        self.boundary_object = boundary_object
-
-
-    def update(self):
-        #FIXME: Maybe only update field_values and dirichlet's the first time
-
-        face = self.face
-        
-        V = Volume.instances[self.volume]
-        q = self.boundary_object.get_conserved_quantities(V, face) #Get value
-
-        #z = self.boundary_object.get_field_values(V, face) #Get field values
-        #Since boundary always returns field values from calling face
-        #we might as well just take those :-)
-        if face == 0:
-            z = V.field_values_face0
-        elif face == 1:
-            z = V.field_values_face1
-        else:
-            z = V.field_values_face2            
-
-        self.conserved_quantities[:] = q
-        self.field_values[:] = z        
-                
-
-        
-    
-#Module functions
-def set_initialcondition(domain, q):
-    """Set conserved quantities to q
-    for all volumes. Argument q must be a list of numbers.
-    """
-
-    from Numeric import array
-    for volume in domain:
-        volume.conserved_quantities = array(q)            
-        
-
-    #FIXME: Some precomputation could be entered here    
-    #E.g. distribute_to_vertices_and_edges(domain, order=2)
-
-
-
-def rectangular_domain(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0), 
-                       symmetry_type='diagonal'):
-    """Setup a rectangular domain of volumes.
-       See mesh.rectangular_mesh for details.
-    """
-    return rectangular_mesh(m, n, len1, len2, origin, Point, Volume, Domain, 
-                            symmetry_type)
-
-
-def circular_domain(m, n, radius=1.0, center = (0.0, 0.0)):    
-    """Setup a circular domain of volumes.
-       See mesh.circular_mesh for details.
-    """
-    return circular_mesh(m, n, radius, center, Point, Volume, Domain)
-
-
-def oblique_domain(m, n, lenx=1.0, leny=1.0, theta=8.95, origin = (0.0, 0.0)):   
-    """Setup a oblique domain of volumes.
-       See mesh.oblique_mesh for details.
-    """
-    return oblique_mesh(m, n, lenx, leny, theta, origin, Point, Volume, Domain)
-
-
-def contracting_channel_domain(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0,
-                               y1 = 0.0, y4 = 1./4., y8 = 1.0,
-                               origin = (0.0, 0.0)):
-    """Setup a contracting_channel domain of volumes.
-       See mesh.contracting_channel_mesh for details.
-    """
-    return contracting_channel_mesh(m, n, x1, x2, x3, x4, y1, y4, y8,
-                                    origin, Point, Volume, Domain)
-
-def Strang_domain(filename):
-    """Read Strang generated mesh.
-    """
-    from mesh import Strang_mesh
-    return Strang_mesh(filename, Point, Volume, Domain)
-
-   
-def pmesh_dict_to_marker_dict(triangles,meshdic,Vertices):
-    """ Convert the pmesh dictionary (meshdic) description of boundary tags
-    to a dictionary of markers, indexed with volume id and face number.
-    """
-    sides = calc_sides(triangles)
-    marker_dict = {}
-    for seg,marker in map(None,meshdic['generatedsegmentlist'],
-                          meshdic['generatedsegmentmarkerlist']):
-        v1 = seg[0]
-        v2 = seg[1]
-        for key in [(v1,v2),(v2,v1)]:
-            if sides.has_key(key):
-                #this creates a dict of lists of faces, indexed by marker
-                #tagged_edges.setdefault(marker,[]).append(sides[key])
-                marker_dict[sides[key]] = marker
-    return marker_dict
-    
-def calc_sides(triangles):
-    #Build dictionary mapping from sides (2-tuple of points)
-    #to left hand side neighbouring triangle        
-    sides = {}
-    for triangle in triangles:
-        a = triangle.vertices[0]
-        b = triangle.vertices[1]
-        c = triangle.vertices[2]       
-        sides[a,b] = (triangle.id, 2) #(id, face)
-        sides[b,c] = (triangle.id, 0) #(id, face)
-        sides[c,a] = (triangle.id, 1) #(id, face)
-    return sides        
-
-
-
-def initialise_consecutive_datastructure(points=0, elements=0,
-                                         point_class=Point, element_class=Volume,
-                                         mesh_class=Domain):
-
-    mesh_init(points, elements,
-              point_class, element_class, mesh_class)    
-
-
-
-############################
-# Main components of evolve
-############################
-
-#FIXME: consider passing only domain into C extensions.
-        
-
-
-def update_conserved_quantities_c(domain):        
-
-
-    from domain_ext import update_conserved_quantities
-    update_conserved_quantities(domain.time,
-                                domain.timestep,
-                                Volume,
-                                Point.coordinates,
-                                domain.gravity,
-                                domain.friction,
-                                domain.explicit_forcing_terms,
-                                domain.semi_implicit_forcing_terms)
-
-
-
-
-def distribute_to_vertices_and_edges_c(domain):
-    """Extrapolate conserved quantities from centroid to
-    vertices and edge-midpoints for each volume
-    This version works directly with consecutive data structure
-    """
-
-    from domain_ext import distribute_to_vertices_and_edges as func
-    func(domain.limiter,
-         Volume.instances,          #only used if limiter is specified 
-         Volume.number_of_instances,
-         Volume.number_of_boundaries,
-         Volume.neighbours,         
-         Volume.field_values_centroid,
-         Volume.field_values_vertex0,
-         Volume.field_values_vertex1,
-         Volume.field_values_vertex2,
-         Volume.conserved_quantities_centroid,
-         Volume.conserved_quantities_vertex0,
-         Volume.conserved_quantities_vertex1,
-         Volume.conserved_quantities_vertex2,
-         Volume.conserved_quantities_face0,
-         Volume.conserved_quantities_face1,
-         Volume.conserved_quantities_face2,
-         Volume.gradient_info,
-         Volume.points,
-         Point.coordinates,
-         domain.order)
-
-
-def update_boundary_values():
-    """Go through list of boundary objects and update values
-    """
-
-    for B in Boundary_value.instances:
-        if B is not None:
-            B.update()
-        #else:
-        #    raise 'Boundary was None ?'
-            
-
-def compute_flux_c(domain, max_timestep):
-    """Compute all fluxes and the timestep suitable for all volumes
-    This version works directly with consecutive data structure
-    and calls C-extension.
-    Note that flux function is hardwired into C-extension.
-    """
-    
-    import sys
-
-    neighbours = Volume.neighbours
-    neighbour_faces = Volume.neighbour_faces
-    normals = Volume.normals
-    flux = Volume.explicit_update    
-
-    area = Volume.geometric[:,0]            
-    radius = Volume.geometric[:,1]            
-    edgelengths = Volume.geometric[:,2:]
-
-
-    flux[:] = 0.0 #Reset stored fluxes to zero        
-    from domain_ext import compute_flux
-    timestep = compute_flux(domain.flux_function,
-                            Volume.conserved_quantities_face0,
-                            Volume.conserved_quantities_face1,
-                            Volume.conserved_quantities_face2,
-                            Volume.field_values_face0,
-                            Volume.field_values_face1,
-                            Volume.field_values_face2,
-                            Boundary_value.conserved_quantities,
-                            Boundary_value.field_values,
-                            Vector.coordinates,
-                            neighbours,
-                            neighbour_faces,
-                            normals,
-                            flux, area, radius, edgelengths,
-                            max_timestep)
-
-    domain.timestep = timestep
-
-
-
-
-def update_timestep(domain, yieldstep, finaltime):
-
-    from config import min_timestep, max_smallsteps
-    
-    timestep = domain.timestep
-    
-    #Record maximal and minimal values of timestep for reporting    
-    domain.max_timestep = max(timestep, domain.max_timestep)
-    domain.min_timestep = min(timestep, domain.min_timestep)
-    
-    #Protect against degenerate time steps
-    if timestep < min_timestep:
-        domain.smallsteps += 1 #Number of consecutive small steps taken b4 taking action
-
-        if domain.smallsteps > max_smallsteps:
-            domain.smallsteps = 0 #Reset
-           
-            if domain.order == 1:
-                msg = 'Minimal timestep %.16f reached using 1 order scheme'\
-                      %domain.min_timestep  
-                raise msg
-            else:
-                #Try to overcome situation by switching to 1 order
-                #print 'Switching temporarily to first order: t = %.8f, delta t=%.8f'\
-                #      %(domain.time, timestep)
-                domain.order = 1
-
-    else:
-        domain.smallsteps = 0
-        if domain.order == 1:
-            if domain.order != domain.default_order:
-                #print 'Switching back to second order: t = %.8f, delta t=%.8f: '\
-                #      %(domain.time, timestep)
-                domain.order = 2
-                
-
-
-    #Ensure that final time is not exceeded
-    if finaltime is not None and domain.time + timestep > finaltime:
-        timestep = finaltime-domain.time
-
-    #Ensure that model time is aligned with yieldsteps
-    if domain.yieldtime + timestep > yieldstep:
-        timestep = yieldstep-domain.yieldtime
-
-    domain.timestep = timestep
-
-    
-
-##############################################
-#Initialise module
-
-#C-extensions
-import util
-if util.can_use_C_extension('domain_ext.c'):
-    compute_flux = compute_flux_c
-    distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
-    update_conserved_quantities = update_conserved_quantities_c
-else:
-    from python_versions import compute_flux
-    from python_versions import distribute_to_vertices_and_edges
-    from python_versions import update_conserved_quantities
-
-
-#Optimisation with psyco
-from config import use_psyco
-if use_psyco:
-    try:
-        import psyco
-    except:
-        msg = 'WARNING: psyco (speedup) could not import'+\
-              ', you may want to consider installing it'
-        print msg
-    else:
-        psyco.bind(distribute_to_vertices_and_edges)
-        psyco.bind(compute_flux)  
-        psyco.bind(update_boundary_values) 
-        #psyco.bind(Domain.update_timestep)     #Not worth it
-        psyco.bind(update_conserved_quantities)
-        
-        psyco.bind(rectangular_domain)
-
-if __name__ == "__main__":
-    pass
Index: undation/ga/storm_surge/analytical solutions/mesh.py
===================================================================
--- /inundation/ga/storm_surge/analytical solutions/mesh.py	(revision 538)
+++ 	(revision )
@@ -1,1094 +1,0 @@
-"""Classes implementing general 2D geometrical mesh of triangles.
-
-   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-   Geoscience Australia, 2004   
-"""
-
-
-from consecutive import Consecutive        
-
-class XY:
-    """A generic coordinate set in the 2D plane
-    """
-
-    def __init__(self, x, y):
-        from Numeric import array
-        self.coordinates = array([x, y])
-
-    def get_x(self):
-        return self.coordinates[0]
-
-    def get_y(self):
-        return self.coordinates[1]
-
-    def set_x(self, x):
-        self.coordinates[0] = x
-
-    def set_y(self, x):
-        self.coordinates[1] = y
-
-    def set_coordinates(self, x, y):
-        self.coordinates[0] = x
-        self.coordinates[1] = y        
-
-    def distance_to(self, other):
-        """Returns the Euclidian distance from this point to another
-        """
-
-        from math import sqrt
-        return sqrt( (self.coordinates[0] - other.coordinates[0])**2 +
-                     (self.coordinates[1] - other.coordinates[1])**2 )
-
-
-
-
-class Point(Consecutive, XY):
-    """A generic point in the 2D plane with consecutive storage 
-    """
-
-    arrays = {'coordinates': ('Float', 2)}
-    
-    def __init__(self, x, y, transient = False):
-        Consecutive.__init__(self, transient)
-
-        #Store coordinates 
-        self.coordinates[0] = x
-        self.coordinates[1] = y        
-
-    def __repr__(self):
-        return 'Point(%s, %s)' %(self.coordinates[0], self.coordinates[1])
-
-
-class Vector(Point):
-    """A vector in the plane
-    """
-
-    def __init__(self, x, y, transient = False):
-        Point.__init__(self, x, y, transient)
-        
-
-    def __repr__(self):
-        return 'Vector(%s, %s)' %(self.coordinates[0], self.coordinates[1])
-
-
-    def length(self):
-        from math import sqrt
-        return sqrt( self.coordinates[0]**2 + self.coordinates[1]**2 )        
-
-    def __div__(self, other):
-        """Divides x and y component by other. If other is not a
-        number an exception is raised
-        """
-
-        try:
-            self.coordinates[0] /= other
-            self.coordinates[1] /= other
-        except:
-            raise TypeError, 'Vectors must be divided by scalars only'
-
-        return self
-
-    def __sub__(self, other):
-        """Forms difference vector between self and other.
-        """
-
-        return Vector (self.coordinates[0] - other.coordinates[0],
-                       self.coordinates[1] - other.coordinates[1])
-
-    def normalise(self):
-        """Divide x and y component by length.
-        """
-
-        l = self.length()
-
-        self.coordinates[0] /= l
-        self.coordinates[1] /= l
-
-
-
-    def inner_product(self, other):
-        """Forms the dot product between self and other.
-        """
-
-        return self.coordinates[0]*other.coordinates[0] +\
-               self.coordinates[1]*other.coordinates[1]
-
-    def angle(self):
-        """Compute angle between e1 (the unit vector in the x-direction)
-        and this vector
-        """
-
-        from math import acos, pi
-
-        v1 = self.coordinates[0]/self.length()
-        v2 = self.coordinates[1]/self.length()
-        
-        theta = acos(v1)
-
-        if v2 < 0:
-            #Quadrant 3 or 4
-            theta = 2*pi-theta
-
-        return theta
-
-
-    def anglediff(self, other):
-        """Compute difference between angle of vector and angle of other.
-        This is used for determining the ordering of vertices,
-        e.g. for checking if they are counter clockwise. 
-
-        Always return a positive value
-        """
-        
-        from math import pi
-
-        a0 = other.angle()
-        a1 = self.angle()
-
-        #Ensure that difference will be positive
-        if a1 < a0:
-            a1 += 2*pi
-            
-        return a1-a0
-
-
-
-class Triangle(Consecutive):
-    """A triangle defined in terms of three ve 
-       vertices assumed to be
-       ordered counter-clockwise.
-       Triangle inherits from class Point which represents
-       its centroid and provide the method distance_to.
-       Triangle is purely geometric.
-    """
-
-    arrays = {
-        #Integer array of indices into 
-        #vertices and centroid
-        'points': ('Int', 4),
-	
-        #Integer array of indices into normal vectors
-        'normals': ('Int', 3),        
-        
-        #Integer array of indices into 
-        #neighbours and associated edge numbers 
-        'neighbours': ('Int', 3),
-        'neighbour_faces': ('Int', 3),
-        'number_of_boundaries': ('Int', 1),                
-
-        #Float array of triangle data 
-        #area, radius, edgelength0 - 2
-        'geometric': ('Float', 5)
-        }
-
-    
-    def __init__(self, vertex0, vertex1, vertex2, transient = False):
-
-        Consecutive.__init__(self, transient)        
-
-        #Record vertex ids in common data array
-        self.vertices = self.points[:3]
-        self.vertices[0] = vertex0.id
-        self.vertices[1] = vertex1.id
-        self.vertices[2] = vertex2.id
-        
-        #Temporary variables
-        x0, y0, x1, y1, x2, y2 = self.get_instance_vertex_coordinates()
-       
-        #Compute centroid
-        x = (x0 + x1 + x2)/3
-        y = (y0 + y1 + y2)/3
-        centroid = Point(x, y, transient)
-        self.points[3] = centroid.id
-        self.centroid = self.points[3:] #FIXME: remove
-        
-        #Area
-        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-        self.geometric[0] = area
-        
-
-        msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
-        msg += ' is degenerate:  area == %f' %area
-        assert area > 0.0, msg
-        
-        #Midpoints        
-        x = (x1 + x2)/2
-        y = (y1 + y2)/2
-        m0 = Point(x,y,transient=True)
-
-        x = (x0 + x2)/2
-        y = (y0 + y2)/2
-        m1 = Point(x,y,transient=True)
-
-        x = (x1 + x0)/2
-        y = (y1 + y0)/2
-        m2 = Point(x,y,transient=True)
-
-        #The radius is the distance from the centroid of
-        #a triangle to the midpoint of the side of the triangle
-        #closest to the centroid
-        d0 = centroid.distance_to(m0)
-        d1 = centroid.distance_to(m1)        
-        d2 = centroid.distance_to(m2)
-        
-        radius = min(d0, d1, d2)
-        self.geometric[1] = radius
-
-        #Normals
-        #The normal vectors
-        # - point outward from each edge
-        # - are orthogonal to the edge
-        # - have unit length
-        # - Are enumerated according to the opposite corner:
-        #   (First normal is associated with the edge opposite
-        #    the first vertex, etc)
-
-        x = x2 - x1
-        y = y2 - y1
-        n0 = Vector(y, -x, transient)
-        l0 = n0.length()
-        n0.normalise()
-
-        x = x0 - x2
-        y = y0 - y2
-        n1 = Vector(y, -x, transient)
-        l1 = n1.length()
-        n1.normalise()        
-
-        x = x1 - x0
-        y = y1 - y0
-        n2 = Vector(y, -x, transient)
-        l2 = n2.length()
-        n2.normalise()                
-
-        self.normals[0] = n0.id
-        self.normals[1] = n1.id
-        self.normals[2] = n2.id
-        #print 'N', self.normals, 'V', self.vertices, 'C', self.centroid
-
-        #Segmentlengths
-        self.edgelengths = self.geometric[2:5]        
-        self.edgelengths[0] = l0
-        self.edgelengths[1] = l1
-        self.edgelengths[2] = l2        
-
-        #Initialise Neighbours (-1 means that it is a boundary neighbour)
-        self.neighbours[0] = -1
-        self.neighbours[1] = -1
-        self.neighbours[2] = -1
-
-        #Initialise face ids of neighbours
-        #In case of boundaries this slot is not used
-        self.neighbour_faces[0] = -1
-        self.neighbour_faces[1] = -1
-        self.neighbour_faces[2] = -1                       
-        
-
-    def get_vertex_coordinates(cls):
-        """Convert vertex information to arrays of vertex coordinates for
-        all triangles: X0, Y0, X1, Y1, X2, Y2
-        """
-
-        from Numeric import zeros, Float
-
-        M = cls.number_of_instances
-        
-        X0 = zeros(M, Float)  #x of vertex 0 of all volumes
-        X1 = zeros(M, Float)  #x of vertex 1 of all volumes
-        X2 = zeros(M, Float)  #x of vertex 2 of all volumes        
-
-        Y0 = zeros(M, Float)  #x of vertex 0 of all volumes
-        Y1 = zeros(M, Float)  #x of vertex 1 of all volumes
-        Y2 = zeros(M, Float)  #x of vertex 2 of all volumes        
-
-        for i in range(M):
-            k = cls.points[i,0]  #Index of vertex 0
-            X0[i] = Point.coordinates[k,0]
-            Y0[i] = Point.coordinates[k,1]
-
-            k = cls.points[i,1]  #Index of vertex 0
-            X1[i] = Point.coordinates[k,0]
-            Y1[i] = Point.coordinates[k,1]
-
-            k = cls.points[i,2]  #Index of vertex 0
-            X2[i] = Point.coordinates[k,0]
-            Y2[i] = Point.coordinates[k,1]                        
-            
-        return X0, Y0, X1, Y1, X2, Y2            
-    
-    get_vertex_coordinates = classmethod(get_vertex_coordinates)
-
-    
-    def get_instance_vertex_coordinates(self):
-        """Convert vertex information to x0,y0,x1,y1,x2,y2
-        for current instance
-        """
-
-        coords = []
-        for i in range(3):
-            k = self.vertices[i] #Vertex index
-            x = Point.coordinates[k,0]
-            y = Point.coordinates[k,1]            
-            
-            coords.append(x)
-            coords.append(y)            
-
-        return coords
-
-    def get_vertex_objects(self):
-        vertexobjs = []
-        for i in range(3):
-            k = self.vertices[i] #Vertex index
-            v = Point.instances[k]
-            vertexobjs.append(v)
-        
-        return vertexobjs
-
-    def get_normals(self, face=None):
-        if face is None:
-            normals = []
-            for k in range(3):
-                normals.append( Vector.instances[self.normals[k]] )
-
-            return normals
-        else:
-            return Vector.instances[self.normals[face]]           
-            
-        
-    def get_radius(self):
-        return self.geometric[1]
-    
-    def get_area(self):
-        return self.geometric[0]
-
-    def calc_area(self):
-        #Area
-        x0, y0, x1, y1, x2, y2 = self.get_instance_vertex_coordinates()
-        return abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
-
-    def get_centroid(self):
-        if self.transient:
-            raise 'get_centroid not implemented for transient points'
-        else:
-            return Point.instances[self.points[3]]        
-        
-    def set_number_of_boundaries(self, n):
-        self.number_of_boundaries[0] = n
-
-    def get_number_of_boundaries(self):
-        return self.number_of_boundaries[0]
-
-
-    def check_integrity(self):
-        """Check that triangle is internally consistent e.g.
-        that area corresponds to edgelengths, that vertices
-        are arranged in a counter-clockwise order, etc etc
-        """
-        
-        from config import epsilon
-        from math import pi
-
-        #Temporary variables
-        x0, y0, x1, y1, x2, y2 = self.get_instance_vertex_coordinates()
-        area = self.geometric[0]
-
-        #Check area
-        msg = '%s, ' %str(self.vertices)
-        msg += 'vertex coordinates: %f %f %f %f %f %f' %(x0,y0,x1,y1,x2,y2)
-        assert abs((area - self.calc_area())/area) < epsilon, msg
-        
-      
-        #Check that points are arranged in counter clock-wise order
-        v0 = Vector(x1-x0, y1-y0, transient = True)
-        v1 = Vector(x2-x1, y2-y1, transient = True)
-        v2 = Vector(x0-x2, y0-y2, transient = True)
-
-        a0 = v1.anglediff(v0)
-        a1 = v2.anglediff(v1)
-        a2 = v0.anglediff(v2)        
-
-        msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
-        in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
-        assert a0 < pi and a1 < pi and a2 < pi, msg
-        
-        #Check that normals are orthogonal to edge vectors
-        #Note that normal[k] lies opposite vertex k
-
-
-        normal0, normal1, normal2 = self.get_normals()
-        #normal0 = Vector.instances[self.normals[0]]
-        #normal1 = Vector.instances[self.normals[1]]
-        #normal2 = Vector.instances[self.normals[2]]
-
-        #print v0, normal2
-        assert abs(v0.inner_product(normal2))/v0.length() < epsilon
-        #print v1, normal0        
-        assert abs(v1.inner_product(normal0))/v1.length() < epsilon
-        #print v2, normal1
-        assert abs(v2.inner_product(normal1))/v2.length() < epsilon        
-
-
-        
-    def __repr__(self):
-        v = []
-        for id in self.vertices:
-            v.append(Point.instances[id])
-        return 'Triangle(%s, %s, %s)' %tuple(v)
-
-
-class Mesh:
-    """Representation of a 2D triangular mesh.
-    """
-    
-    def __init__(self, elements = None,  marker_dict = None):
-        
-        if elements is None:
-            self.elements = []
-            self.vertexdict = {} #Dictionary of all vertices to be built
-                                 #when neighbour structure is set up
-        
-        else:
-            self.elements = elements
-            self.vertexdict, self.number_of_boundaries =\
-               self.build_neighbour_and_boundary(marker_dict = marker_dict)
-
-
-    def __repr__(self):
-        return 'Mesh(%s)' %self.elements
-    
-    def __len__(self):
-        return len(self.elements)
-
-    def __getitem__(self, i):
-        return self.elements[i]
-
-
-    def build_neighbour_structure(self):
-        """Update all registered triangles to point to their neighbours.
-        Also, build a list of all vertices
-        """
-        #Step 1:
-        #Build dictionary mapping from segments (2-tuple of points)
-        #to left hand side face (facing neighbouring triangle)
-        #Also build dictionary of vertices
-        vertexdict = {}
-        neighbour_dict = {} 
-        for triangle in self.elements:
-            a = triangle.vertices[0]
-            b = triangle.vertices[1]
-            c = triangle.vertices[2]            
-            neighbour_dict[a,b] = (triangle.id, 2) #(id, face)
-            neighbour_dict[b,c] = (triangle.id, 0) #(id, face)
-            neighbour_dict[c,a] = (triangle.id, 1) #(id, face)
-
-            #Register the vertices as keys mapping to 
-            #(triangle, face) tuples associated with them
-            #This is used for smoothing
-            for vertex_id, v in enumerate([a,b,c]):
-                if not vertexdict.has_key(v):
-                    vertexdict[v] = []
-                    
-                vertexdict[v].append( (triangle.id, vertex_id) )
-
-        #Step 2:
-        #Go through triangles again, but this time
-        #reverse direction of segments and lookup neighbours.
-        for triangle in self.elements:            
-            a = triangle.vertices[0]
-            b = triangle.vertices[1]
-            c = triangle.vertices[2]
-
-            number_of_boundaries = 3
-            if neighbour_dict.has_key((b,a)):
-                triangle.neighbours[2] = neighbour_dict[b,a][0]
-                triangle.neighbour_faces[2] = neighbour_dict[b,a][1]
-                number_of_boundaries -= 1
-                
-            if neighbour_dict.has_key((c,b)):
-                triangle.neighbours[0] = neighbour_dict[c,b][0]
-                triangle.neighbour_faces[0] = neighbour_dict[c,b][1]
-                number_of_boundaries -= 1               
-                
-            if neighbour_dict.has_key((a,c)):
-                triangle.neighbours[1] = neighbour_dict[a,c][0]
-                triangle.neighbour_faces[1] = neighbour_dict[a,c][1]
-                number_of_boundaries -= 1              
-                
-            triangle.set_number_of_boundaries(number_of_boundaries)
-
-        return vertexdict
-    
-    def build_marker_dictionary(self,  marker_dict = None):
-        """Build or input the marker dictionary.
-         self.boundary is a dictionary of markers,
-         keyed by volume id and face.
-         Postcondition: self.boundary is defined.
-        """
-        if  marker_dict is None:
-            boundary = {}
-            #self.boundary = {}
-            for triangle in self.elements:
-                for i in range(0,3):
-                    if triangle.neighbours[i] == -1:
-                        boundary[(triangle.id, i)] = 0 #Default tag 
-            self.boundary = boundary
-        else:
-            self.boundary =  marker_dict
-          
-    def build_neighbour_and_boundary(self,  marker_dict = None):
-        vertexdict = self.build_neighbour_structure()
-        self.build_marker_dictionary( marker_dict =  marker_dict)
-        total_number_of_boundaries = len(self.boundary)
-        return vertexdict, total_number_of_boundaries
-
-    def get_boundarytags(self):
-        """Return list of available boundary tags
-        """
-        
-        tags = {}
-        for v in self.boundary.values():
-            tags[v] = 1
-
-        return tags.keys()    
-            
-
-    def check_integrity(self, Element=Triangle):
-        """Check integrity of3 mesh
-        """
-
-        #Check triangle integrity separately
-        for triangle in self:
-            triangle.check_integrity()
-
-
-        #Check mesh integrity
-        for triangle in self:
-            for v in triangle.vertices:
-                #Check that all vertices have been registered
-                assert self.vertexdict.has_key(v)
-
-                #Check that this triangle is listed with at least one vertex
-                assert (triangle.id, 0) in self.vertexdict[v] or\
-                       (triangle.id, 1) in self.vertexdict[v] or\
-                       (triangle.id, 2) in self.vertexdict[v]
-                
-                
-
-            #Check neighbour structure
-            for i, id in enumerate(triangle.neighbours):
-                #Boundaries need not fulfill this
-                if id >= 0:
-                    neighbour = Element.instances[id]
-                    face = triangle.neighbour_faces[i]
-
-                    assert neighbour.neighbours[face] == triangle.id
-
-
-        #Check integrity of boundary list
-        assert len(self.boundary) == self.number_of_boundaries
-
-        #NOTE: This assert doesn't hold true if ther are internal boundaries
-        #FIXME: Look into this further.
-        #for id, face in self.boundary:
-        #    assert Element.neighbours[id][face] < 0
-        
-
-
-# Module functions and constants
-def rectangular_mesh(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0),
-                     point_class=Point, element_class=Triangle, 
-		     mesh_class=Mesh, symmetry_type = 'diagonal'):
-    """Setup a rectangular grid of triangles
-    with m+1 by n+1 grid points
-    and side lengths len1, len2. If side lengths are omitted
-    the mesh defaults to the unit square.
-    
-    len1: x direction (left to right)
-    len2: y direction (bottom to top)
-
-    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
-      rectangular_mesh(...,Triangle=Volume, Mesh=Domain)
-    """
-
-    #FIXME: symmetry_type only works for diagonal    
-    
-    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)
-    delta1 = float(len1)/m
-    delta2 = float(len2)/n    
-
-    #Dictionary of vertex objects
-    vertices = {}
-    for i in range(m+1):
-        for j in range(n+1):
-            vertices[i,j] = Point(i*delta1 + origin[0], j*delta2 + origin[1])
-
-    #Construct 2 triangles per rectangular 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]            
-	    
-            if symmetry_type == 'diagonal':
-        	elements.append(element_class(v4,v3,v2)) #Lower                
-        	elements.append(element_class(v1,v2,v3)) #Upper    
-            else:
-	        if j < n//2:
-                    #print 'Low', i, j
-                    elements.append(element_class(v4,v3,v2)) #Lower
-        	    elements.append(element_class(v1,v2,v3)) #Upper    	    
-                else:  
-                    #print 'Hi', i, j		
-                    elements.append(element_class(v2,v4,v1)) #Lower
-        	    elements.append(element_class(v3,v1,v4)) #Upper		
-
-    M = mesh_class(elements)
-
-    #Set a default tagging
-    for id, face in M.boundary:
-    
-	
-        e = element_class.instances[id]
-        x0, y0, x1, y1, x2, y2 = e.get_instance_vertex_coordinates()
-
-        if symmetry_type == 'diagonal':
-            if face==2:
-                #Left or right#
-                if abs(x0-origin[0]) < epsilon:
-                    M.boundary[(id,face)] = 'left'            
-                elif abs(origin[0]+len1-x0) < epsilon:
-                    M.boundary[(id,face)] = 'right'
-                else:
-                    raise 'Unknown boundary'                            
-            elif face==1:
-                #Top or bottom
-                if abs(y0-origin[1]) < epsilon:            
-                    M.boundary[(id,face)] = 'bottom'
-                elif abs(origin[1]+len2-y0) < epsilon:
-                    M.boundary[(id,face)] = 'top'
-                else:
-                    raise 'Unknown boundary'  
-            else:
-                print face, id, id%m, m, n
-                raise 'Unknown boundary'          
-        else:
-	
-            #Identify lower half of rectangle
-	    if id//n % 2 == 0:
-	        #Lower half 
-            	if face==2:
-            	    #Left or right#
-            	    if abs(x0-origin[0]) < epsilon:
-            	        M.boundary[(id,face)] = 'left'            
-            	    elif abs(origin[0]+len1-x0) < epsilon:
-            	        M.boundary[(id,face)] = 'right'
-            	    else:
-		        print id, face, x0
-            	        raise 'Unknown boundary'                            
-            	elif face==1:
-            	    #Top or bottom
-            	    if abs(y0-origin[1]) < epsilon:            
-            	        M.boundary[(id,face)] = 'bottom'
-            	    elif abs(origin[1]+len2-y0) < epsilon:
-            	        M.boundary[(id,face)] = 'top'
-            	    else:
-		        print id, face, x0		    
-            	        raise 'Unknown boundary'  
-            	else:
-            	    print face, id, id%m, m, n
-            	    raise 'Unknown boundary'          
-	    else:    
-                #Upper half 
-            	if face==1:
-            	    #Left or right#
-            	    if abs(x0-origin[0]) < epsilon:
-            	        M.boundary[(id,face)] = 'left'            
-            	    elif abs(origin[0]+len1-x0) < epsilon:
-            	        M.boundary[(id,face)] = 'right'
-            	    else:
-		        print id, face, x0		    
-            	        raise 'Unknown boundary'                            
-            	elif face==2:
-            	    #Top or bottom
-            	    if abs(y0-origin[1]) < epsilon:            
-            	        M.boundary[(id,face)] = 'bottom'
-            	    elif abs(origin[1]+len2-y0) < epsilon:
-            	        M.boundary[(id,face)] = 'top'
-            	    else:
-		        print id, face, x0		    
-            	        raise 'Unknown boundary'  
-            	else:
-            	    print face, id, id%m, m, n
-            	    raise 'Unknown boundary'          		       
-    return M
-
-def Strang_mesh(filename,point_class=Point, element_class=Triangle, 
-		     mesh_class=Mesh):
-    """Read Strang generated mesh.
-    """
-
-    from math import pi
-    
-    fid = open(filename)
-    points = [] # empty list
-    triangles = [] #empty list
-    
-    for line in fid.readlines():
-#        print line #testing only
-        fields = line.split()
-        if len(fields) == 2:
-            # we are reading vertex coordinates
-            points.append((float(fields[0]),float(fields[1])))
-        elif len(fields) == 3:
-            # we are reading tiangle point id's
-            triangles.append((int(float(fields[0]))-1,
-                              int(float(fields[1]))-1,
-                              int(float(fields[2]))-1))
-        else:
-            raise 'wrong format in ' + filename
-#    print len(points),len(triangles)
-
-    E = len(triangles)        #Number of triangular elements
-    P = len(points)  #Number of initial vertices
-
-    initialise_consecutive_datastructure(P+E, E,
-                                         point_class,
-                                         element_class,
-                                         mesh_class)                                      
-    vertices = []
-    elements = []
-
-    for p in points:
-        vertices.append(point_class(p[0],p[1]))
-
-    for t in triangles:
-        v0 = vertices[t[0]]
-        v1 = vertices[t[1]]
-        v2 = vertices[t[2]]
-        x0 = v0.get_x()
-        y0 = v0.get_y()
-        x1 = v1.get_x()
-        y1 = v1.get_y()
-        x2 = v2.get_x()
-        y2 = v2.get_y()
-        
-#Check that points are arranged in counter clock-wise order
-        vec0 = Vector(x1-x0, y1-y0, transient = True)
-        vec1 = Vector(x2-x1, y2-y1, transient = True)
-        vec2 = Vector(x0-x2, y0-y2, transient = True)
-
-        a0 = vec1.anglediff(vec0)
-        a1 = vec2.anglediff(vec1)
-        a2 = vec0.anglediff(vec2)        
-
-        if a0 < pi and a1 < pi and a2 < pi:
-            elements.append(element_class(v0,v1,v2))
-        else:
-            elements.append(element_class(v0,v2,v1))       
-                        
-    M = mesh_class(elements)
-    return M
-
-#FIXME: These haven't been done yet
-def circular_mesh(m, n, radius=1.0, center = (0.0, 0.0),  Triangle=Triangle,
-                  Mesh=Mesh, Point=Point):
-    """Setup a circular grid of triangles
-    with m+1 radial segments by n+1 points around the circuference
-    and radius. If radius is are omitted the mesh defaults to the unit circle
-    radius.
- 
-    radius: radius of circle
-
-    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
-    rectangular_mesh(...,Triangle=Volume, Mesh=Domain)"""
-    
-
-    from Numeric import array
-    from visual import rate
-    import math
-
-    delta = radius/float(m)
-
-    #Dictionary of vertex objects
-    vertices = {}
-    for i in range(n+1):
-        theta = float(i)*(2*math.pi)/float(n)  
-        for j in range(m+1):
-            delta = float(j)*radius/float(m)
-            vertices[i,j] = Point(delta*math.cos(theta),delta*math.sin(theta))
-
-    #Construct 2 triangles per element        
-    elements = []
-    for i in range(n):
-        for j in range(1,m):
-            v1 = vertices[i,j+1]
-            v2 = vertices[i,j]            
-            v3 = vertices[i+1,j+1]            
-            v4 = vertices[i+1,j]
-
-            elements.append(Triangle(v4,v2,v3)) #Lower                
-            elements.append(Triangle(v1,v3,v2)) #Upper    
-
-    #Do the center
-    v1 = vertices[0,0]
-    for i in range(n):
-        v2 = vertices[i,1]
-        v3 = vertices[i+1,1]
-
-        T = Triangle(v1,v2,v3) #center
-        elements.append(T)
-
-    return Mesh(elements)
-
-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 an 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
-                
-            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
-
-    for id, face in M.boundary:
-
-        e = element_class.instances[id]
-        x0, y0, x1, y1, x2, y2 = e.get_instance_vertex_coordinates()
-
-        if face==2:
-            #Left or right#
-            if abs(x0-origin[0]) < epsilon:
-                M.boundary[(id,face)] = 'left'            
-            elif abs(origin[0]+lenx-x0) < epsilon:
-                M.boundary[(id,face)] = 'right'
-            else:
-                print face, id, id%m, m, n
-                raise 'Left or Right Unknown boundary'                            
-        elif face==1:
-            #Top or bottom
-            if x0 > 0.25*lenx and abs(y0-a1-a2*x0-origin[1]) < epsilon or\
-               x0 <= 0.25*lenx and abs(y0-origin[1]) < epsilon:            
-                   M.boundary[(id,face)] = 'bottom'
-            elif abs(origin[1]+leny-y0) < epsilon:
-                M.boundary[(id,face)] = 'top'
-            else:
-                print face, id, id%m, m, n
-                raise 'Top or Bottom Unknown boundary'  
-        else:
-            print face, id, id%m, m, n
-            raise 'Unknown boundary'          
-    
-    return M
-
-def contracting_channel_mesh(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0,
-                            y1 =0.0, y4 = -1./4., y8 = 1.0,
-                             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 = (x4 - x1)/float(m)
-    deltay = (y8 - y1)/float(n)
-    a = y4 - y1
-    
-    if y8 - a <= y1 + a:
-        print a,y1,y4
-        raise 'Constriction is too large reduce a'
-            
-    y2 = y1
-    y3 = y4
-    x5 = x4
-    y5 = y8 - a
-    x6 = x3
-    y6 = y5
-    x7 = x2
-    y7 = y8
-    x8 = x1
-    
-    a2 = a/(x3 - x2)
-    a1 = a - a*x3/(x3 - x2)
-    a4 = (-a + a2*(x7 - x6))/(x6 - x7)/y7
-    a3 = (y7 - a1 - x7*a2 - a4*x7*y7)/y7
-    
-    # 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 > x2:
-                if x < x3:
-                    y = a1 + a2*x + a3*y + a4*x*y
-                else:
-                    y = a + y*(y5 - y4)/(y8 - y1)   
-                
-            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
-
-    for id, face in M.boundary:
-
-        e = element_class.instances[id]
-        x_0, y_0, x_1, y_1, x_2, y_2 = e.get_instance_vertex_coordinates()
-        lenx = x4 - x1
-        if face==2:
-            #Left or right#
-            if abs(x_0-origin[0]) < epsilon:
-                M.boundary[(id,face)] = 'left'            
-            elif abs(origin[0]+lenx-x_0) < epsilon:
-                M.boundary[(id,face)] = 'right'
-            else:
-                print face, id, id%m, m, n
-                raise 'Left or Right Unknown boundary'                            
-        elif face==1:
-            #Top or bottom
-            if x_0 <= x2 and abs(y_0-y1-origin[1]) < epsilon or\
-               x_0 > x3 and abs(y_0-y3-origin[1]) < epsilon or\
-               x_0 >  x2 and x_0 <= x3 and abs(y_0-(y2+(y3-y2)*(x_0-x2)/(x3-x2)+origin[1])) < epsilon:            
-                M.boundary[(id,face)] = 'bottom'
-            elif x_0 <= x7 and abs(y_0-y8-origin[1]) < epsilon or\
-               x_0 > x6 and abs(y_0-y6-origin[1]) < epsilon or\
-               x_0 > x7 and x_0 <= x6 and abs(y_0-(y7+(y6-y7)*(x_0-x7)/(x6-x7)+origin[1])) < epsilon:            
-                M.boundary[(id,face)] = 'top'
-            else:
-                print face, id, id%m, m, n 
-                raise 'Top or Bottom Unknown boundary'  
-        else:
-            print face, id, id%m, m, n
-            raise 'Unknown boundary'          
-
-            
-      #  print id, face, M.boundary[(id,face)],x_0,y_0,x_1,y_1,x_2,y_2
-    
-    return M
-
-
-
-#Map from edge number to indices of associated vertices
-edge_map = ((1,2), (0,2), (0,1))
-
-                                         
-def initialise_consecutive_datastructure(points=0, elements=0,
-                                         point_class=Point, 
-					 element_class=Triangle,
-                                         mesh_class=Mesh):
-    """Initialise consecutive datastructures for 
-    vertices, volumes and vectors
-    """
-    
-    point_class.initialise(points)
-    element_class.initialise(elements)
-    Vector.initialise(3*elements)
-                                         
-    
-
