Changeset 1387

Timestamp:

May 13, 2005, 6:15:08 PM (19 years ago)

Author:

steve

Message:

Need to look at uint test for inscribed circle

Location:

inundation/ga/storm_surge

Files:

• parallel/parallel_advection.py (modified) (2 diffs)
• parallel/run_advection.py (modified) (1 diff)
• parallel/run_parallel_advection.py (modified) (2 diffs)
• parallel/small.tsh (modified) (3 diffs)
• pyvolution/domain.py (modified) (1 diff)
• pyvolution/general_mesh.py (modified) (21 diffs)
• pyvolution/mesh.py (modified) (1 diff)
• pyvolution/mesh_factory.py (modified) (23 diffs)
• pyvolution/netherlands.py (modified) (6 diffs)
• pyvolution/realtime_visualisation_new.py (modified) (10 diffs)
• pyvolution/shallow_water.py (modified) (1 diff)
• pyvolution/shallow_water_ext.c (modified) (40 diffs)
• pyvolution/util.py (modified) (76 diffs)
• zeus/anuga-workspace.zwi (modified) (1 diff)
• zeus/parallel.zpi (modified) (1 diff)
• zeus/pyvolution.zbi (modified) (previous)
• zeus/steve.zpi (modified) (1 diff)

inundation/ga/storm_surge/parallel/parallel_advection.py

@@ -14 +14 @@
 """
 
-#Subversion keywords:
-#
-#$LastChangedDate:  $
-#$LastChangedRevision:  $
-#$LastChangedBy: $
+from advection import *
+Advection_Domain = Domain
+from Numeric import zeros, Float, Int
 
-class Parallel_Domain(Domain):
+class Parallel_Domain(Advection_Domain):
 
-    def __init__(self, coordinates, vertices, boundary = None, velocity = None):
+    def __init__(self, coordinates, vertices, boundary = None, velocity = None,
+                 processor = 0, global_ids = None):
 
+        Advection_Domain.__init__(self, coordinates, vertices, boundary, velocity)
 
-        Generic_domain.__init__(self, coordinates, vertices, boundary,
-                                ['stage'])
+        self.processor = processor
 
-        if velocity is None:
-            self.velocity = [1,0]
+        N = self.number_of_elements
+
+        if  global_ids == None:
+            self.global_ids = global_ids
         else:
-            self.velocity = velocity
-
-        #Only first is implemented for advection
-        self.default_order = self.order = 1
-
-        #Realtime visualisation
-        self.visualise = False
-        self.visualise_color_stage = False
-        self.visualise_timer = True
-        self.visualise_range_z = None
-        self.smooth = True
-
+            self.global_ids = zeros(N, Int)
 
 
     def check_integrity(self):
-        Generic_domain.check_integrity(self)
+        Advection_Domain.check_integrity(self)
 
-        msg = 'Conserved quantity must be "stage"'
+        msg = 'Will need to check global and local numbering'
         assert self.conserved_quantities[0] == 'stage', msg
 
-
-    def flux_function(self, normal, ql, qr, zl=None, zr=None):
-        """Compute outward flux as inner product between velocity
-        vector v=(v_1, v_2) and normal vector n.
-
-        if <n,v> > 0 flux direction is outward bound and its magnitude is
-        determined by the quantity inside volume: ql.
-        Otherwise it is inbound and magnitude is determined by the
-        quantity outside the volume: qr.
-        """
-
-        v1 = self.velocity[0]
-        v2 = self.velocity[1]
-
-
-        normal_velocity = v1*normal[0] + v2*normal[1]
-
-        if normal_velocity < 0:
-            flux = qr * normal_velocity
-        else:
-            flux = ql * normal_velocity
-
-        max_speed = abs(normal_velocity)
-        return flux, max_speed
-
-
-    def compute_fluxes(self):
-        """Compute all fluxes and the timestep suitable for all volumes
-        in domain.
-
-        Compute total flux for each conserved quantity using "flux_function"
-
-        Fluxes across each edge are scaled by edgelengths and summed up
-        Resulting flux is then scaled by area and stored in
-        domain.explicit_update
-
-        The maximal allowable speed computed by the flux_function
-        for each volume
-        is converted to a timestep that must not be exceeded. The minimum of
-        those is computed as the next overall timestep.
-
-        Post conditions:
-        domain.explicit_update is reset to computed flux values
-        domain.timestep is set to the largest step satisfying all volumes.
-        """
-
-        import sys
-        from Numeric import zeros, Float
-        from config import max_timestep
-
-        N = self.number_of_elements
-
-        neighbours = self.neighbours
-        neighbour_edges = self.neighbour_edges
-        normals = self.normals
-
-        areas = self.areas
-        radii = self.radii
-        edgelengths = self.edgelengths
-
-        timestep = max_timestep #FIXME: Get rid of this
-
-        #Shortcuts
-        Stage = self.quantities['stage']
-
-        #Arrays
-        stage = Stage.edge_values
-
-        stage_bdry = Stage.boundary_values
-
-        flux = zeros(1, Float) #Work array for summing up fluxes
-
-        #Loop
-        for k in range(N):
-            optimal_timestep = float(sys.maxint)
-
-            flux[:] = 0.  #Reset work array
-            for i in range(3):
-                #Quantities inside volume facing neighbour i
-                ql = stage[k, i]
-
-                #Quantities at neighbour on nearest face
-                n = neighbours[k,i]
-                if n < 0:
-                    m = -n-1 #Convert neg flag to index
-                    qr = stage_bdry[m]
-                else:
-                    m = neighbour_edges[k,i]
-                    qr = stage[n, m]
-
-
-                #Outward pointing normal vector
-                normal = normals[k, 2*i:2*i+2]
-
-                #Flux computation using provided function
-                edgeflux, max_speed = self.flux_function(normal, ql, qr)
-                flux -= edgeflux * edgelengths[k,i]
-
-                #Update optimal_timestep
-                try:
-                    optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
-                except ZeroDivisionError:
-                    pass
-
-            #Normalise by area and store for when all conserved
-            #quantities get updated
-            flux /= areas[k]
-            Stage.explicit_update[k] = flux[0]
-
-            timestep = min(timestep, optimal_timestep)
-
-        self.timestep = timestep
 
 
@@ -173 +51 @@
         if self.visualise is True and self.time == 0.0:
             import realtime_visualisation_new as visualise
-            visualise.create_surface(self)
+            self.visualiser = visualise.Visualiser(self)
 
         #Call basic machinery from parent class
-        for t in Generic_domain.evolve(self, yieldstep, finaltime):
+        for t in Advection_Domain.evolve(self, yieldstep, finaltime):
             #Real time viz
             if self.visualise is True:
-                visualise.update(self)
+                self.visualiser.update_quantity('stage')
 
             #Pass control on to outer loop for more specific actions
             yield(t)
+
+
+
+def rectangular_with_ghosts(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+    """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)
+
+    Also returns a list of
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements)
+    """
+
+    from config import epsilon
+
+    #E = m*n*2        #Number of triangular elements
+    #P = (m+1)*(n+1)  #Number of initial vertices
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n
+
+    #Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        for j in range(n+1):
+            vertices[i,j] = len(points)
+            points.append([i*delta1 + origin[0], j*delta2 + origin[1]])
+
+
+
+    #Construct 2 triangles per rectangular element and assign tags to boundary
+    elements = []
+    boundary = {}
+    for i in range(m):
+        for j in range(n):
+            v1 = vertices[i,j+1]
+            v2 = vertices[i,j]
+            v3 = vertices[i+1,j+1]
+            v4 = vertices[i+1,j]
+
+            #Update boundary dictionary and create elements
+            if i == m-1:
+                boundary[(len(elements), 2)] = 'right'
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v3,v2]) #Lower element
+
+            if i == 0:
+                boundary[(len(elements), 2)] = 'left'
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v2,v3]) #Upper element
+
+    ghosts = {}
+    i=0
+    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]
+        ghosts.append(elements.index([
+
+    return points, elements, boundary

inundation/ga/storm_surge/parallel/run_advection.py

@@ -19 +19 @@
 
 domain.visualise = True
-#domain.visualise_range_z = 1.0
-domain.visualise_color_stage = True
+domain.visualise_range_z = 0.5
+#domain.visualise_color_stage = True
 
 

inundation/ga/storm_surge/parallel/run_parallel_advection.py

@@ -6 +6 @@
 from config import g, epsilon
 from Numeric import allclose, array, zeros, ones, Float
-from advection import *
+from parallel_advection import *
 from Numeric import array
 
 from mesh_factory import rectangular
 
-points, vertices, boundary = rectangular(60, 60)
+points, vertices, boundary = rectangular(30, 30)
 
 #Create advection domain with direction (1,-1)
-domain = Domain(points, vertices, boundary, velocity=[1.0, 0.0])
+domain = Parallel_Domain(points, vertices, boundary, velocity=[1.0, 0.0])
 
 # Initial condition is zero by default
@@ -29 +29 @@
 
 
-domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
+domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T} )
 domain.check_integrity()
+
+class Set_Stage:
+    """Set an initial condition with constant water height, for x<x0
+    """
+
+    def __init__(self, x0=0.25, x1=0.5, h=1.0):
+        self.x0 = x0
+        self.x1 = x1
+        self.h  = h
+
+    def __call__(self, x, y):
+        return self.h*((x>self.x0)&(x<self.x1))
+
+domain.set_quantity('stage', Set_Stage(0.2,0.4,1.0))
 
 #Check that the boundary value gets propagated to all elements

inundation/ga/storm_surge/parallel/small.tsh

@@ -1 +1 @@
 31 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices...
-0 554.0 -242.0 
-1 443.0 -190.0 
-2 340.0 -184.0 
-3 276.0 -237.0 
-4 281.0 -347.0 
-5 282.0 -380.0 
-6 312.0 -425.0 
-7 387.0 -472.0 
-8 465.0 -487.0 
-9 549.0 -469.0 
-10 566.0 -413.0 
-11 589.0 -359.0 
-12 472.0 -258.0 
-13 377.0 -411.0 
-14 428.0 -257.0 
-15 386.0 -278.0 
-16 372.0 -325.0 
-17 369.0 -371.0 
-18 433.0 -422.0 
-19 476.0 -361.0 
-20 502.0 -304.0 
-21 333.0 -225.0 
-22 305.0 -244.0 
-23 333.0 -270.0 
-24 369.0 -237.0 
-25 278.5 -292.0 
-26 499.90712817 -423.516449623 
-27 532.621862616 -366.885237781 
-28 498.5 -216.0 
-29 322.359287054 -318.872505543 
-30 571.5 -300.5 
+0 554.0 -242.0
+1 443.0 -190.0
+2 340.0 -184.0
+3 276.0 -237.0
+4 281.0 -347.0
+5 282.0 -380.0
+6 312.0 -425.0
+7 387.0 -472.0
+8 465.0 -487.0
+9 549.0 -469.0
+10 566.0 -413.0
+11 589.0 -359.0
+12 472.0 -258.0
+13 377.0 -411.0
+14 428.0 -257.0
+15 386.0 -278.0
+16 372.0 -325.0
+17 369.0 -371.0
+18 433.0 -422.0
+19 476.0 -361.0
+20 502.0 -304.0
+21 333.0 -225.0
+22 305.0 -244.0
+23 333.0 -270.0
+24 369.0 -237.0
+25 278.5 -292.0
+26 499.90712817 -423.516449623
+27 532.621862616 -366.885237781
+28 498.5 -216.0
+29 322.359287054 -318.872505543
+30 571.5 -300.5
 # attribute column titles ...Triangulation Vertex Titles...
 38 # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles...
-0 29 5 17 1 5 26  
-1 17 5 6 -1 3 0  
-2 10 11 27 28 21 -1  
-3 13 17 6 1 6 -1  
-4 15 23 16 24 -1 16  
-5 16 29 17 0 -1 24  
-6 6 7 13 37 3 -1  
-7 29 23 25 9 8 24  
-8 25 4 29 26 7 -1  
-9 25 23 22 10 18 7  
-10 22 23 21 13 11 9  
-11 3 22 21 10 15 18  
-12 21 24 2 14 15 13  
-13 24 21 23 10 16 12  
-14 2 24 1 27 -1 12  
-15 21 2 3 -1 11 12  
-16 24 23 15 4 30 13  
-17 27 30 20 31 36 28  
-18 22 3 25 -1 9 11  
-19 7 8 18 35 37 -1  
-20 8 9 26 23 35 -1  
-21 26 10 27 2 25 23  
-22 26 19 18 -1 35 25  
-23 10 26 9 20 -1 21  
-24 16 23 29 7 5 4  
-25 26 27 19 36 22 21  
-26 29 4 5 -1 0 8  
-27 14 1 24 14 30 29  
-28 27 11 30 -1 17 2  
-29 1 14 12 -1 33 27  
-30 14 24 15 16 -1 27  
-31 20 30 0 -1 34 17  
-32 12 0 28 -1 33 34  
-33 12 28 1 -1 29 32  
-34 0 12 20 -1 31 32  
-35 8 26 18 22 19 20  
-36 27 20 19 -1 25 17  
-37 18 13 7 6 19 -1  
+0 29 5 17 1 5 26 building
+1 17 5 6 -1 3 0
+2 10 11 27 28 21 -1
+3 13 17 6 1 6 -1
+4 15 23 16 24 -1 16
+5 16 29 17 0 -1 24
+6 6 7 13 37 3 -1
+7 29 23 25 9 8 24
+8 25 4 29 26 7 -1
+9 25 23 22 10 18 7
+10 22 23 21 13 11 9 building
+11 3 22 21 10 15 18
+12 21 24 2 14 15 13
+13 24 21 23 10 16 12
+14 2 24 1 27 -1 12
+15 21 2 3 -1 11 12
+16 24 23 15 4 30 13
+17 27 30 20 31 36 28 house
+18 22 3 25 -1 9 11
+19 7 8 18 35 37 -1
+20 8 9 26 23 35 -1
+21 26 10 27 2 25 23
+22 26 19 18 -1 35 25 building
+23 10 26 9 20 -1 21
+24 16 23 29 7 5 4
+25 26 27 19 36 22 21
+26 29 4 5 -1 0 8
+27 14 1 24 14 30 29 house
+28 27 11 30 -1 17 2
+29 1 14 12 -1 33 27 oval
+30 14 24 15 16 -1 27
+31 20 30 0 -1 34 17
+32 12 0 28 -1 33 34
+33 12 28 1 -1 29 32
+34 0 12 20 -1 31 32
+35 8 26 18 22 19 20
+36 27 20 19 -1 25 17
+37 18 13 7 6 19 -1
 28 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments...
 0 14 12 exterior
@@ -81 +81 @@
 7 19 20 exterior
 8 12 20 exterior
-9 21 22 
-10 23 22 
-11 23 24 
-12 21 24 
+9 21 22
+10 23 22
+11 23 24
+12 21 24
 13 4 25 exterior
 14 7 6 exterior
@@ -101 +101 @@
 27 30 11 exterior
 25 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices...
-0 554.0 -242.0 
-1 443.0 -190.0 
-2 340.0 -184.0 
-3 276.0 -237.0 
-4 281.0 -347.0 
-5 282.0 -380.0 
-6 312.0 -425.0 
-7 387.0 -472.0 
-8 465.0 -487.0 
-9 549.0 -469.0 
-10 566.0 -413.0 
-11 589.0 -359.0 
-12 472.0 -258.0 
-13 377.0 -411.0 
-14 428.0 -257.0 
-15 386.0 -278.0 
-16 372.0 -325.0 
-17 369.0 -371.0 
-18 433.0 -422.0 
-19 476.0 -361.0 
-20 502.0 -304.0 
-21 333.0 -225.0 
-22 305.0 -244.0 
-23 333.0 -270.0 
-24 369.0 -237.0 
+0 554.0 -242.0
+1 443.0 -190.0
+2 340.0 -184.0
+3 276.0 -237.0
+4 281.0 -347.0
+5 282.0 -380.0
+6 312.0 -425.0
+7 387.0 -472.0
+8 465.0 -487.0
+9 549.0 -469.0
+10 566.0 -413.0
+11 589.0 -359.0
+12 472.0 -258.0
+13 377.0 -411.0
+14 428.0 -257.0
+15 386.0 -278.0
+16 372.0 -325.0
+17 369.0 -371.0
+18 433.0 -422.0
+19 476.0 -361.0
+20 502.0 -304.0
+21 333.0 -225.0
+22 305.0 -244.0
+23 333.0 -270.0
+24 369.0 -237.0
 25 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments...
-0 14 12 
-1 15 14 
-2 16 15 
-3 17 16 
-4 13 17 
-5 18 13 
-6 19 18 
-7 20 19 
-8 12 20 
-9 22 21 
-10 23 22 
-11 24 23 
-12 21 24 
-13 3 4 
-14 6 7 
-15 5 6 
-16 2 3 
-17 4 5 
-18 7 8 
-19 9 10 
-20 8 9 
-21 0 1 
-22 1 2 
-23 11 0 
-24 10 11 
+0 14 12
+1 15 14
+2 16 15
+3 17 16
+4 13 17
+5 18 13
+6 19 18
+7 20 19
+8 12 20
+9 22 21
+10 23 22
+11 24 23
+12 21 24
+13 3 4
+14 6 7
+15 5 6
+16 2 3
+17 4 5
+18 7 8
+19 9 10
+20 8 9
+21 0 1
+22 1 2
+23 11 0
+24 10 11
 1 # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...
 0 418.0 -332.0
 1 # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions...
-0 334.0 -243.0 
+0 334.0 -243.0
 1 # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions...
-0 
+0
 # Geo reference
 56

inundation/ga/storm_surge/pyvolution/domain.py

@@ -260 +260 @@
         x = self.boundary.keys()
         x.sort()
+
         for k, (vol_id, edge_id) in enumerate(x):
             tag = self.boundary[ (vol_id, edge_id) ]

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -6 +6 @@
 
     A triangular element is defined in terms of three vertex ids,
-    ordered counter clock-wise, 
+    ordered counter clock-wise,
     each corresponding to a given coordinate set.
     Vertices from different elements can point to the same
@@ -12 +12 @@
 
     Coordinate sets are implemented as an N x 2 Numeric array containing
-    x and y coordinates. 
-    
+    x and y coordinates.
+
 
     To instantiate:
@@ -33 +33 @@
         c = [2.0,0.0]
         e = [2.0, 2.0]
-        
+
         points = [a, b, c, e]
-        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
         mesh = Mesh(points, triangles)
-    
+
         #creates two triangles: bac and bce
 
     Other:
-    
+
       In addition mesh computes an Nx6 array called vertex_coordinates.
-      This structure is derived from coordinates and contains for each 
+      This structure is derived from coordinates and contains for each
       triangle the three x,y coordinates at the vertices.
-                                
-
-        This is a cut down version of mesh from pyvolution\mesh.py
+
+
+        This is a cut down version of mesh from pyvolution mesh.py
     """
 
@@ -58 +58 @@
 
         origin is a 3-tuple consisting of UTM zone, easting and northing.
-        If specified coordinates are assumed to be relative to this origin. 
+        If specified coordinates are assumed to be relative to this origin.
         """
 
@@ -67 +67 @@
 
         if geo_reference is None:
-            self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) 
+            self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0)
         else:
             self.geo_reference = geo_reference
 
         #Input checks
-        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
+        msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples'
         assert len(self.triangles.shape) == 2, msg
 
         msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
-        assert len(self.coordinates.shape) == 2, msg        
+        assert len(self.coordinates.shape) == 2, msg
 
         msg = 'Vertex indices reference non-existing coordinate sets'
@@ -84 +84 @@
         #Register number of elements (N)
         self.number_of_elements = N = self.triangles.shape[0]
-    
-        #Allocate space for geometric quantities        
+
+        #Allocate space for geometric quantities
         self.normals = zeros((N, 6), Float)
         self.areas = zeros(N, Float)
@@ -92 +92 @@
         #Get x,y coordinates for all triangles and store
         self.vertex_coordinates = V = self.compute_vertex_coordinates()
-        
+
         #Initialise each triangle
         for i in range(N):
             #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
-            
+
             x0 = V[i, 0]; y0 = V[i, 1]
             x1 = V[i, 2]; y1 = V[i, 3]
-            x2 = V[i, 4]; y2 = V[i, 5]            
+            x2 = V[i, 4]; y2 = V[i, 5]
 
             #Area
@@ -107 +107 @@
             msg += ' is degenerate:  area == %f' %self.areas[i]
             assert self.areas[i] > 0.0, msg
-        
+
 
             #Normals
@@ -120 +120 @@
 
             n0 = array([x2 - x1, y2 - y1])
-            l0 = sqrt(sum(n0**2))
+            l0 = sqrt(sum(n0**2))
 
             n1 = array([x0 - x2, y0 - y2])
@@ -137 +137 @@
                                   n1[1], -n1[0],
                                   n2[1], -n2[0]]
-            
+
             #Edgelengths
             self.edgelengths[i, :] = [l0, l1, l2]
 
         #Build vertex list
-        self.build_vertexlist()        
-            
+        self.build_vertexlist()
+
     def __len__(self):
         return self.number_of_elements
@@ -154 +154 @@
         """Return all normal vectors.
         Return normal vectors for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
         """
         return self.normals
@@ -163 +163 @@
         Return value is the numeric array slice [x, y]
         """
-        return self.normals[i, 2*j:2*j+2]     
-    
-
-            
+        return self.normals[i, 2*j:2*j+2]
+
+
+
     def get_vertex_coordinates(self):
         """Return all vertex coordinates.
         Return all vertex coordinates for all triangles as an Nx6 array
-        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
         """
         return self.vertex_coordinates
@@ -179 +179 @@
         Return value is the numeric array slice [x, y]
         """
-        return self.vertex_coordinates[i, 2*j:2*j+2]     
+        return self.vertex_coordinates[i, 2*j:2*j+2]
 
 
@@ -189 +189 @@
         #FIXME: Perhaps they should be ordered as in obj files??
         #See quantity.get_vertex_values
-        
+
         from Numeric import zeros, Float
 
         N = self.number_of_elements
-        vertex_coordinates = zeros((N, 6), Float) 
+        vertex_coordinates = zeros((N, 6), Float)
 
         for i in range(N):
@@ -203 +203 @@
 
         return vertex_coordinates
-    
+
     def get_vertices(self,  indexes=None):
         """Get connectivity
         indexes is the set of element ids of interest
         """
-        
+
         from Numeric import take
-        
+
         if (indexes ==  None):
             indexes = range(len(self))  #len(self)=number of elements
-            
+
         return  take(self.triangles, indexes)
 
@@ -220 +220 @@
         unique_verts = {}
         for triangle in triangles:
-           unique_verts[triangle[0]] = 0 
-           unique_verts[triangle[1]] = 0 
+           unique_verts[triangle[0]] = 0
+           unique_verts[triangle[1]] = 0
            unique_verts[triangle[2]] = 0
-        return unique_verts.keys()   
-        
+        return unique_verts.keys()
+
     def build_vertexlist(self):
         """Build vertexlist index by vertex ids and for each entry (point id)
@@ -231 +231 @@
 
         Preconditions:
-          self.coordinates and self.triangles are defined  
-        
+          self.coordinates and self.triangles are defined
+
         Postcondition:
           self.vertexlist is built
@@ -242 +242 @@
             a = self.triangles[i, 0]
             b = self.triangles[i, 1]
-            c = self.triangles[i, 2]            
+            c = self.triangles[i, 2]
 
             #Register the vertices v as lists of
@@ -251 +251 @@
                     vertexlist[v] = []
 
-                vertexlist[v].append( (i, vertex_id) )    
+                vertexlist[v].append( (i, vertex_id) )
 
         self.vertexlist = vertexlist
-    
+
 
     def get_extent(self):
@@ -264 +264 @@
         X = C[:,0:6:2].copy()
         Y = C[:,1:6:2].copy()
-        
+
         xmin = min(X.flat)
         xmax = max(X.flat)
         ymin = min(Y.flat)
-        ymax = max(Y.flat)                
-        
+        ymax = max(Y.flat)
+
         return xmin, xmax, ymin, ymax
 
@@ -277 +277 @@
         """
 
-        return sum(self.areas)          
-        
+        return sum(self.areas)

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -130 +130 @@
                 c = sqrt((x2-x0)**2+(y2-y0)**2)
 
-                self.radii[i]=2.0*self.areas[i]/(a+b+c)
+                self.radii[i]=self.areas[i]/(a+b+c)
 
 

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -4 +4 @@
 
 
-def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+def rectangular_old(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
 
     """Setup a rectangular grid of triangles
@@ -10 +10 @@
     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)
 
     Return to lists: points and elements suitable for creating a Mesh or
-    FVMesh object, e.g. Mesh(points, elements) 
+    FVMesh object, e.g. Mesh(points, elements)
     """
 
     from config import epsilon
-    
+
     #E = m*n*2        #Number of triangular elements
     #P = (m+1)*(n+1)  #Number of initial vertices
 
     delta1 = float(len1)/m
-    delta2 = float(len2)/n    
+    delta2 = float(len2)/n
 
     #Dictionary of vertex objects
@@ -43 +43 @@
         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]            
+            v2 = vertices[i,j]
+            v3 = vertices[i+1,j+1]
+            v4 = vertices[i+1,j]
 
             #Update boundary dictionary and create elements
@@ -51 +51 @@
                 boundary[(len(elements), 2)] = 'right'
             if j == 0:
-                boundary[(len(elements), 1)] = 'bottom'                
+                boundary[(len(elements), 1)] = 'bottom'
             elements.append([v4,v3,v2]) #Lower element
 
@@ -57 +57 @@
                 boundary[(len(elements), 2)] = 'left'
             if j == n-1:
-                boundary[(len(elements), 1)] = 'top' 
-            elements.append([v1,v2,v3]) #Upper element   
-
-    return points, elements, boundary       
-
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v2,v3]) #Upper element
+
+    return points, elements, boundary
+
+def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+    """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)
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements)
+    """
+
+    from config import epsilon
+    from Numeric import zeros, Float, Int
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n
+
+    #Calculate number of points
+    Np = (m+1)*(n+1)
+
+    class Index:
+
+        def __init__(self, n,m):
+            self.n = n
+            self.m = m
+
+        def __call__(self, i,j):
+            return j+i*(self.n+1)
+
+
+    index = Index(n,m)
+
+    points = zeros( (Np,2), Float)
+
+    for i in range(m+1):
+        for j in range(n+1):
+
+            points[index(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
+
+    #Construct 2 triangles per rectangular element and assign tags to boundary
+    #Calculate number of triangles
+    Nt = 2*m*n
+
+
+    elements = zeros( (Nt,3), Int)
+    boundary = {}
+    nt = -1
+    for i in range(m):
+        for j in range(n):
+            nt = nt + 1
+            i1 = index(i,j+1)
+            i2 = index(i,j)
+            i3 = index(i+1,j+1)
+            i4 = index(i+1,j)
+
+            #Update boundary dictionary and create elements
+            if i == m-1:
+                boundary[nt, 2] = 'right'
+            if j == 0:
+                boundary[nt, 1] = 'bottom'
+            elements[nt,:] = [i4,i3,i2] #Lower element
+            nt = nt + 1
+
+            if i == 0:
+                boundary[nt, 2] = 'left'
+            if j == n-1:
+                boundary[nt, 1] = 'top'
+            elements[nt,:] = [i1,i2,i3] #Upper element
+
+    return points, elements, boundary
 
 def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
@@ -97 +170 @@
     for i in range(m+1):
         x = deltax*i
-        for j in range(n+1):     
+        for j in range(n+1):
             y = deltay*j
             if x > 0.25*lenx:
                 y = a1 + a2*x + a3*y + a4*x*y
-                
+
             vertices[i,j] = len(points)
             points.append([x + origin[0], y + origin[1]])
 
-    # Construct 2 triangles per element        
+    # Construct 2 triangles per element
     elements = []
     boundary = {}
@@ -111 +184 @@
         for j in range(n):
             v1 = vertices[i,j+1]
-            v2 = vertices[i,j]            
-            v3 = vertices[i+1,j+1]            
+            v2 = vertices[i,j]
+            v3 = vertices[i+1,j+1]
             v4 = vertices[i+1,j]
-            
+
             #Update boundary dictionary and create elements
             if i == m-1:
                 boundary[(len(elements), 2)] = 'right'
             if j == 0:
-                boundary[(len(elements), 1)] = 'bottom'                
+                boundary[(len(elements), 1)] = 'bottom'
             elements.append([v4,v3,v2]) #Lower
 
@@ -125 +198 @@
                 boundary[(len(elements), 2)] = 'left'
             if j == n-1:
-                boundary[(len(elements), 1)] = 'top' 
-             
+                boundary[(len(elements), 1)] = 'top'
+
             elements.append([v1,v2,v3]) #Upper
 
-    return points, elements, boundary       
+    return points, elements, boundary
 
 
@@ -136 +209 @@
     with n radial segments. If radius is are omitted the mesh defaults to
     the unit circle radius.
- 
+
     radius: radius of circle
-    
+
     #FIXME: The triangles become degenerate for large values of m or n.
     """
 
 
-    
+
     from math import pi, cos, sin
 
@@ -152 +225 @@
     points = [[0.0, 0.0]] #Center point
     vertices[0, 0] = 0
-    
+
     for i in range(n):
         theta = 2*i*pi/n
@@ -162 +235 @@
             points.append([delta*x, delta*y])
 
-    #Construct 2 triangles per element        
+    #Construct 2 triangles per element
     elements = []
     for i in range(n):
         for j in range(1,m):
 
-            i1 = (i + 1) % n  #Wrap around  
-            
+            i1 = (i + 1) % n  #Wrap around
+
             v1 = vertices[i,j+1]
-            v2 = vertices[i,j]            
-            v3 = vertices[i1,j+1]            
+            v2 = vertices[i,j]
+            v3 = vertices[i1,j+1]
             v4 = vertices[i1,j]
 
-            elements.append([v4,v2,v3]) #Lower                
+            elements.append([v4,v2,v3]) #Lower
             elements.append([v1,v3,v2]) #Upper
 
@@ -181 +254 @@
     v1 = vertices[0,0]
     for i in range(n):
-        i1 = (i + 1) % n  #Wrap around          
+        i1 = (i + 1) % n  #Wrap around
         v2 = vertices[i,1]
         v3 = vertices[i1,1]
 
         elements.append([v1,v2,v3]) #center
-         
+
     return points, elements
 
@@ -205 +278 @@
         filename = name + '.poly'
 
-    
+
     fid = open(filename)
 
     points = []    #x, y
     values = []    #z
-    ##vertex_values = []    #Repeated z    
+    ##vertex_values = []    #Repeated z
     triangles = [] #v0, v1, v2
-    
+
     lines = fid.readlines()
 
@@ -237 +310 @@
         x = float(xyz[0])
         y = float(xyz[1])
-        z = float(xyz[2])        
+        z = float(xyz[2])
 
         points.append([x, y])
         values.append(z)
-        
-
-    k = i    
+
+
+    k = i
     while keyword == 'POLYS':
         line = lines[i].strip()
@@ -251 +324 @@
             keyword = line
             break
-        
+
 
         fields = line.split(':')
@@ -269 +342 @@
         y1 = points[i1][1]
         x2 = points[i2][0]
-        y2 = points[i2][1]                
+        y2 = points[i2][1]
 
         area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
@@ -281 +354 @@
             a0 = anglediff(v1, v0)
             a1 = anglediff(v2, v1)
-            a2 = anglediff(v0, v2)        
-
-            
+            a2 = anglediff(v0, v2)
+
+
             if a0 < pi and a1 < pi and a2 < pi:
                 #all is well
@@ -294 +367 @@
                 j1 = i0
                 j2 = i2
-                
+
             triangles.append([j0, j1, j2])
             ##vertex_values.append([values[j0], values[j1], values[j2]])
         else:
-            offending +=1 
+            offending +=1
 
     print 'Removed %d offending triangles out of %d' %(offending, len(lines))
@@ -311 +384 @@
     from math import pi
     from util import anglediff
-    
-    
+
+
     fid = open(filename)
     points = []    # List of x, y coordinates
     triangles = [] # List of vertex ids as listed in the file
-    
+
     for line in fid.readlines():
         fields = line.split()
@@ -337 +410 @@
         v1 = t[1]
         v2 = t[2]
-        
+
         x0 = points[v0][0]
         y0 = points[v0][1]
@@ -343 +416 @@
         y1 = points[v1][1]
         x2 = points[v2][0]
-        y2 = points[v2][1]                
+        y2 = points[v2][1]
 
         #Check that points are arranged in counter clock-wise order
@@ -352 +425 @@
         a0 = anglediff(vec1, vec0)
         a1 = anglediff(vec2, vec1)
-        a2 = anglediff(vec0, vec2)        
+        a2 = anglediff(vec0, vec2)
 
         if a0 < pi and a1 < pi and a2 < pi:
             elements.append([v0, v1, v2])
         else:
-            elements.append([v0, v2, v1])       
-
-    return points, elements                             
+            elements.append([v0, v2, v1])
+
+    return points, elements
 
 # #Map from edge number to indices of associated vertices
@@ -365 +438 @@
 
 
-                                         
-    

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -10 +10 @@
 #
 from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
-     Transmissive_boundary, Constant_height
+     Transmissive_boundary, Constant_height, Constant_stage
 
 from mesh_factory import rectangular
@@ -80 +80 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 60
+N = 50
+M = 40
 
 #N = 15
@@ -86 +87 @@
 print 'Creating domain'
 #Create basic mesh
-points, vertices, boundary = rectangular(N, N)
+points, vertices, boundary = rectangular(N, M)
 
 #Create shallow water domain
@@ -102 +103 @@
 
 #Set bed-slope and friction
-inflow_stage = 0.1
-manning = 0.02
+inflow_stage = 0.2
+manning = 0.00
 Z = Weir(inflow_stage)
 
@@ -148 +149 @@
 #
 print 'Initial condition'
-domain.set_quantity('stage', Constant_height(Z, 0.))
+domain.set_quantity('stage', Constant_height(Z, 0.0))
+#domain.set_quantity('stage', Constant_stage(-1.0))
 
 #Evolve
@@ -154 +156 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0):
+for t in domain.evolve(yieldstep = 0.01, finaltime = 5.0):
     domain.write_time()
+    #domain.visualiser.scale_z = 1.0
+    domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
 
 

inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

@@ -3 +3 @@
 elevation_color = (0.3,0.4,0.5)
 stage_color = (0.1,0.4,0.99)
+other_color = (0.99,0.4,0.1)
 
 class Visualiser:
@@ -117 +118 @@
 
 
-    def update_quantity(self,qname,qcolor=None):
+    def update_quantity(self,qname,qcolor=None,scale_z=None):
 
         #print 'update '+qname+' arrays'
 
         if qcolor is None:
+            qcolor = other_color
+
             if qname=='elevation':
                 qcolor = elevation_color
@@ -128 +131 @@
                 qcolor = stage_color
 
+        if scale_z is None:
+            scale_z = self.scale_z
 
         try:
-            self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor)
+            self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor, scale_z)
 
             #print 'update bed image'
@@ -140 +145 @@
             pass
 
-    def update_quantity_color(self,qname,qcolor=None):
+    def update_quantity_color(self,qname,qcolor=None,scale_z=None):
 
         #print 'update '+qname+' arrays'
@@ -151 +156 @@
                 qcolor = stage_color
 
+        if scale_z is None:
+            scale_z = self.scale_z
 
         #try:
-        self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor)
+        self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor, scale_z)
 
             #print 'update bed image'
@@ -164 +171 @@
 
 
-    def update_arrays(self,quantity,qcolor):
+    def update_arrays(self,quantity,qcolor,scale_z):
 
         col = asarray(qcolor)
@@ -170 +177 @@
         N = len(self.domain)
 
-        scale_z = self.scale_z
         min_x = self.min_x
         min_y = self.min_y
@@ -189 +195 @@
 
 
-    def update_arrays_color(self,quantity,qcolor):
+    def update_arrays_color(self,quantity,qcolor,scale_z):
 
         col = asarray(qcolor)
@@ -195 +201 @@
         N = len(self.domain)
 
-        scale_z = self.scale_z
         min_x = self.min_x
         min_y = self.min_y
@@ -308 +313 @@
             normal(1) = normal(1)/norm;
             normal(2) = normal(2)/norm;
+
 
             for (int j=0; j<3; j++) {

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -1507 +1507 @@
 
 
+class Constant_stage:
+    """Set an initial condition with constant stage
+    """
+    def __init__(self, s):
+        self.s = s
+
+    def __call__(self, x, y):
+        return self.s
+
 class Constant_height:
     """Set an initial condition with constant water height, e.g
     stage s = z+h
     """
+
     def __init__(self, W, h):
         self.W = W

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -2 +2 @@
 //
 // To compile (Python2.3):
-//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O 
-//  gcc -shared domain_ext.o  -o domain_ext.so  
+//  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
+//  gcc -shared domain_ext.o  -o domain_ext.so
 //
 // or use python compile.py
@@ -11 +11 @@
 //
 // Ole Nielsen, GA 2004
-    
-    
+
+
 #include "Python.h"
 #include "Numeric/arrayobject.h"
 #include "math.h"
+#include <stdio.h>
 
 //Shared code snippets
@@ -26 +27 @@
   /*Rotate the momentum component q (q[1], q[2])
     from x,y coordinates to coordinates based on normal vector (n1, n2).
-    
-    Result is returned in array 3x1 r     
+
+    Result is returned in array 3x1 r
     To rotate in opposite direction, call rotate with (q, n1, -n2)
-    
-    Contents of q are changed by this function */    
+
+    Contents of q are changed by this function */
 
 
   double q1, q2;
-  
+
   //Shorthands
   q1 = q[1];  //uh momentum
@@ -41 +42 @@
   //Rotate
   q[1] =  n1*q1 + n2*q2;
-  q[2] = -n2*q1 + n1*q2; 
+  q[2] = -n2*q1 + n1*q2;
 
   return 0;
-}  
-
-
-        
+}
+
+
+
 // Computational function for flux computation (using stage w=z+h)
-int flux_function(double *q_left, double *q_right, 
-                  double z_left, double z_right, 
-                  double n1, double n2, 
-                  double epsilon, double g, 
+int flux_function(double *q_left, double *q_right,
+                  double z_left, double z_right,
+                  double n1, double n2,
+                  double epsilon, double g,
                   double *edgeflux, double *max_speed) {
-  
+
   /*Compute fluxes between volumes for the shallow water wave equation
-    cast in terms of the 'stage', w = h+z using 
+    cast in terms of the 'stage', w = h+z using
     the 'central scheme' as described in
-    
+
     Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
     Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
     Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
-    
-    The implemented formula is given in equation (3.15) on page 714 
+
+    The implemented formula is given in equation (3.15) on page 714
   */
-       
+
   int i;
-       
+
   double w_left, h_left, uh_left, vh_left, u_left;
   double w_right, h_right, uh_right, vh_right, u_right;
   double s_min, s_max, soundspeed_left, soundspeed_right;
   double denom, z;
-  double q_left_copy[3], q_right_copy[3];    
+  double q_left_copy[3], q_right_copy[3];
   double flux_right[3], flux_left[3];
-  
+
   //Copy conserved quantities to protect from modification
   for (i=0; i<3; i++) {
     q_left_copy[i] = q_left[i];
     q_right_copy[i] = q_right[i];
-  } 
-    
+  }
+
   //Align x- and y-momentum with x-axis
   _rotate(q_left_copy, n1, n2);
-  _rotate(q_right_copy, n1, n2);    
+  _rotate(q_right_copy, n1, n2);
 
   z = (z_left+z_right)/2; //Take average of field values
@@ -95 +96 @@
     h_left = 0.0;  //Could have been negative
     u_left = 0.0;
-  } else { 
+  } else {
     u_left = uh_left/h_left;
   }
-  
+
   w_right = q_right_copy[0];
   h_right = w_right-z;
@@ -106 +107 @@
     h_right = 0.0; //Could have been negative
     u_right = 0.0;
-  } else { 
+  } else {
     u_right = uh_right/h_right;
   }
 
-  //Momentum in y-direction              
+  //Momentum in y-direction
   vh_left  = q_left_copy[2];
-  vh_right = q_right_copy[2];   
-       
+  vh_right = q_right_copy[2];
+
 
   //Maximal and minimal wave speeds
-  soundspeed_left  = sqrt(g*h_left); 
+  soundspeed_left  = sqrt(g*h_left);
   soundspeed_right = sqrt(g*h_right);
-    
+
   s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
-  if (s_max < 0.0) s_max = 0.0; 
-  
+  if (s_max < 0.0) s_max = 0.0;
+
   s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
-  if (s_min > 0.0) s_min = 0.0;   
-  
-  //Flux formulas  
+  if (s_min > 0.0) s_min = 0.0;
+
+  //Flux formulas
   flux_left[0] = u_left*h_left;
   flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
   flux_left[2] = u_left*vh_left;
-  
+
   flux_right[0] = u_right*h_right;
   flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
   flux_right[2] = u_right*vh_right;
-    
-
-  //Flux computation    
+
+
+  //Flux computation
   denom = s_max-s_min;
   if (denom == 0.0) {
     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     *max_speed = 0.0;
-  } else {    
+  } else {
     for (i=0; i<3; i++) {
       edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
       edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]);
       edgeflux[i] /= denom;
-    }  
-        
+    }
+
     //Maximal wavespeed
     *max_speed = max(fabs(s_max), fabs(s_min));
-    
-    //Rotate back       
+
+    //Rotate back
     _rotate(edgeflux, n1, -n2);
   }
   return 0;
 }
-        
-void _manning_friction(double g, double eps, int N, 
-                       double* w, double* z, 
-                       double* uh, double* vh, 
-                       double* eta, double* xmom, double* ymom) {     
+
+void _manning_friction(double g, double eps, int N,
+                       double* w, double* z,
+                       double* uh, double* vh,
+                       double* eta, double* xmom, double* ymom) {
 
   int k;
   double S, h;
-  
+
   for (k=0; k<N; k++) {
     if (eta[k] > eps) {
@@ -180 +181 @@
   }
 }
-            
+
 
 
 int _balance_deep_and_shallow(int N,
                               double* wc,
-                              double* zc, 
-                              double* hc,                             
-                              double* wv, 
-                              double* zv, 
+                              double* zc,
+                              double* hc,
+                              double* wv,
+                              double* zv,
                               double* hv,
-                              double* hvbar,                          
-                              double* xmomc, 
-                              double* ymomc, 
-                              double* xmomv, 
-                              double* ymomv) { 
-  
+                              double* hvbar,
+                              double* xmomc,
+                              double* ymomc,
+                              double* xmomv,
+                              double* ymomv) {
+
   int k, k3, i;
   double dz, hmin, alpha;
-  
+
   //Compute linear combination between w-limited stages and
-  //h-limited stages close to the bed elevation.     
-  
+  //h-limited stages close to the bed elevation.
+
   for (k=0; k<N; k++) {
     // Compute maximal variation in bed elevation
@@ -211 +212 @@
 
     k3 = 3*k;
-    
+
     //FIXME: Try with this one precomputed
     dz = 0.0;
@@ -220 +221 @@
     }
 
-    
-    //Create alpha in [0,1], where alpha==0 means using the h-limited 
+
+    //Create alpha in [0,1], where alpha==0 means using the h-limited
     //stage and alpha==1 means using the w-limited stage as
     //computed by the gradient limiter (both 1st or 2nd order)
@@ -227 +228 @@
     //If hmin > dz/2 then alpha = 1 and the bed will have no effect
     //If hmin < 0 then alpha = 0 reverting to constant height above bed.
-   
+
 
     if (dz > 0.0)
-      //if (hmin<0.0) 
+      //if (hmin<0.0)
       //        alpha = 0.0;
       //else
       //  alpha = max( min( hc[k]/dz, 1.0), 0.0 );
       alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 );
-   
     else
       alpha = 1.0;  //Flat bed
 
+    //alpha = 1.0;
+
     //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
 
-    //  Let 
-    //  
-    //    wvi be the w-limited stage (wvi = zvi + hvi)        
+    //  Let
+    //
+    //    wvi be the w-limited stage (wvi = zvi + hvi)
     //    wvi- be the h-limited state (wvi- = zvi + hvi-)
-    //      
-    //      
+    //
+    //
     //  where i=0,1,2 denotes the vertex ids
-    //  
-    //  Weighted balance between w-limited and h-limited stage is  
-    //    
-    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
-    //  
+    //
+    //  Weighted balance between w-limited and h-limited stage is
+    //
+    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
+    //
     //  It follows that the updated wvi is
     //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
-    //   
+    //
     //   Momentum is balanced between constant and limited
 
-    if (alpha < 1) {         
+    if (alpha < 1) {
       for (i=0; i<3; i++) {
         wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
-            
-        //Update momentum as a linear combination of 
+
+        //Update momentum as a linear combination of
         //xmomc and ymomc (shallow) and momentum
         //from extrapolator xmomv and ymomv (deep).
-        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];            
+        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
         ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
-      } 
+      }
     }
-  }         
+  }
   return 0;
 }
@@ -275 +277 @@
 
 
-int _protect(int N, 
-             double minimum_allowed_height,        
+int _protect(int N,
+             double minimum_allowed_height,
              double* wc,
-             double* zc, 
-             double* xmomc, 
+             double* zc,
+             double* xmomc,
              double* ymomc) {
-  
-  int k; 
+
+  int k;
   double hc;
-  
+
   //Protect against initesimal and negative heights
-  
+
   for (k=0; k<N; k++) {
     hc = wc[k] - zc[k];
-    
-    if (hc < minimum_allowed_height) {    
-      wc[k] = zc[k];       
+
+    if (hc < minimum_allowed_height) {
+      wc[k] = zc[k];
       xmomc[k] = 0.0;
-      ymomc[k] = 0.0;      
+      ymomc[k] = 0.0;
     }
-    
+
   }
   return 0;
@@ -309 +311 @@
 int _assign_wind_field_values(int N,
                               double* xmom_update,
-                              double* ymom_update, 
+                              double* ymom_update,
                               double* s_vec,
                               double* phi_vec,
@@ -318 +320 @@
   int k;
   double S, s, phi, u, v;
-  
+
   for (k=0; k<N; k++) {
-    
+
     s = s_vec[k];
     phi = phi_vec[k];
@@ -330 +332 @@
     u = s*cos(phi);
     v = s*sin(phi);
-        
+
     //Compute wind stress
     S = cw * sqrt(u*u + v*v);
     xmom_update[k] += S*u;
-    ymom_update[k] += S*v;        
-  }
-  return 0;  
-}            
+    ymom_update[k] += S*v;
+  }
+  return 0;
+}
 
 
@@ -348 +350 @@
   //  gravity(g, h, v, x, xmom, ymom)
   //
-  
-  
+
+
   PyArrayObject *h, *v, *x, *xmom, *ymom;
   int k, i, N, k3, k6;
   double g, avg_h, zx, zy;
   double x0, y0, x1, y1, x2, y2, z0, z1, z2;
-    
+
   if (!PyArg_ParseTuple(args, "dOOOOO",
-                        &g, &h, &v, &x, 
-                        &xmom, &ymom))  
-    return NULL;  
+                        &g, &h, &v, &x,
+                        &xmom, &ymom))
+    return NULL;
 
   N = h -> dimensions[0];
   for (k=0; k<N; k++) {
     k3 = 3*k;  // base index
-    k6 = 6*k;  // base index    
-    
+    k6 = 6*k;  // base index
+
     avg_h = 0.0;
     for (i=0; i<3; i++) {
       avg_h += ((double *) h -> data)[k3+i];
-    }   
+    }
     avg_h /= 3;
-        
-    
+
+
     //Compute bed slope
     x0 = ((double*) x -> data)[k6 + 0];
-    y0 = ((double*) x -> data)[k6 + 1];    
+    y0 = ((double*) x -> data)[k6 + 1];
     x1 = ((double*) x -> data)[k6 + 2];
-    y1 = ((double*) x -> data)[k6 + 3];        
+    y1 = ((double*) x -> data)[k6 + 3];
     x2 = ((double*) x -> data)[k6 + 4];
-    y2 = ((double*) x -> data)[k6 + 5];           
+    y2 = ((double*) x -> data)[k6 + 5];
 
 
     z0 = ((double*) v -> data)[k3 + 0];
     z1 = ((double*) v -> data)[k3 + 1];
-    z2 = ((double*) v -> data)[k3 + 2];        
+    z2 = ((double*) v -> data)[k3 + 2];
 
     _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
@@ -389 +391 @@
     //Update momentum
     ((double*) xmom -> data)[k] += -g*zx*avg_h;
-    ((double*) ymom -> data)[k] += -g*zy*avg_h;        
-  }
-    
+    ((double*) ymom -> data)[k] += -g*zy*avg_h;
+  }
+
   return Py_BuildValue("");
 }
@@ -400 +402 @@
   // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
   //
-  
-  
+
+
   PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
   int N;
   double g, eps;
-    
+
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
-                        &g, &eps, &w, &z, &uh, &vh, &eta, 
-                        &xmom, &ymom))  
-    return NULL;  
-
-  N = w -> dimensions[0];    
+                        &g, &eps, &w, &z, &uh, &vh, &eta,
+                        &xmom, &ymom))
+    return NULL;
+
+  N = w -> dimensions[0];
   _manning_friction(g, eps, N,
                     (double*) w -> data,
-                    (double*) z -> data,                    
-                    (double*) uh -> data, 
-                    (double*) vh -> data, 
+                    (double*) z -> data,
+                    (double*) uh -> data,
+                    (double*) vh -> data,
                     (double*) eta -> data,
-                    (double*) xmom -> data, 
+                    (double*) xmom -> data,
                     (double*) ymom -> data);
 
   return Py_BuildValue("");
-}                   
+}
 
 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
@@ -431 +433 @@
   // normal a Float numeric array of length 2.
 
-  
+
   PyObject *Q, *Normal;
   PyArrayObject *q, *r, *normal;
-  
+
   static char *argnames[] = {"q", "normal", "direction", NULL};
   int dimensions[1], i, direction=1;
   double n1, n2;
 
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
-                                   &Q, &Normal, &direction))  
-    return NULL;  
+  // Convert Python arguments to C
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
+                                   &Q, &Normal, &direction))
+    return NULL;
 
   //Input checks (convert sequences into numeric arrays)
-  q = (PyArrayObject *) 
+  q = (PyArrayObject *)
     PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
-  normal = (PyArrayObject *) 
+  normal = (PyArrayObject *)
     PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
 
-  
+
   if (normal -> dimensions[0] != 2) {
     PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
@@ -456 +458 @@
   }
 
-  //Allocate space for return vector r (don't DECREF) 
+  //Allocate space for return vector r (don't DECREF)
   dimensions[0] = 3;
   r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
@@ -462 +464 @@
   //Copy
   for (i=0; i<3; i++) {
-    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];  
-  }
-  
+    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
+  }
+
   //Get normal and direction
-  n1 = ((double *) normal -> data)[0]; 
-  n2 = ((double *) normal -> data)[1];   
+  n1 = ((double *) normal -> data)[0];
+  n2 = ((double *) normal -> data)[1];
   if (direction == -1) n2 = -n2;
 
@@ -474 +476 @@
 
   //Release numeric arrays
-  Py_DECREF(q);    
+  Py_DECREF(q);
   Py_DECREF(normal);
-        
+
   //return result using PyArray to avoid memory leak
   return PyArray_Return(r);
-}    
+}
 
 
@@ -490 +492 @@
     Fluxes across each edge are scaled by edgelengths and summed up
     Resulting flux is then scaled by area and stored in
-    explicit_update for each of the three conserved quantities 
+    explicit_update for each of the three conserved quantities
     stage, xmomentum and ymomentum
 
@@ -496 +498 @@
     is converted to a timestep that must not be exceeded. The minimum of
     those is computed as the next overall timestep.
-    
+
     Python call:
-    domain.timestep = compute_fluxes(timestep, 
-                                     domain.epsilon, 
+    domain.timestep = compute_fluxes(timestep,
+                                     domain.epsilon,
                                      domain.g,
                                      domain.neighbours,
                                      domain.neighbour_edges,
                                      domain.normals,
-                                     domain.edgelengths,                       
+                                     domain.edgelengths,
                                      domain.radii,
                                      domain.areas,
-                                     Stage.edge_values, 
-                                     Xmom.edge_values, 
-                                     Ymom.edge_values,  
-                                     Bed.edge_values,   
+                                     Stage.edge_values,
+                                     Xmom.edge_values,
+                                     Ymom.edge_values,
+                                     Bed.edge_values,
                                      Stage.boundary_values,
                                      Xmom.boundary_values,
@@ -517 +519 @@
                                      Xmom.explicit_update,
                                      Ymom.explicit_update)
-        
+
 
     Post conditions:
       domain.explicit_update is reset to computed flux values
-      domain.timestep is set to the largest step satisfying all volumes. 
-
-            
+      domain.timestep is set to the largest step satisfying all volumes.
+
+
   */
 
-  
+
   PyArrayObject *neighbours, *neighbour_edges,
     *normals, *edgelengths, *radii, *areas,
-    *stage_edge_values, 
-    *xmom_edge_values, 
-    *ymom_edge_values,  
-    *bed_edge_values,   
+    *stage_edge_values,
+    *xmom_edge_values,
+    *ymom_edge_values,
+    *bed_edge_values,
     *stage_boundary_values,
     *xmom_boundary_values,
@@ -540 +542 @@
     *ymom_explicit_update;
 
-    
-  //Local variables  
+
+  //Local variables
   double timestep, max_speed, epsilon, g;
   double normal[2], ql[3], qr[3], zl, zr;
@@ -548 +550 @@
   int number_of_elements, k, i, j, m, n;
   int ki, nm, ki2; //Index shorthands
-  
- 
-  // Convert Python arguments to C  
+
+
+  // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
                         &timestep,
                         &epsilon,
                         &g,
-                        &neighbours, 
+                        &neighbours,
                         &neighbour_edges,
-                        &normals, 
+                        &normals,
                         &edgelengths, &radii, &areas,
-                        &stage_edge_values, 
-                        &xmom_edge_values, 
-                        &ymom_edge_values,  
-                        &bed_edge_values,   
+                        &stage_edge_values,
+                        &xmom_edge_values,
+                        &ymom_edge_values,
+                        &bed_edge_values,
                         &stage_boundary_values,
                         &xmom_boundary_values,
@@ -574 +576 @@
 
   number_of_elements = stage_edge_values -> dimensions[0];
-  
-    
+
+
   for (k=0; k<number_of_elements; k++) {
 
     //Reset work array
     for (j=0; j<3; j++) flux[j] = 0.0;
-                          
+
     //Loop through neighbours and compute edge flux for each
     for (i=0; i<3; i++) {
@@ -586 +588 @@
       ql[0] = ((double *) stage_edge_values -> data)[ki];
       ql[1] = ((double *) xmom_edge_values -> data)[ki];
-      ql[2] = ((double *) ymom_edge_values -> data)[ki];            
-      zl =    ((double *) bed_edge_values -> data)[ki];                  
+      ql[2] = ((double *) ymom_edge_values -> data)[ki];
+      zl =    ((double *) bed_edge_values -> data)[ki];
 
       //Quantities at neighbour on nearest face
@@ -593 +595 @@
       if (n < 0) {
         m = -n-1; //Convert negative flag to index
-        qr[0] = ((double *) stage_boundary_values -> data)[m];  
-        qr[1] = ((double *) xmom_boundary_values -> data)[m];   
-        qr[2] = ((double *) ymom_boundary_values -> data)[m];   
+        qr[0] = ((double *) stage_boundary_values -> data)[m];
+        qr[1] = ((double *) xmom_boundary_values -> data)[m];
+        qr[2] = ((double *) ymom_boundary_values -> data)[m];
         zr = zl; //Extend bed elevation to boundary
-      } else {    
+      } else {
         m = ((long *) neighbour_edges -> data)[ki];
-        
-        nm = n*3+m;     
+
+        nm = n*3+m;
         qr[0] = ((double *) stage_edge_values -> data)[nm];
         qr[1] = ((double *) xmom_edge_values -> data)[nm];
-        qr[2] = ((double *) ymom_edge_values -> data)[nm];            
-        zr =    ((double *) bed_edge_values -> data)[nm];                  
+        qr[2] = ((double *) ymom_edge_values -> data)[nm];
+        zr =    ((double *) bed_edge_values -> data)[nm];
       }
 
       //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
-      //             ql[0], ql[1], ql[2]);       
-
-      
-      // Outward pointing normal vector   
+      //             ql[0], ql[1], ql[2]);
+
+
+      // Outward pointing normal vector
       // normal = domain.normals[k, 2*i:2*i+2]
       ki2 = 2*ki; //k*6 + i*2
       normal[0] = ((double *) normals -> data)[ki2];
-      normal[1] = ((double *) normals -> data)[ki2+1];      
+      normal[1] = ((double *) normals -> data)[ki2+1];
 
       //Edge flux computation
-      flux_function(ql, qr, zl, zr, 
+      flux_function(ql, qr, zl, zr,
                     normal[0], normal[1],
-                    epsilon, g, 
+                    epsilon, g,
                     edgeflux, &max_speed);
 
-                    
+
       //flux -= edgeflux * edgelengths[k,i]
-      for (j=0; j<3; j++) { 
+      for (j=0; j<3; j++) {
         flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
       }
-      
+
       //Update timestep
       //timestep = min(timestep, domain.radii[k]/max_speed)
@@ -634 +636 @@
       if (max_speed > epsilon) {
         timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-      }   
+      }
     } // end for i
-    
+
     //Normalise by area and store for when all conserved
     //quantities get updated
     // flux /= areas[k]
-    for (j=0; j<3; j++) { 
+    for (j=0; j<3; j++) {
       flux[j] /= ((double *) areas -> data)[k];
     }
@@ -646 +648 @@
     ((double *) stage_explicit_update -> data)[k] = flux[0];
     ((double *) xmom_explicit_update -> data)[k] = flux[1];
-    ((double *) ymom_explicit_update -> data)[k] = flux[2];        
+    ((double *) ymom_explicit_update -> data)[k] = flux[2];
 
   } //end for k
 
   return Py_BuildValue("d", timestep);
-}    
+}
 
 
@@ -658 +660 @@
   //
   //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
-  
-
-  PyArrayObject 
+
+
+  PyArrayObject
   *wc,            //Stage at centroids
-  *zc,            //Elevation at centroids    
+  *zc,            //Elevation at centroids
   *xmomc,         //Momentums at centroids
-  *ymomc; 
-
-    
-  int N; 
+  *ymomc;
+
+
+  int N;
   double minimum_allowed_height;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "dOOOO",  
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOO",
                         &minimum_allowed_height,
                         &wc, &zc, &xmomc, &ymomc))
@@ -677 +679 @@
 
   N = wc -> dimensions[0];
-    
+
   _protect(N,
            minimum_allowed_height,
            (double*) wc -> data,
-           (double*) zc -> data, 
-           (double*) xmomc -> data, 
+           (double*) zc -> data,
+           (double*) xmomc -> data,
            (double*) ymomc -> data);
- 
-  return Py_BuildValue("");  
+
+  return Py_BuildValue("");
 }
 
@@ -694 +696 @@
   //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
   //                             xmomc, ymomc, xmomv, ymomv)
-  
-
-  PyArrayObject 
+
+
+  PyArrayObject
     *wc,            //Stage at centroids
-    *zc,            //Elevation at centroids    
-    *hc,            //Height at centroids        
+    *zc,            //Elevation at centroids
+    *hc,            //Height at centroids
     *wv,            //Stage at vertices
     *zv,            //Elevation at vertices
-    *hv,            //Depths at vertices    
-    *hvbar,         //h-Limited depths at vertices        
+    *hv,            //Depths at vertices
+    *hvbar,         //h-Limited depths at vertices
     *xmomc,         //Momentums at centroids and vertices
-    *ymomc, 
-    *xmomv, 
-    *ymomv;    
-    
+    *ymomc,
+    *xmomv,
+    *ymomv;
+
   int N; //, err;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",  
-                        &wc, &zc, &hc, 
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
+                        &wc, &zc, &hc,
                         &wv, &zv, &hv, &hvbar,
                         &xmomc, &ymomc, &xmomv, &ymomv))
@@ -719 +721 @@
 
   N = wc -> dimensions[0];
-    
+
   _balance_deep_and_shallow(N,
                             (double*) wc -> data,
-                            (double*) zc -> data, 
-                            (double*) hc -> data,                           
-                            (double*) wv -> data, 
-                            (double*) zv -> data, 
+                            (double*) zc -> data,
+                            (double*) hc -> data,
+                            (double*) wv -> data,
+                            (double*) zv -> data,
                             (double*) hv -> data,
-                            (double*) hvbar -> data,                        
-                            (double*) xmomc -> data, 
-                            (double*) ymomc -> data, 
-                            (double*) xmomv -> data, 
-                            (double*) ymomv -> data);  
-  
- 
-  return Py_BuildValue("");  
+                            (double*) hvbar -> data,
+                            (double*) xmomc -> data,
+                            (double*) ymomc -> data,
+                            (double*) xmomv -> data,
+                            (double*) ymomv -> data);
+
+
+  return Py_BuildValue("");
 }
 
@@ -740 +742 @@
 
 PyObject *h_limiter(PyObject *self, PyObject *args) {
-  
+
   PyObject *domain, *Tmp;
-  PyArrayObject 
-    *hv, *hc, //Depth at vertices and centroids        
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
     *hvbar,   //Limited depth at vertices (return values)
     *neighbours;
-    
+
   int k, i, n, N, k3;
   int dimensions[2];
   double beta_h; //Safety factor (see config.py)
   double *hmin, *hmax, hn;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) 
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
     return NULL;
 
@@ -759 +761 @@
 
   //Get safety factor beta_h
-  Tmp = PyObject_GetAttrString(domain, "beta_h"); 
-  if (!Tmp) 
-    return NULL;
-      
-  beta_h = PyFloat_AsDouble(Tmp);  
-  
-  Py_DECREF(Tmp);    
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
 
   N = hc -> dimensions[0];
-  
+
   //Create hvbar
   dimensions[0] = N;
-  dimensions[1] = 3;  
+  dimensions[1] = 3;
   hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
- 
+
 
   //Find min and max of this and neighbour's centroid values
   hmin = malloc(N * sizeof(double));
-  hmax = malloc(N * sizeof(double));  
+  hmax = malloc(N * sizeof(double));
   for (k=0; k<N; k++) {
     k3=k*3;
-    
+
     hmin[k] = ((double*) hc -> data)[k];
     hmax[k] = hmin[k];
-    
+
     for (i=0; i<3; i++) {
       n = ((long*) neighbours -> data)[k3+i];
-      
+
       //Initialise hvbar with values from hv
-      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];  
-      
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+
       if (n >= 0) {
         hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
-          
+
         hmin[k] = min(hmin[k], hn);
         hmax[k] = max(hmax[k], hn);
@@ -798 +800 @@
     }
   }
-  
+
   // Call underlying standard routine
   _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
-         
-  // // //Py_DECREF(domain); //FIXME: NEcessary?         
+
+  // // //Py_DECREF(domain); //FIXME: NEcessary?
   free(hmin);
-  free(hmax);  
-  
-  //return result using PyArray to avoid memory leak  
+  free(hmax);
+
+  //return result using PyArray to avoid memory leak
   return PyArray_Return(hvbar);
-  //return Py_BuildValue("");    
+  //return Py_BuildValue("");
 }
 
@@ -819 +821 @@
   //                              s_vec, phi_vec, self.const)
 
-  
+
 
   PyArrayObject   //(one element per triangle)
-  *s_vec,         //Speeds 
-  *phi_vec,       //Bearings 
+  *s_vec,         //Speeds
+  *phi_vec,       //Bearings
   *xmom_update,   //Momentum updates
-  *ymom_update; 
-
-    
-  int N; 
+  *ymom_update;
+
+
+  int N;
   double cw;
-  
-  // Convert Python arguments to C  
-  if (!PyArg_ParseTuple(args, "OOOOd",  
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOd",
                         &xmom_update,
-                        &ymom_update,                   
-                        &s_vec, &phi_vec, 
+                        &ymom_update,
+                        &s_vec, &phi_vec,
                         &cw))
     return NULL;
 
   N = xmom_update -> dimensions[0];
-    
+
   _assign_wind_field_values(N,
            (double*) xmom_update -> data,
-           (double*) ymom_update -> data, 
-           (double*) s_vec -> data, 
+           (double*) ymom_update -> data,
+           (double*) s_vec -> data,
            (double*) phi_vec -> data,
            cw);
- 
-  return Py_BuildValue("");  
-}
-
-
-
-
-//////////////////////////////////////////      
+
+  return Py_BuildValue("");
+}
+
+
+
+
+//////////////////////////////////////////
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
@@ -861 +863 @@
    * three.
    */
-  
+
   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
-  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},    
-  {"gravity", gravity, METH_VARARGS, "Print out"},      
-  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},        
-  {"balance_deep_and_shallow", balance_deep_and_shallow, 
-   METH_VARARGS, "Print out"},          
-  {"h_limiter", h_limiter, 
-   METH_VARARGS, "Print out"},             
-  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},   
-  {"assign_windfield_values", assign_windfield_values, 
-   METH_VARARGS | METH_KEYWORDS, "Print out"},     
-  //{"distribute_to_vertices_and_edges", 
-  // distribute_to_vertices_and_edges, METH_VARARGS},    
-  //{"update_conserved_quantities", 
-  // update_conserved_quantities, METH_VARARGS},       
-  //{"set_initialcondition", 
-  // set_initialcondition, METH_VARARGS},    
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {"gravity", gravity, METH_VARARGS, "Print out"},
+  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
+  {"balance_deep_and_shallow", balance_deep_and_shallow,
+   METH_VARARGS, "Print out"},
+  {"h_limiter", h_limiter,
+   METH_VARARGS, "Print out"},
+  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"assign_windfield_values", assign_windfield_values,
+   METH_VARARGS | METH_KEYWORDS, "Print out"},
+  //{"distribute_to_vertices_and_edges",
+  // distribute_to_vertices_and_edges, METH_VARARGS},
+  //{"update_conserved_quantities",
+  // update_conserved_quantities, METH_VARARGS},
+  //{"set_initialcondition",
+  // set_initialcondition, METH_VARARGS},
   {NULL, NULL}
 };
-        
-// Module initialisation   
+
+// Module initialisation
 void initshallow_water_ext(void){
   Py_InitModule("shallow_water_ext", MethodTable);
-  
-  import_array();     //Necessary for handling of NumPY structures  
-}
+
+  import_array();     //Necessary for handling of NumPY structures
+}

inundation/ga/storm_surge/pyvolution/util.py

@@ -2 +2 @@
 
 It is also a clearing house for functions that may later earn a module
-of their own. 
+of their own.
 """
 
@@ -17 +17 @@
     v2 = v[1]/l
 
-    try:        
+    try:
         theta = acos(v1)
     except ValueError, e:
@@ -31 +31 @@
             theta = 3.1415926535897931
         print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
-              %(v1, theta)   
-              
+              %(v1, theta)
+
     if v2 < 0:
         #Quadrant 3 or 4
@@ -43 +43 @@
     """Compute difference between angle of vector x0, y0 and x1, y1.
     This is used for determining the ordering of vertices,
-    e.g. for checking if they are counter clockwise. 
-    
+    e.g. for checking if they are counter clockwise.
+
     Always return a positive value
     """
-        
+
     from math import pi
-    
+
     a0 = angle(v0)
     a1 = angle(v1)
@@ -56 +56 @@
     if a0 < a1:
         a0 += 2*pi
-            
+
     return a0-a1
 
@@ -65 +65 @@
 
 
-def point_on_line(x, y, x0, y0, x1, y1): 
-    """Determine whether a point is on a line segment 
-   
+def point_on_line(x, y, x0, y0, x1, y1):
+    """Determine whether a point is on a line segment
+
     Input: x, y, x0, x0, x1, y1: where
         point is given by x, y
         line is given by (x0, y0) and (x1, y1)
-        
-    """
-        
+
+    """
+
     from Numeric import array, dot, allclose
     from math import sqrt
 
-    a = array([x - x0, y - y0]) 
+    a = array([x - x0, y - y0])
     a_normal = array([a[1], -a[0]])
-                        
+
     b = array([x1 - x0, y1 - y0])
-    
+
     if dot(a_normal, b) == 0:
         #Point is somewhere on the infinite extension of the line
 
-        len_a = sqrt(sum(a**2))                
-        len_b = sqrt(sum(b**2))                                 
+        len_a = sqrt(sum(a**2))
+        len_b = sqrt(sum(b**2))
         if dot(a, b) >= 0 and len_a <= len_b:
            return True
-        else:   
+        else:
            return False
     else:
-      return False  
-    
+      return False
+
 def ensure_numeric(A, typecode = None):
     """Ensure that sequence is a Numeric array.
@@ -117 +117 @@
         if type(A) == ArrayType:
             if A.typecode == typecode:
-                return array(A)            
+                return array(A)
             else:
                 return A.astype(typecode)
         else:
             return array(A).astype(typecode)
-        
-    
+
+
 
 def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False):
@@ -133 +133 @@
     #In fact, this is where origin should be converted to that of domain
     #Also, check that file covers domain fully.
-    
+
     assert type(filename) == type(''),\
                'First argument to File_function must be a string'
@@ -148 +148 @@
 
     if domain is not None:
-        quantities = domain.conserved_quantities   
+        quantities = domain.conserved_quantities
     else:
         quantities = None
-        
-    #Choose format 
+
+    #Choose format
     #FIXME: Maybe these can be merged later on
     if line[:3] == 'CDF':
         return File_function_NetCDF(filename, domain, quantities, interpolation_points, verbose = verbose)
     else:
-        return File_function_ASCII(filename, domain, quantities, interpolation_points)    
-           
-             
+        return File_function_ASCII(filename, domain, quantities, interpolation_points)
+
+
 class File_function_NetCDF:
-    """Read time history of spatial data from NetCDF sww file and 
-    return a callable object f(t,x,y)  
-    which will return interpolated values based on the input file.    
-    
+    """Read time history of spatial data from NetCDF sww file and
+    return a callable object f(t,x,y)
+    which will return interpolated values based on the input file.
+
     x, y may be either scalars or vectors
-    
+
     #FIXME: More about format, interpolation  and order of quantities
-    
+
     The quantities returned by the callable objects are specified by
-    the list quantities which must contain the names of the 
+    the list quantities which must contain the names of the
     quantities to be returned and also reflect the order, e.g. for
-    the shallow water wave equation, on would have 
+    the shallow water wave equation, on would have
     quantities = ['stage', 'xmomentum', 'ymomentum']
-    
-    interpolation_points decides at which points interpolated 
+
+    interpolation_points decides at which points interpolated
     quantities are to be computed whenever object is called.
-    
-    
-    
+
+
+
     If None, return average value
     """
-    
+
     def  __init__(self, filename, domain=None, quantities=None,
                   interpolation_points=None, verbose = False):
@@ -189 +189 @@
         If domain is specified, model time (domain.starttime)
         will be checked and possibly modified
-        
+
         All times are assumed to be in UTC
         """
@@ -196 +196 @@
         #(both in UTM coordinates)
         #If not - modify those from file to match domain
-        
+
         import time, calendar
         from config import time_format
-        from Scientific.IO.NetCDF import NetCDFFile     
+        from Scientific.IO.NetCDF import NetCDFFile
 
         #Open NetCDF file
@@ -213 +213 @@
              if not fid.variables.has_key(quantity):
                  missing.append(quantity)
-                 
+
         if len(missing) > 0:
             msg = 'Quantities %s could not be found in file %s'\
                   %(str(missing), filename)
             raise msg
-                 
+
 
         #Get first timestep
         try:
             self.starttime = fid.starttime[0]
-        except ValueError:    
+        except ValueError:
             msg = 'Could not read starttime from file %s' %filename
             raise msg
@@ -229 +229 @@
         #Get origin
         self.xllcorner = fid.xllcorner[0]
-        self.yllcorner = fid.yllcorner[0]            
-    
+        self.yllcorner = fid.yllcorner[0]
+
         self.number_of_values = len(quantities)
         self.fid = fid
@@ -249 +249 @@
                 if verbose: print 'Domain starttime is now set to %f' %domain.starttime
 
-        #Read all data in and produce values for desired data points at each timestep 
-        self.spatial_interpolation(interpolation_points, verbose = verbose) 
+        #Read all data in and produce values for desired data points at each timestep
+        self.spatial_interpolation(interpolation_points, verbose = verbose)
         fid.close()
 
     def spatial_interpolation(self, interpolation_points, verbose = False):
-        """For each timestep in fid: read surface, interpolate to desired points 
+        """For each timestep in fid: read surface, interpolate to desired points
         and store values for use when object is called.
         """
-        
+
         from Numeric import array, zeros, Float, alltrue, concatenate, reshape
 
@@ -263 +263 @@
         from least_squares import Interpolation
         import time, calendar
-        
+
         fid = self.fid
 
@@ -273 +273 @@
         triangles = fid.variables['volumes'][:]
         time = fid.variables['time'][:]
-        
+
         #Check
         msg = 'File %s must list time as a monotonuosly ' %self.filename
         msg += 'increasing sequence'
-        assert alltrue(time[1:] - time[:-1] > 0 ), msg  
-        
-
-        if interpolation_points is not None: 
-        
+        assert alltrue(time[1:] - time[:-1] > 0 ), msg
+
+
+        if interpolation_points is not None:
+
             try:
-                P = ensure_numeric(interpolation_points) 
+                P = ensure_numeric(interpolation_points)
             except:
                 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
                 msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
                 raise msg
-            
-            
+
+
             self.T = time[:]  #Time assumed to be relative to starttime
-            self.index = 0    #Initial time index       
-
-            self.values = {}                
+            self.index = 0    #Initial time index
+
+            self.values = {}
             for name in self.quantities:
-                self.values[name] = zeros( (len(self.T), 
-                                            len(interpolation_points)), 
-                                            Float)                  
+                self.values[name] = zeros( (len(self.T),
+                                            len(interpolation_points)),
+                                            Float)
 
             #Build interpolator
             if verbose: print 'Build interpolation matrix'
             x = reshape(x, (len(x),1))
-            y = reshape(y, (len(y),1))                
+            y = reshape(y, (len(y),1))
             vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
 
-            interpol = Interpolation(vertex_coordinates, 
+            interpol = Interpolation(vertex_coordinates,
                                      triangles,
                                      point_coordinates = P,
-                                     alpha = 0, 
+                                     alpha = 0,
                                      verbose = verbose)
 
@@ -318 +318 @@
                 for name in self.quantities:
                     self.values[name][i, :] =\
-                    interpol.interpolate(fid.variables[name][i,:])   
+                    interpol.interpolate(fid.variables[name][i,:])
 
             #Report
@@ -327 +327 @@
                 print '  Reference:'
                 print '    Lower left corner: [%f, %f]'\
-                      %(self.xllcorner, self.yllcorner) 
+                      %(self.xllcorner, self.yllcorner)
                 print '    Start time (file):   %f' %self.starttime
                 #print '    Start time (domain): %f' %domain.starttime
@@ -342 +342 @@
                     print '    %s in [%f, %f]' %(name, min(q), max(q))
 
-                    
+
                 print '  Interpolation points (xi, eta):'\
                       'number of points == %d ' %P.shape[0]
@@ -349 +349 @@
                 print '  Interpolated quantities (over all timesteps):'
                 for name in self.quantities:
-                    q = self.values[name][:].flat                    
+                    q = self.values[name][:].flat
                     print '    %s at interpolation points in [%f, %f]'\
-                          %(name, min(q), max(q))                
-                
-                
-                
+                          %(name, min(q), max(q))
+
+
+
 
                 print '------------------------------------------------'
@@ -367 +367 @@
     def __call__(self, t, x=None, y=None, point_id = None):
         """Evaluate f(t, point_id)
-        
+
         Inputs:
-          t: time - Model time (tau = domain.starttime-self.starttime+t)  
-                    must lie within existing timesteps   
-          point_id: index of one of the preprocessed points. 
+          t: time - Model time (tau = domain.starttime-self.starttime+t)
+                    must lie within existing timesteps
+          point_id: index of one of the preprocessed points.
                     If point_id is None all preprocessed points are computed
 
-          FIXME: point_id could also be a slice     
-          FIXME: One could allow arbitrary x, y coordinates 
+          FIXME: point_id could also be a slice
+          FIXME: One could allow arbitrary x, y coordinates
           (requires computation of a new interpolator)
           Maybe not,.,.
-          FIXME: What if x and y are vectors? 
+          FIXME: What if x and y are vectors?
         """
 
@@ -393 +393 @@
         #
         # domain.starttime + t == self.starttime + tau
-        
+
         if self.domain is not None:
             tau = self.domain.starttime-self.starttime+t
@@ -403 +403 @@
         #print 'S start', self.starttime
         #print 't', t
-        #print 'tau', tau             
-                
+        #print 'tau', tau
+
         msg = 'Time interval derived from file %s [%s:%s]'\
               %(self.filename, self.T[0], self.T[1])
@@ -411 +411 @@
 
         if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg        
+        if tau > self.T[-1]: raise msg
 
 
@@ -424 +424 @@
             #Protect against case where tau == T[-1] (last time)
             # - also works in general when tau == T[i]
-            ratio = 0 
+            ratio = 0
         else:
-            #t is now between index and index+1            
+            #t is now between index and index+1
             ratio = (tau - self.T[self.index])/\
                     (self.T[self.index+1] - self.T[self.index])
@@ -432 +432 @@
 
 
-                
+
         #Compute interpolated values
         q = zeros( len(self.quantities), Float)
 
         for i, name in enumerate(self.quantities):
-            Q = self.values[name]    
+            Q = self.values[name]
 
             if ratio > 0:
@@ -450 +450 @@
         #Return vector of interpolated values
         return q
-                
-             
+
+
 class File_function_ASCII:
     """Read time series from file and return a callable object:
@@ -468 +468 @@
     interpolated values, one each value stated in the file.
 
-    
+
     E.g
 
@@ -477 +477 @@
     ignoring x and y, which returns three values per call.
 
-    
+
     NOTE: At this stage the function is assumed to depend on
     time only, i.e  no spatial dependency!!!!!
     When that is needed we can use the least_squares interpolation.
-    
-    #FIXME: This should work with netcdf (e.g. sww) and thus render the 
+
+    #FIXME: This should work with netcdf (e.g. sww) and thus render the
     #spatio-temporal boundary condition in shallow water fully general
-    
-    #FIXME: Specified quantites not used here - 
+
+    #FIXME: Specified quantites not used here -
     #always return whatever is in the file
     """
 
-    
+
     def __init__(self, filename, domain=None, quantities = None, interpolation_points=None):
         """Initialise callable object from a file.
@@ -507 +507 @@
         line = fid.readline()
         fid.close()
-    
+
         fields = line.split(',')
         msg = 'File %s must have the format date, value0 value1 value2 ...'
@@ -516 +516 @@
         try:
             starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:    
+        except ValueError:
             msg = 'First field in file %s must be' %filename
             msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0] 
+            msg += 'I got %s instead.' %fields[0]
             raise msg
 
@@ -529 +529 @@
 
         q = ensure_numeric(values)
-        
+
         msg = 'ERROR: File must contain at least one independent value'
         assert len(q.shape) == 1, msg
-            
+
         self.number_of_values = len(q)
 
@@ -549 +549 @@
                 ##if verbose: print msg
                 domain.starttime = self.starttime #Modifying model time
-        
+
             #if domain.starttime is None:
             #    domain.starttime = self.starttime
@@ -563 +563 @@
             #        domain.starttime = self.starttime #Modifying model time
 
-        if mode == 2:        
+        if mode == 2:
             self.read_times() #Now read all times in.
         else:
             raise 'x,y dependency not yet implemented'
 
-        
+
     def read_times(self):
-        """Read time and values 
+        """Read time and values
         """
         from Numeric import zeros, Float, alltrue
         from config import time_format
         import time, calendar
-        
-        fid = open(self.filename)        
+
+        fid = open(self.filename)
         lines = fid.readlines()
         fid.close()
-        
+
         N = len(lines)
         d = self.number_of_values
@@ -585 +585 @@
         T = zeros(N, Float)       #Time
         Q = zeros((N, d), Float)  #Values
-        
+
         for i, line in enumerate(lines):
             fields = line.split(',')
@@ -603 +603 @@
         self.index = 0 #Initial index
 
-        
+
     def __repr__(self):
         return 'File function'
@@ -613 +613 @@
         result is a vector of same length as x and y
         FIXME: Naaaa
-        
+
         FIXME: Rethink semantics when x,y are vectors.
         """
 
         from math import pi, cos, sin, sqrt
-        
+
 
         #Find time tau such that
         #
         # domain.starttime + t == self.starttime + tau
-        
+
         if self.domain is not None:
             tau = self.domain.starttime-self.starttime+t
         else:
             tau = t
-        
-                
+
+
         msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
               %(self.filename, self.T[0], self.T[1], tau)
         if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg        
+        if tau > self.T[-1]: raise msg
 
         oldindex = self.index
@@ -647 +647 @@
         #Compute interpolated values
         q = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
 
         #Return vector of interpolated values
@@ -662 +662 @@
                 N = len(x)
                 assert len(y) == N, 'x and y must have same length'
-             
+
                 res = []
                 for col in q:
                     res.append(col*ones(N, Float))
-            
+
                 return res
-            
-            
-#FIXME: TEMP        
+
+
+#FIXME: TEMP
 class File_function_copy:
     """Read time series from file and return a callable object:
@@ -687 +687 @@
     interpolated values, one each value stated in the file.
 
-    
+
     E.g
 
@@ -696 +696 @@
     ignoring x and y, which returns three values per call.
 
-    
+
     NOTE: At this stage the function is assumed to depend on
     time only, i.e  no spatial dependency!!!!!
     When that is needed we can use the least_squares interpolation.
-    
-    #FIXME: This should work with netcdf (e.g. sww) and thus render the 
+
+    #FIXME: This should work with netcdf (e.g. sww) and thus render the
     #spatio-temporal boundary condition in shallow water fully general
     """
 
-    
+
     def __init__(self, filename, domain=None):
         """Initialise callable object from a file.
@@ -742 +742 @@
         try:
             starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:    
+        except ValueError:
             msg = 'First field in file %s must be' %filename
             msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0] 
+            msg += 'I got %s instead.' %fields[0]
             raise msg
 
@@ -755 +755 @@
 
         q = ensure_numeric(values)
-        
+
         msg = 'ERROR: File must contain at least one independent value'
         assert len(q.shape) == 1, msg
-            
+
         self.number_of_values = len(q)
 
@@ -779 +779 @@
                     domain.starttime = self.starttime #Modifying model time
 
-        if mode == 2:        
+        if mode == 2:
             self.read_times() #Now read all times in.
         else:
             raise 'x,y dependency not yet implemented'
 
-        
+
     def read_times(self):
-        """Read time and values 
+        """Read time and values
         """
         from Numeric import zeros, Float, alltrue
         from config import time_format
         import time, calendar
-        
-        fid = open(self.filename)        
+
+        fid = open(self.filename)
         lines = fid.readlines()
         fid.close()
-        
+
         N = len(lines)
         d = self.number_of_values
@@ -801 +801 @@
         T = zeros(N, Float)       #Time
         Q = zeros((N, d), Float)  #Values
-        
+
         for i, line in enumerate(lines):
             fields = line.split(',')
@@ -819 +819 @@
         self.index = 0 #Initial index
 
-        
+
     def __repr__(self):
         return 'File function'
 
-        
+
 
     def __call__(self, t, x=None, y=None):
@@ -834 +834 @@
 
         from math import pi, cos, sin, sqrt
-        
+
 
         #Find time tau such that
         #
         # domain.starttime + t == self.starttime + tau
-        
+
         if self.domain is not None:
             tau = self.domain.starttime-self.starttime+t
         else:
             tau = t
-        
-                
+
+
         msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
               %(self.filename, self.T[0], self.T[1], tau)
         if tau < self.T[0]: raise msg
-        if tau > self.T[-1]: raise msg        
+        if tau > self.T[-1]: raise msg
 
         oldindex = self.index
@@ -863 +863 @@
         #Compute interpolated values
         q = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:])
 
         #Return vector of interpolated values
@@ -878 +878 @@
                 N = len(x)
                 assert len(y) == N, 'x and y must have same length'
-             
+
                 res = []
                 for col in q:
                     res.append(col*ones(N, Float))
-            
+
                 return res
-            
+
 
 def read_xya(filename, format = 'netcdf'):
@@ -892 +892 @@
 
     Format can be either ASCII or NetCDF
-   
+
     Return list of points, list of attributes and
     attribute names if present otherwise None
     """
     #FIXME: Probably obsoleted by similar function in load_ASCII
-    
+
     from Scientific.IO.NetCDF import NetCDFFile
 
-    if format.lower() == 'netcdf': 
+    if format.lower() == 'netcdf':
         #Get NetCDF
-        
-        fid = NetCDFFile(filename, 'r') 
-      
+
+        fid = NetCDFFile(filename, 'r')
+
         # Get the variables
         points = fid.variables['points']
@@ -913 +913 @@
         #Don't close - arrays are needed outside this function,
         #alternatively take a copy here
-    else:    
+    else:
         #Read as ASCII file assuming that it is separated by whitespace
         fid = open(filename)
@@ -930 +930 @@
             attribute_names = ['elevation']  #HACK
 
-        attributes = {}    
+        attributes = {}
         for key in attribute_names:
             attributes[key] = []
@@ -940 +940 @@
             for i, key in enumerate(attribute_names):
                 attributes[key].append( float(fields[2+i]) )
-        
+
     return points, attributes
-    
+
 
 #####################################
@@ -950 +950 @@
 
 
-             
-#Redefine these to use separate_by_polygon which will put all inside indices 
+
+#Redefine these to use separate_by_polygon which will put all inside indices
 #in the first part of the indices array and outside indices in the last
 
@@ -957 +957 @@
     """See separate_points_by_polygon for documentation
     """
-    
-    from Numeric import array, Float, reshape    
-    
-    if verbose: print 'Checking input to inside_polygon'    
+
+    from Numeric import array, Float, reshape
+
+    if verbose: print 'Checking input to inside_polygon'
     try:
-        points = ensure_numeric(points, Float)             
+        points = ensure_numeric(points, Float)
     except:
         msg = 'Points could not be converted to Numeric array'
-        raise msg        
-        
+        raise msg
+
     try:
-        polygon = ensure_numeric(polygon, Float)           
+        polygon = ensure_numeric(polygon, Float)
     except:
         msg = 'Polygon could not be converted to Numeric array'
-        raise msg               
-    
-    
-    
+        raise msg
+
+
+
     if len(points.shape) == 1:
         one_point = True
         points = reshape(points, (1,2))
-    else:       
+    else:
         one_point = False
 
-    indices, count = separate_points_by_polygon(points, polygon, 
-                                                closed, verbose)     
-    
+    indices, count = separate_points_by_polygon(points, polygon,
+                                                closed, verbose)
+
     if one_point:
         return count == 1
     else:
-        return indices[:count]    
-        
+        return indices[:count]
+
 def outside_polygon(points, polygon, closed = True, verbose = False):
     """See separate_points_by_polygon for documentation
     """
-    
-    from Numeric import array, Float, reshape    
-    
-    if verbose: print 'Checking input to outside_polygon'    
+
+    from Numeric import array, Float, reshape
+
+    if verbose: print 'Checking input to outside_polygon'
     try:
-        points = ensure_numeric(points, Float)             
+        points = ensure_numeric(points, Float)
     except:
         msg = 'Points could not be converted to Numeric array'
-        raise msg        
-        
+        raise msg
+
     try:
-        polygon = ensure_numeric(polygon, Float)           
+        polygon = ensure_numeric(polygon, Float)
     except:
         msg = 'Polygon could not be converted to Numeric array'
-        raise msg               
-    
-    
-    
+        raise msg
+
+
+
     if len(points.shape) == 1:
         one_point = True
         points = reshape(points, (1,2))
-    else:       
+    else:
         one_point = False
 
-    indices, count = separate_points_by_polygon(points, polygon, 
-                                                closed, verbose)     
-    
+    indices, count = separate_points_by_polygon(points, polygon,
+                                                closed, verbose)
+
     if one_point:
         return count != 1
     else:
-        return indices[count:][::-1]  #return reversed    
-        
-
-def separate_points_by_polygon(points, polygon, 
+        return indices[count:][::-1]  #return reversed
+
+
+def separate_points_by_polygon(points, polygon,
                                closed = True, verbose = False):
     """Determine whether points are inside or outside a polygon
-    
+
     Input:
        points - Tuple of (x, y) coordinates, or list of tuples
        polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be 
-       regarded as belonging to the polygon (closed = True) 
-       or not (closed = False) 
-    
+       closed - (optional) determine whether points on boundary should be
+       regarded as belonging to the polygon (closed = True)
+       or not (closed = False)
+
     Outputs:
-       indices: array of same length as points with indices of points falling 
-       inside the polygon listed from the beginning and indices of points 
+       indices: array of same length as points with indices of points falling
+       inside the polygon listed from the beginning and indices of points
        falling outside listed from the end.
-       
+
        count: count of points falling inside the polygon
-       
+
        The indices of points inside are obtained as indices[:count]
-       The indices of points outside are obtained as indices[count:]       
-       
-       
+       The indices of points outside are obtained as indices[count:]
+
+
     Examples:
        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
-       
+
        separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
-       will return the indices [0, 2, 1] and count == 2 as only the first 
+       will return the indices [0, 2, 1] and count == 2 as only the first
        and the last point are inside the unit square
-       
+
     Remarks:
        The vertices may be listed clockwise or counterclockwise and
-       the first point may optionally be repeated.    
+       the first point may optionally be repeated.
        Polygons do not need to be convex.
-       Polygons can have holes in them and points inside a hole is 
-       regarded as being outside the polygon.          
-               
-    Algorithm is based on work by Darel Finley, 
-    http://www.alienryderflex.com/polygon/ 
-    
-    """    
+       Polygons can have holes in them and points inside a hole is
+       regarded as being outside the polygon.
+
+    Algorithm is based on work by Darel Finley,
+    http://www.alienryderflex.com/polygon/
+
+    """
 
     from Numeric import array, Float, reshape, Int, zeros
-    
-    
-    #Input checks
-    try:
-        points = ensure_numeric(points, Float)             
-    except:
-        msg = 'Points could not be converted to Numeric array'
-        raise msg        
-        
-    try:
-        polygon = ensure_numeric(polygon, Float)           
-    except:
-        msg = 'Polygon could not be converted to Numeric array'
-        raise msg               
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices'
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns'
-
-    assert len(points.shape) == 2,\
-       'Points array must be a 2d array'
-       
-    assert points.shape[1] == 2,\
-       'Points array must have two columns'
-       
-    N = polygon.shape[0] #Number of vertices in polygon 
-    M = points.shape[0]  #Number of points
-    
-    px = polygon[:,0]
-    py = polygon[:,1]   
-
-    #Used for an optimisation when points are far away from polygon
-    minpx = min(px); maxpx = max(px)
-    minpy = min(py); maxpy = max(py)    
-
-
-    #Begin main loop
-    indices = zeros(M, Int)
-    
-    inside_index = 0    #Keep track of points inside
-    outside_index = M-1 #Keep track of points outside (starting from end) 
-    
-    for k in range(M):
-
-        if verbose:
-            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
-        
-        #for each point   
-        x = points[k, 0]
-        y = points[k, 1]        
-
-        inside = False
-
-        if not x > maxpx or x < minpx or y > maxpy or y < minpy: 
-            #Check polygon 
-            for i in range(N):
-                j = (i+1)%N
-                
-                #Check for case where point is contained in line segment        
-                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
-                    if closed:
-                        inside = True
-                    else:        
-                        inside = False
-                    break
-                else:    
-                    #Check if truly inside polygon                  
-                    if py[i] < y and py[j] >= y or\
-                       py[j] < y and py[i] >= y:
-                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
-                            inside = not inside
-        
-        if inside: 
-            indices[inside_index] = k
-            inside_index += 1
-        else:
-            indices[outside_index] = k
-            outside_index -= 1    
-            
-    return indices, inside_index
-
-
-def separate_points_by_polygon_c(points, polygon, 
-                                 closed = True, verbose = False):
-    """Determine whether points are inside or outside a polygon
-
-    C-wrapper
-    """    
-
-    from Numeric import array, Float, reshape, zeros, Int
-    
-
-    if verbose: print 'Checking input to separate_points_by_polygon'    
+
+
     #Input checks
     try:
@@ -1169 +1076 @@
         raise msg
 
-    #if verbose: print 'Checking input to separate_points_by_polygon 2' 
     try:
         polygon = ensure_numeric(polygon, Float)
     except:
         msg = 'Polygon could not be converted to Numeric array'
-        raise msg               
-
-    if verbose: print 'check'
-    
+        raise msg
+
     assert len(polygon.shape) == 2,\
        'Polygon array must be a 2d array of vertices'
@@ -1186 +1090 @@
     assert len(points.shape) == 2,\
        'Points array must be a 2d array'
-       
+
     assert points.shape[1] == 2,\
        'Points array must have two columns'
-       
-    N = polygon.shape[0] #Number of vertices in polygon 
+
+    N = polygon.shape[0] #Number of vertices in polygon
     M = points.shape[0]  #Number of points
 
+    px = polygon[:,0]
+    py = polygon[:,1]
+
+    #Used for an optimisation when points are far away from polygon
+    minpx = min(px); maxpx = max(px)
+    minpy = min(py); maxpy = max(py)
+
+
+    #Begin main loop
+    indices = zeros(M, Int)
+
+    inside_index = 0    #Keep track of points inside
+    outside_index = M-1 #Keep track of points outside (starting from end)
+
+    for k in range(M):
+
+        if verbose:
+            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
+
+        #for each point
+        x = points[k, 0]
+        y = points[k, 1]
+
+        inside = False
+
+        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
+            #Check polygon
+            for i in range(N):
+                j = (i+1)%N
+
+                #Check for case where point is contained in line segment
+                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
+                    if closed:
+                        inside = True
+                    else:
+                        inside = False
+                    break
+                else:
+                    #Check if truly inside polygon
+                    if py[i] < y and py[j] >= y or\
+                       py[j] < y and py[i] >= y:
+                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
+                            inside = not inside
+
+        if inside:
+            indices[inside_index] = k
+            inside_index += 1
+        else:
+            indices[outside_index] = k
+            outside_index -= 1
+
+    return indices, inside_index
+
+
+def separate_points_by_polygon_c(points, polygon,
+                                 closed = True, verbose = False):
+    """Determine whether points are inside or outside a polygon
+
+    C-wrapper
+    """
+
+    from Numeric import array, Float, reshape, zeros, Int
+
+
+    if verbose: print 'Checking input to separate_points_by_polygon'
+    #Input checks
+    try:
+        points = ensure_numeric(points, Float)
+    except:
+        msg = 'Points could not be converted to Numeric array'
+        raise msg
+
+    #if verbose: print 'Checking input to separate_points_by_polygon 2'
+    try:
+        polygon = ensure_numeric(polygon, Float)
+    except:
+        msg = 'Polygon could not be converted to Numeric array'
+        raise msg
+
+    if verbose: print 'check'
+
+    assert len(polygon.shape) == 2,\
+       'Polygon array must be a 2d array of vertices'
+
+    assert polygon.shape[1] == 2,\
+       'Polygon array must have two columns'
+
+    assert len(points.shape) == 2,\
+       'Points array must be a 2d array'
+
+    assert points.shape[1] == 2,\
+       'Points array must have two columns'
+
+    N = polygon.shape[0] #Number of vertices in polygon
+    M = points.shape[0]  #Number of points
+
     from util_ext import separate_points_by_polygon
 
@@ -1199 +1199 @@
     indices = zeros( M, Int )
 
-    if verbose: print 'Calling C-version of inside poly' 
+    if verbose: print 'Calling C-version of inside poly'
     count = separate_points_by_polygon(points, polygon, indices,
                                        int(closed), int(verbose))
 
-    return indices, count 
+    return indices, count
 
 
 
 class Polygon_function:
-    """Create callable object f: x,y -> z, where a,y,z are vectors and 
-    where f will return different values depending on whether x,y belongs 
+    """Create callable object f: x,y -> z, where a,y,z are vectors and
+    where f will return different values depending on whether x,y belongs
     to specified polygons.
-    
+
     To instantiate:
-      
+
        Polygon_function(polygons)
-       
+
     where polygons is a list of tuples of the form
-    
+
       [ (P0, v0), (P1, v1), ...]
-      
-      with Pi being lists of vertices defining polygons and vi either 
+
+      with Pi being lists of vertices defining polygons and vi either
       constants or functions of x,y to be applied to points with the polygon.
-      
-    The function takes an optional argument, default which is the value 
+
+    The function takes an optional argument, default which is the value
     (or function) to used for points not belonging to any polygon.
     For example:
-    
-       Polygon_function(polygons, default = 0.03)       
-      
-    If omitted the default value will be 0.0  
-    
+
+       Polygon_function(polygons, default = 0.03)
+
+    If omitted the default value will be 0.0
+
     Note: If two polygons overlap, the one last in the list takes precedence
 
@@ -1236 +1236 @@
 
     """
-    
+
     def __init__(self, regions, default = 0.0):
-        
+
         try:
             len(regions)
         except:
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
+            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
             raise msg
 
 
-        T = regions[0]                                   
+        T = regions[0]
         try:
             a = len(T)
-        except:    
-            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
-            raise msg   
-                   
-        assert a == 2, 'Must have two component each: %s' %T       
-
-        self.regions = regions             
+        except:
+            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
+            raise msg
+
+        assert a == 2, 'Must have two component each: %s' %T
+
+        self.regions = regions
         self.default = default
-         
-         
+
+
     def __call__(self, x, y):
         from util import inside_polygon
         from Numeric import ones, Float, concatenate, array, reshape, choose
-        
+
         x = array(x).astype(Float)
-        y = array(y).astype(Float)      
-        
+        y = array(y).astype(Float)
+
         N = len(x)
         assert len(y) == N
-        
-        points = concatenate( (reshape(x, (N, 1)), 
+
+        points = concatenate( (reshape(x, (N, 1)),
                                reshape(y, (N, 1))), axis=1 )
 
         if callable(self.default):
             z = self.default(x,y)
-        else:                   
+        else:
             z = ones(N, Float) * self.default
-            
+
         for polygon, value in self.regions:
             indices = inside_polygon(points, polygon)
-            
-            #FIXME: This needs to be vectorised 
+
+            #FIXME: This needs to be vectorised
             if callable(value):
                 for i in indices:
@@ -1286 +1286 @@
                     yy = array([y[i]])
                     z[i] = value(xx, yy)[0]
-            else:    
+            else:
                 for i in indices:
-                    z[i] = value 
-
-        return z                 
+                    z[i] = value
+
+        return z
 
 def read_polygon(filename):
     """Read points assumed to form a polygon
        There must be exactly two numbers in each line
-    """              
-    
+    """
+
     #Get polygon
     fid = open(filename)
@@ -1306 +1306 @@
         polygon.append( [float(fields[0]), float(fields[1])] )
 
-    return polygon      
-    
+    return polygon
+
 def populate_polygon(polygon, number_of_points, seed = None):
     """Populate given polygon with uniformly distributed points.
@@ -1313 +1313 @@
     Input:
        polygon - list of vertices of polygon
-       number_of_points - (optional) number of points 
+       number_of_points - (optional) number of points
        seed - seed for random number generator (default=None)
-       
+
     Output:
        points - list of points inside polygon
-       
+
     Examples:
        populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
-       will return five randomly selected points inside the unit square 
+       will return five randomly selected points inside the unit square
     """
 
@@ -1334 +1334 @@
     max_y = min_y = polygon[0][1]
     for point in polygon[1:]:
-        x = point[0] 
+        x = point[0]
         if x > max_x: max_x = x
-        if x < min_x: min_x = x           
-        y = point[1] 
+        if x < min_x: min_x = x
+        y = point[1]
         if y > max_y: max_y = y
-        if y < min_y: min_y = y           
-
-
-    while len(points) < number_of_points:     
+        if y < min_y: min_y = y
+
+
+    while len(points) < number_of_points:
         x = uniform(min_x, max_x)
-        y = uniform(min_y, max_y)        
+        y = uniform(min_y, max_y)
 
         if inside_polygon( [x,y], polygon ):
@@ -1353 +1353 @@
 def tsh2sww(filename, verbose=True):
     """
-    to check if a tsh/msh file 'looks' good. 
+    to check if a tsh/msh file 'looks' good.
     """
     from shallow_water import Domain
     from pmesh2domain import pmesh_to_domain_instance
-    import time, os 
-    from data_manager import get_dataobject 
-    from os import sep, path  
+    import time, os
+    from data_manager import get_dataobject
+    from os import sep, path
     #from util import mean
-    
+
     if verbose == True:print 'Creating domain from', filename
     domain = pmesh_to_domain_instance(filename, Domain)
@@ -1378 +1378 @@
     if verbose == True:
         print "Output written to " + domain.get_datadir() + sep + \
-              domain.filename + "." + domain.format  
+              domain.filename + "." + domain.format
     sww = get_dataobject(domain)
     sww.store_connectivity()
@@ -1390 +1390 @@
     """
     """
-    
-    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
+
+    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
     a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
     a /= det
 
     b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
-    b /= det            
+    b /= det
 
     return a, b
@@ -1416 +1416 @@
     pass
 
-    
-    
-    
-    
-
-#OBSOLETED STUFF    
+
+
+
+
+
+#OBSOLETED STUFF
 def inside_polygon_old(point, polygon, closed = True, verbose = False):
     #FIXME Obsoleted
     """Determine whether points are inside or outside a polygon
-    
+
     Input:
        point - Tuple of (x, y) coordinates, or list of tuples
        polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be 
-       regarded as belonging to the polygon (closed = True) 
-       or not (closed = False) 
-    
+       closed - (optional) determine whether points on boundary should be
+       regarded as belonging to the polygon (closed = True)
+       or not (closed = False)
+
     Output:
        If one point is considered, True or False is returned.
-       If multiple points are passed in, the function returns indices 
-       of those points that fall inside the polygon     
-       
+       If multiple points are passed in, the function returns indices
+       of those points that fall inside the polygon
+
     Examples:
        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
        inside_polygon( [0.5, 0.5], U)
-       will evaluate to True as the point 0.5, 0.5 is inside the unit square 
-       
+       will evaluate to True as the point 0.5, 0.5 is inside the unit square
+
        inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
-       will return the indices [0, 2] as only the first and the last point 
+       will return the indices [0, 2] as only the first and the last point
        is inside the unit square
-       
+
     Remarks:
        The vertices may be listed clockwise or counterclockwise and
-       the first point may optionally be repeated.    
+       the first point may optionally be repeated.
        Polygons do not need to be convex.
-       Polygons can have holes in them and points inside a hole is 
-       regarded as being outside the polygon.          
-               
-       
-    Algorithm is based on work by Darel Finley, 
-    http://www.alienryderflex.com/polygon/ 
-    
-    """    
+       Polygons can have holes in them and points inside a hole is
+       regarded as being outside the polygon.
+
+
+    Algorithm is based on work by Darel Finley,
+    http://www.alienryderflex.com/polygon/
+
+    """
 
     from Numeric import array, Float, reshape
-    
-    
+
+
     #Input checks
     try:
-        point = array(point).astype(Float)         
+        point = array(point).astype(Float)
     except:
         msg = 'Point %s could not be converted to Numeric array' %point
-        raise msg        
-        
+        raise msg
+
     try:
-        polygon = array(polygon).astype(Float)             
+        polygon = array(polygon).astype(Float)
     except:
         msg = 'Polygon %s could not be converted to Numeric array' %polygon
-        raise msg               
-    
+        raise msg
+
     assert len(polygon.shape) == 2,\
        'Polygon array must be a 2d array of vertices: %s' %polygon
 
-        
+
     assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns: %s' %polygon    
+       'Polygon array must have two columns: %s' %polygon
 
 
@@ -1487 +1487 @@
         one_point = True
         points = reshape(point, (1,2))
-    else:       
+    else:
         one_point = False
         points = point
 
-    N = polygon.shape[0] #Number of vertices in polygon 
+    N = polygon.shape[0] #Number of vertices in polygon
     M = points.shape[0]  #Number of points
-    
+
     px = polygon[:,0]
-    py = polygon[:,1]   
+    py = polygon[:,1]
 
     #Used for an optimisation when points are far away from polygon
     minpx = min(px); maxpx = max(px)
-    minpy = min(py); maxpy = max(py)    
+    minpy = min(py); maxpy = max(py)
 
 
@@ -1508 +1508 @@
         if verbose:
             if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
-        
-        #for each point   
-        x = points[k, 0]
-        y = points[k, 1]        
+
+        #for each point
+        x = points[k, 0]
+        y = points[k, 1]
 
         inside = False
@@ -1517 +1517 @@
         #Optimisation
         if x > maxpx or x < minpx: continue
-        if y > maxpy or y < minpy: continue        
-
-        #Check polygon 
+        if y > maxpy or y < minpy: continue
+
+        #Check polygon
         for i in range(N):
             j = (i+1)%N
-            
-            #Check for case where point is contained in line segment    
+
+            #Check for case where point is contained in line segment
             if point_on_line(x, y, px[i], py[i], px[j], py[j]):
                 if closed:
                     inside = True
-                else:        
+                else:
                     inside = False
                 break
-            else:        
-                #Check if truly inside polygon              
+            else:
+                #Check if truly inside polygon
                 if py[i] < y and py[j] >= y or\
                    py[j] < y and py[i] >= y:
                     if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
                         inside = not inside
-                    
+
         if inside: indices.append(k)
 
-    if one_point: 
+    if one_point:
         return inside
     else:
         return indices
-             
+
 
 #def remove_from(A, B):
@@ -1552 +1552 @@
 #    ind = search_sorted(A, B)
 
-    
+
 
 def outside_polygon_old(point, polygon, closed = True, verbose = False):
     #OBSOLETED
     """Determine whether points are outside a polygon
-    
+
     Input:
        point - Tuple of (x, y) coordinates, or list of tuples
        polygon - list of vertices of polygon
-       closed - (optional) determine whether points on boundary should be 
-       regarded as belonging to the polygon (closed = True) 
-       or not (closed = False) 
-    
+       closed - (optional) determine whether points on boundary should be
+       regarded as belonging to the polygon (closed = True)
+       or not (closed = False)
+
     Output:
        If one point is considered, True or False is returned.
-       If multiple points are passed in, the function returns indices 
-       of those points that fall outside the polygon     
-       
+       If multiple points are passed in, the function returns indices
+       of those points that fall outside the polygon
+
     Examples:
-       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
        outside_polygon( [0.5, 0.5], U )
        will evaluate to False as the point 0.5, 0.5 is inside the unit square
@@ -1577 +1577 @@
        ouside_polygon( [1.5, 0.5], U )
        will evaluate to True as the point 1.5, 0.5 is outside the unit square
-       
+
        outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
-       will return the indices [1] as only the first and the last point 
+       will return the indices [1] as only the first and the last point
        is inside the unit square
     """
 
     #FIXME: This is too slow
-    
+
     res  = inside_polygon(point, polygon, closed, verbose)
 
@@ -1602 +1602 @@
 
     C-wrapper
-    """    
+    """
 
     from Numeric import array, Float, reshape, zeros, Int
-    
-
-    if verbose: print 'Checking input to inside_polygon'    
+
+
+    if verbose: print 'Checking input to inside_polygon'
     #Input checks
     try:
-        point = array(point).astype(Float)         
+        point = array(point).astype(Float)
     except:
         msg = 'Point %s could not be converted to Numeric array' %point
-        raise msg        
-        
+        raise msg
+
     try:
-        polygon = array(polygon).astype(Float)             
+        polygon = array(polygon).astype(Float)
     except:
         msg = 'Polygon %s could not be converted to Numeric array' %polygon
-        raise msg               
-    
+        raise msg
+
     assert len(polygon.shape) == 2,\
        'Polygon array must be a 2d array of vertices: %s' %polygon
 
-        
+
     assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns: %s' %polygon    
+       'Polygon array must have two columns: %s' %polygon
 
 
@@ -1632 +1632 @@
         one_point = True
         points = reshape(point, (1,2))
-    else:       
+    else:
         one_point = False
         points = point
@@ -1642 +1642 @@
     indices = zeros( points.shape[0], Int )
 
-    if verbose: print 'Calling C-version of inside poly' 
+    if verbose: print 'Calling C-version of inside poly'
     count = inside_polygon(points, polygon, indices,
                            int(closed), int(verbose))
@@ -1650 +1650 @@
     else:
         if verbose: print 'Got %d points' %count
-        
+
         return indices[:count]   #Return those indices that were inside
-

inundation/ga/storm_surge/zeus/anuga-workspace.zwi

@@ -2 +2 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>parallel</active>
+    <active>steve</active>
     <project name="Merimbula">Merimbula.zpi</project>
     <project name="parallel">parallel.zpi</project>

inundation/ga/storm_surge/zeus/parallel.zpi

@@ -74 +74 @@
     <ReleaseProjectReload>Off</ReleaseProjectReload>
     <file>..\parallel\advection.py</file>
+    <file>..\parallel\domain.py</file>
     <file>..\parallel\parallel_advection.py</file>
     <file>..\parallel\run_advection.py</file>

inundation/ga/storm_surge/zeus/steve.zpi

@@ -73 +73 @@
     <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
     <ReleaseProjectReload>Off</ReleaseProjectReload>
+    <file>..\steve\get_triangle_attribute.py</file>
     <file>..\steve\get_triangle_data.py</file>
     <file>..\steve\get_triangle_data_new.py</file>
+    <file>..\hobart\run_hobart_buildings.py</file>
     <file>..\pyvolution-parallel\vtk-pylab.py</file>
     <folder name="Header Files" />