Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 3689)
@@ -318,4 +318,6 @@
 
         See methods inside the quantity object for more options
+
+        FIXME: clean input args
         """
 
@@ -498,18 +500,4 @@
         print self.timestepping_statistics()
 
-        #Old version
-        #if self.min_timestep == self.max_timestep:
-        #    print 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
-        #          %(self.time, self.min_timestep, self.number_of_steps,
-        #            self.number_of_first_order_steps)
-        #elif self.min_timestep > self.max_timestep:
-        #    print 'Time = %.4f, steps=%d (%d)'\
-        #          %(self.time, self.number_of_steps,
-        #            self.number_of_first_order_steps)
-        #else:
-        #    print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
-        #          %(self.time, self.min_timestep,
-        #            self.max_timestep, self.number_of_steps,
-        #            self.number_of_first_order_steps)
 
     def timestepping_statistics(self):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 3689)
@@ -178,5 +178,5 @@
 
     def __repr__(self):
-        return 'Mesh: %d vertex coordinates, %d triangles'\
+        return 'Mesh: %d vertices, %d triangles'\
                %(self.coordinates.shape[0], len(self))
 
@@ -197,5 +197,5 @@
 
 
-    def get_vertex_coordinates(self, obj=False, absolute=False):
+    def get_vertex_coordinates(self, unique=False, obj=False, absolute=False):
         """Return all vertex coordinates.
         Return all vertex coordinates for all triangles as an Nx6 array
@@ -212,4 +212,14 @@
         (To see which, switch to default absolute=True and run tests).
         """
+
+        if unique is True:
+            V = self.coordinates
+            if absolute is True:
+                if not self.geo_reference.is_absolute():
+                    V = self.geo_reference.get_absolute(V)
+
+            return V
+
+                
 
         V = self.vertex_coordinates
@@ -226,5 +236,4 @@
                 
         if obj is True:
-
             N = V.shape[0]
             return reshape(V, (3*N, 2))
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 3689)
@@ -175,6 +175,6 @@
 
     def __repr__(self):
-        return 'Mesh: %d vertex_coordinates, %d triangles, %d boundary segments'\
-               %(self.coordinates.shape[0], len(self), len(self.boundary))
+        return General_mesh.__repr__(self) + ', %d boundary segments'\
+               %(len(self.boundary))
 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 3689)
@@ -15,4 +15,6 @@
 """
 
+from Numeric import array, zeros, Float, concatenate, NewAxis, argmax, allclose
+from anuga.utilities.numerical_tools import ensure_numeric
 
 class Quantity:
@@ -21,5 +23,4 @@
 
         from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
-        from Numeric import array, zeros, Float
 
         msg = 'First argument in Quantity.__init__ '
@@ -619,7 +620,4 @@
 
 
-        from Numeric import Float
-        from anuga.utilities.numerical_tools import ensure_numeric
-        #from anuga.abstract_2d_finite_volumes.least_squares import fit_to_mesh
         from anuga.fit_interpolate.fit import fit_to_mesh
         from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -754,19 +752,136 @@
 
 
-    def get_values(self, location='vertices', indices = None):
+   
+    
+    def get_maximum_index(self, indices=None):
+        """Return index for maximum value of quantity (on centroids)
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            i = get_maximum_index()
+
+        Notes:
+            We do not seek the maximum at vertices as each vertex can
+            have multiple values - one for each triangle sharing it.
+
+            If there are multiple cells with same maximum value, the first cell
+            encountered in the triangle array is returned.
+        """
+
+        V = self.get_values(location='centroids', indices=indices)
+
+        # Always return absolute indices
+        i = argmax(V)
+
+        if indices is None:
+            return i
+        else:
+            return indices[i]
+
+        
+    def get_maximum_value(self, indices=None):
+        """Return maximum value of quantity (on centroids)
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            v = get_maximum_value()
+
+        Note, we do not seek the maximum at vertices as each vertex can
+        have multiple values - one for each triangle sharing it            
+        """
+
+
+        i = self.get_maximum_index(indices)
+        V = self.get_values(location='centroids') #, indices=indices)
+        
+        return V[i]
+        
+
+    def get_maximum_location(self, indices=None):
+        """Return location of maximum value of quantity (on centroids)
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            x, y = get_maximum_location()
+
+
+        Notes:
+            We do not seek the maximum at vertices as each vertex can
+            have multiple values - one for each triangle sharing it.
+
+            If there are multiple cells with same maximum value, the first cell
+            encountered in the triangle array is returned.            
+        """
+
+        i = self.get_maximum_index(indices)
+        x, y = self.domain.get_centroid_coordinates()[i]
+
+        return x, y
+
+
+
+
+    def get_interpolated_values(self, interpolation_points):
+
+        # Interpolation object based on internal (discontinuous triangles)
+        x, y, vertex_values, triangles = self.get_vertex_values(xy=True, smooth=False)
+        # FIXME: This concat should roll into get_vertex_values
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+
+        can_reuse = False
+        if hasattr(self, 'interpolation_object'):
+            # Reuse to save time
+            I = self.interpolation_object
+
+            if allclose(interpolation_points, I._point_coordinates):
+                can_reuse = True
+                
+
+        if can_reuse is True:        
+            result = I.interpolate(vertex_values) # Use absence to indicate reuse
+        else:    
+            from anuga.fit_interpolate.interpolate import Interpolate
+
+            # Create interpolation object with matrix
+            I = Interpolate(vertex_coordinates, triangles)
+            self.interpolation_object = I
+
+            # Call interpolate with points the first time
+            interpolation_points = ensure_numeric(interpolation_points, Float)                        
+            result = I.interpolate(vertex_values, interpolation_points)            
+
+        return result
+
+
+    def get_values(self, interpolation_points=None, location='vertices', indices = None):
         """get values for quantity
 
         return X, Compatible list, Numeric array (see below)
+        interpolation_points: List of x, y coordinates where value is sought (using interpolation)
+                If points are given, values of location and indices are ignored
         location: Where values are to be stored.
                   Permissible options are: vertices, edges, centroid
                   and unique vertices. Default is 'vertices'
 
-        In case of location == 'centroids' the dimension values must
-        be a list of a Numerical array of length N, N being the number
-        of elements. Otherwise it must be of dimension Nx3
 
         The returned values with be a list the length of indices
-        (N if indices = None).  Each value will be a list of the three
-        vertex values for this quantity.
+        (N if indices = None).
+
+        In case of location == 'centroids' the dimension of returned values will
+        be a list or a Numerical array of length N, N being the number
+        of elements.
+        
+        In case of location == 'vertices' or 'edges' the dimension of returned values will
+        be of dimension Nx3
+
+        In case of location == 'unique vertices' the average value at each vertex will be
+        returned and the dimension of returned values will be a 1d array of length "number of vertices" 
+        
 
         Indices is the set of element ids that the operation applies to.
@@ -777,4 +892,9 @@
         """
         from Numeric import take
+
+        if interpolation_points is not None:
+            return self.get_interpolated_values(interpolation_points)
+        
+        
 
         if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
@@ -810,4 +930,7 @@
                 # Go through all triangle, vertex pairs
                 # Average the values
+                
+                # FIXME (Ole): Should we merge this with get_vertex_values
+                # and use the concept of a reduction operator?
                 sum = 0
                 for triangle_id, vertex_id in triangles:
@@ -940,10 +1063,10 @@
             precision = Float
 
-        if reduction is None:
-            reduction = self.domain.reduction
-
         #Create connectivity
 
         if smooth == True:
+            
+            if reduction is None:
+                reduction = self.domain.reduction
 
             V = self.domain.get_vertices()
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 3689)
@@ -18,4 +18,20 @@
     def tearDown(self):
         pass
+
+
+    def test_get_vertex_coordinates(self):
+        from mesh_factory import rectangular
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        domain = General_mesh(points, vertices, boundary)
+
+        assert allclose(domain.get_vertex_coordinates(unique=True), domain.coordinates)
+
+        #assert allclose(domain.get_vertex_coordinates(), ...TODO
+        #assert allclose(domain.get_vertex_coordinates(absolute=True), ...TODO
+        
+        
 
     def test_get_vertex_values(self):
@@ -67,4 +83,7 @@
         assert unique_vertices == [0,2,4,5,6,7]
 
+
+        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 3688)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 3689)
@@ -110,4 +110,89 @@
                                                [4.5, 4.5, 0.],
                                                [3.0, -1.5, -1.5]])
+
+    def test_get_maximum_1(self):
+        quantity = Conserved_quantity(self.mesh4,
+                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
+        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
+
+        v = quantity.get_maximum_value()
+        assert v == 5
+
+        i = quantity.get_maximum_index()
+        assert i == 1
+        
+        x,y = quantity.get_maximum_location()
+        xref, yref = 4.0/3, 4.0/3
+        assert x == xref
+        assert y == yref
+
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 5)
+
+    def test_get_maximum_2(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        quantity = Quantity(domain)
+        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
+        
+        v = quantity.get_maximum_value()
+        assert v == 6
+
+        i = quantity.get_maximum_index()
+        assert i == 3
+        
+        x,y = quantity.get_maximum_location()
+        xref, yref = 2.0/3, 8.0/3
+        assert x == xref
+        assert y == yref
+
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 6)
+
+
+        #Multiple locations for maximum -
+        #Test that the algorithm picks the first occurrence        
+        v = quantity.get_maximum_value(indices=[0,1,2])
+        assert allclose(v, 4)
+
+        i = quantity.get_maximum_index(indices=[0,1,2])
+        assert i == 1
+        
+        x,y = quantity.get_maximum_location(indices=[0,1,2])
+        xref, yref = 4.0/3, 4.0/3
+        assert x == xref
+        assert y == yref
+
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 4)        
+
+        # More test of indices......
+        v = quantity.get_maximum_value(indices=[2,3])
+        assert allclose(v, 6)
+
+        i = quantity.get_maximum_index(indices=[2,3])
+        assert i == 3
+        
+        x,y = quantity.get_maximum_location(indices=[2,3])
+        xref, yref = 2.0/3, 8.0/3
+        assert x == xref
+        assert y == yref
+
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 6)        
+
+        
 
     def test_boundary_allocation(self):
@@ -1389,4 +1474,142 @@
         #      quantity.get_values(location = 'centroids')
 
+
+
+
+    def test_get_values_2(self):
+        """Different mesh (working with domain object) - also check centroids.
+        """
+
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        quantity = Quantity(domain)
+        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
+        
+        assert allclose(quantity.get_values(location='centroids'), [2,4,4,6])
+        assert allclose(quantity.get_values(location='centroids', indices=[1,3]), [4,6])
+
+
+        assert allclose(quantity.get_values(location='vertices'), [[4,0,2],
+                                                                   [4,2,6],
+                                                                   [6,2,4],
+                                                                   [8,4,6]])
+        
+        assert allclose(quantity.get_values(location='vertices', indices=[1,3]), [[4,2,6],
+                                                                                  [8,4,6]])
+
+
+        assert allclose(quantity.get_values(location='edges'), [[1,3,2],
+                                                                [4,5,3],
+                                                                [3,5,4],
+                                                                [5,7,6]])
+        assert allclose(quantity.get_values(location='edges', indices=[1,3]),
+                        [[4,5,3],
+                         [5,7,6]])        
+
+        # Check averaging over vertices
+        #a: 0
+        #b: (4+4+4)/3
+        #c: (2+2+2)/3
+        #d: 8
+        #e: (6+6+6)/3        
+        #f: 4
+        assert(quantity.get_values(location='unique vertices'), [0, 4, 2, 8, 6, 4])        
+                                                                                  
+        
+
+
+
+
+    def test_get_interpolated_values(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        domain = Domain(points, vertices, boundary)
+
+        #Constant values
+        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
+                                    [4,4,4],[5,5,5]])
+
+        
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        
+        #print quantity.get_values(points=interpolation_points)
+        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))
+
+
+        #Arbitrary values
+        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
+                                    [1,4,-9],[2,5,0]])
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+        #print answer
+        #print quantity.get_values(points=interpolation_points)
+        assert allclose(answer, quantity.get_values(interpolation_points=interpolation_points))        
+                        
+
+        #FIXME TODO
+        #indices = [0,5,3]
+        #answer = [0.5,1,5]
+        #assert allclose(answer,
+        #                quantity.get_values(indices=indices, \
+        #                                    location = 'unique vertices'))
+
+
+
+
+    def test_get_interpolated_values_2(self):
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        quantity = Quantity(domain)
+        quantity.set_values(lambda x, y: x+2*y) #2 4 4 6
+
+        #First pick one point
+        x, y = 2.0/3, 8.0/3
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 6)        
+
+        # Then another to test that algorithm won't blindly
+        # reuse interpolation matrix
+        x, y = 4.0/3, 4.0/3
+        v = quantity.get_values(interpolation_points = [[x,y]])
+        assert allclose(v, 4)        
+
+
+
+
     def test_getting_some_vertex_values(self):
         """
Index: /anuga_core/source/anuga/config.py
===================================================================
--- /anuga_core/source/anuga/config.py	(revision 3688)
+++ /anuga_core/source/anuga/config.py	(revision 3689)
@@ -117,3 +117,3 @@
 maximum_allowed_speed = 100.0 #Maximal particle speed of water
 
-minimum_storable_height = 0.001 #Water depth below which it is *stored* as 0
+minimum_storable_height = minimum_allowed_height #Water depth below which it is *stored* as 0
Index: /anuga_core/source/anuga/damage_modelling/test_inundation_damage.py
===================================================================
--- /anuga_core/source/anuga/damage_modelling/test_inundation_damage.py	(revision 3688)
+++ /anuga_core/source/anuga/damage_modelling/test_inundation_damage.py	(revision 3689)
@@ -142,4 +142,6 @@
     def tearDown(self):
         #print "***** tearDown  ********"
+
+        # FIXME (Ole): Sometimes this fails - is the file open or is it sometimes not created?
         os.remove(self.sww.filename)
         os.remove(self.csv_file)
Index: /anuga_core/source/anuga/examples/netherlands.py
===================================================================
--- /anuga_core/source/anuga/examples/netherlands.py	(revision 3688)
+++ /anuga_core/source/anuga/examples/netherlands.py	(revision 3689)
@@ -166,6 +166,5 @@
 #
 print 'Initial condition'
-domain.set_quantity('stage', Constant_height(Z, 0.0))
-#domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
+domain.set_quantity('stage', expression='Z + 0.0')
 
 #Evolve
Index: /anuga_core/source/anuga/examples/runup.py
===================================================================
--- /anuga_core/source/anuga/examples/runup.py	(revision 3688)
+++ /anuga_core/source/anuga/examples/runup.py	(revision 3689)
@@ -9,5 +9,4 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
-
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 from anuga.shallow_water import Domain
@@ -15,5 +14,5 @@
 from anuga.shallow_water import Dirichlet_boundary
 from anuga.shallow_water import Time_boundary
-from anuga.shallow_water import Transmissive_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
 
 
@@ -21,6 +20,6 @@
 # Setup computational domain
 #------------------------------------------------------------------------------
-
-points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
+N=10
+points, vertices, boundary = rectangular_cross(N, N) # Basic mesh
 
 domain = Domain(points, vertices, boundary) # Create domain
@@ -31,11 +30,10 @@
 # Setup initial conditions
 #------------------------------------------------------------------------------
-
 def topography(x,y):
     return -x/2                             # linear bed slope
-    #return x*(-(2.0-x)*.5)                  # curved bed slope
+    
 
 domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0.1)         # Constant friction 
+domain.set_quantity('friction', 0.0)         # Constant friction 
 domain.set_quantity('stage', -.4)            # Constant negative initial stage
 
@@ -44,11 +42,12 @@
 # Setup boundary conditions
 #------------------------------------------------------------------------------
+from math import sin, pi, exp
 
-from math import sin, pi, exp
+def waveform(t): 
+    return (0.1*sin(t*2*pi)-0.8) * exp(-2*t)
+
+Bw = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
 Br = Reflective_boundary(domain)      # Solid reflective wall
-Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
-Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
-Bw = Time_boundary(domain=domain,     # Time dependent boundary  
-                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-2*t), 0.0, 0.0])
+Bd = Dirichlet_boundary([-0.3,0,0])
 
 # Associate boundary tags with boundary objects
@@ -59,6 +58,5 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-
-for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0):
+for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
     domain.write_time()
     
Index: /anuga_core/source/anuga/examples/runup_convergence.py
===================================================================
--- /anuga_core/source/anuga/examples/runup_convergence.py	(revision 3688)
+++ /anuga_core/source/anuga/examples/runup_convergence.py	(revision 3689)
@@ -34,5 +34,5 @@
 sea_level = 0       # Mean sea level
 min_elevation = -20 # Lowest elevation (elevation of offshore flat part)
-offshore_depth=sea_level-min_elevation # offshore water depth
+offshore_depth = sea_level-min_elevation # offshore water depth
 amplitude = 0.5       # Solitary wave height H
 normalized_amplitude = amplitude/offshore_depth 
@@ -52,6 +52,6 @@
 
 # Structured mesh
-dx = 20           # Resolution: Length of subdivisions on x axis (length)
-dy = 20           # Resolution: Length of subdivisions on y axis (width)
+dx = 30           # Resolution: Length of subdivisions on x axis (length)
+dy = 30           # Resolution: Length of subdivisions on y axis (width)
 
 length = east-west
@@ -74,5 +74,5 @@
                          interior_regions=[[interior_polygon,dx*dy/32]])
 domain = Domain(meshname, use_cache=True, verbose = True)
-
+domain.set_minimum_storable_height(0.01)
 
 
@@ -134,11 +134,30 @@
 
 
+w0 = domain.get_maximum_inundation_elevation()
+x0, y0 = domain.get_maximum_inundation_location()
+print 'Coastline elevation = %.2f at (%.2f, %.2f)' %(w0, x0, y0)
+w_i = domain.get_quantity('stage').get_values(interpolation_points=[[x0,y0]])
+print 'Interpolated elevation at (%.2f, %.2f) is %.2f' %(x0, y0, w_i) 
+
 #------------------------------------------------------------------------------
 # Evolve system through time
 #------------------------------------------------------------------------------
 
+w_max = w0
 stagestep = []
 for t in domain.evolve(yieldstep = 1, finaltime = 300):
     domain.write_time()
+
+    w = domain.get_maximum_inundation_elevation()
+    x, y = domain.get_maximum_inundation_location()
+    print '  Coastline elevation = %.2f at (%.2f, %.2f)' %(w, x, y)    
+    #w_i = domain.get_quantity('stage').get_values(interpolation_points=[[x,y]])
+    #print '  Interpolated elevation at (%.2f, %.2f) is %.2f' %(x, y, w_i)
+                                                             
+
+    if w > w_max:
+        w_max = w
+        x_max = x
+        y_max = y
 
     # Let's find the maximum runup here working directly with the quantities,
@@ -149,5 +168,10 @@
     # 3 Find min x where h>0 over all t.
     # 4 Perhaps do this across a number of ys
-    
+
+print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
+
+print 'Run up distance - %.2f' %sqrt( (x-x0)**2 + (y-y0)**2 )
+
+
 
 #-----------------------------------------------------------------------------
Index: /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 3688)
+++ /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 3689)
@@ -77,5 +77,5 @@
         triangles = ensure_numeric(triangles, Int)
         vertex_coordinates = ensure_absolute(vertex_coordinates,
-                                         geo_reference = mesh_origin)
+                                             geo_reference = mesh_origin)
 
         #Don't pass geo_reference to mesh.  It doesn't work.
@@ -83,6 +83,8 @@
         self.mesh.check_integrity()
         self.root = build_quadtree(self.mesh,
-                              max_points_per_cell = max_vertices_per_cell)
+                                   max_points_per_cell = max_vertices_per_cell)
         #print "self.root",self.root.show() 
         
         
+    def __repr__(self):
+        return 'Interpolation object based on: ' + repr(self.mesh)
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 3688)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 3689)
@@ -76,4 +76,6 @@
               a mesh origin, since geospatial has its own mesh origin.
         """
+
+        # FIXME (Ole): Need an input check
         
         # Initialise variabels
@@ -82,11 +84,11 @@
         
         FitInterpolate.__init__(self,
-                 vertex_coordinates,
-                 triangles,
-                 mesh_origin,
-                 verbose,
-                 max_vertices_per_cell)
-
-     # FIXME: What is a good start_blocking_len value?
+                                vertex_coordinates,
+                                triangles,
+                                mesh_origin,
+                                verbose,
+                                max_vertices_per_cell)
+
+    # FIXME: What is a good start_blocking_len value?
     def interpolate(self, f, point_coordinates = None,
                     start_blocking_len = 500000, verbose=False):
@@ -114,13 +116,20 @@
 	  Interpolated values at inputted points (z).
         """
-        #print "point_coordinates interpolate.interpolate",point_coordinates 
-        if isinstance(point_coordinates,Geospatial_data):
+
+        
+        # FIXME (Ole): Need an input check that dimensions are compatible
+
+        # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the method is called
+        # even if interpolation points are unchanged.
+        
+        #print "point_coordinates interpolate.interpolate", point_coordinates 
+        if isinstance(point_coordinates, Geospatial_data):
             point_coordinates = point_coordinates.get_data_points( \
                 absolute = True)
-        
+
         # Can I interpolate, based on previous point_coordinates?
         if point_coordinates is None:
             if self._A_can_be_reused is True and \
-            len(self._point_coordinates) < start_blocking_len: 
+                   len(self._point_coordinates) < start_blocking_len:
                 z = self._get_point_data_z(f,
                                            verbose=verbose)
@@ -134,7 +143,6 @@
                 msg = 'ERROR (interpolate.py): No point_coordinates inputted'
                 raise Exception(msg)
-            
-            
-        if point_coordinates is not None:
+        
+        if point_coordinates is not None:   
             self._point_coordinates = point_coordinates
             if len(point_coordinates) < start_blocking_len or \
@@ -144,7 +152,8 @@
                                            verbose=verbose)
             else:
+                #print 'BLOCKING'
                 #Handle blocking
                 self._A_can_be_reused = False
-                start=0
+                start = 0
                 # creating a dummy array to concatenate to.
                 
@@ -156,7 +165,7 @@
                     z = zeros((0,))
                     
-                for end in range(start_blocking_len
-                                 ,len(point_coordinates)
-                                 ,start_blocking_len):
+                for end in range(start_blocking_len,
+                                 len(point_coordinates),
+                                 start_blocking_len):
                     t = self.interpolate_block(f, point_coordinates[start:end],
                                                verbose=verbose)
@@ -171,4 +180,5 @@
         return z
 
+
     def interpolate_block(self, f, point_coordinates = None, verbose=False):
         """
@@ -184,6 +194,6 @@
                 absolute = True)
         if point_coordinates is not None:
-            self._A =self._build_interpolation_matrix_A(point_coordinates,
-                                                       verbose=verbose)
+            self._A = self._build_interpolation_matrix_A(point_coordinates,
+                                                         verbose=verbose)
         return self._get_point_data_z(f)
 
@@ -206,6 +216,6 @@
 
     def _build_interpolation_matrix_A(self,
-                                     point_coordinates,
-                                     verbose = False):
+                                      point_coordinates,
+                                      verbose = False):
         """Build n x m interpolation matrix, where
         n is the number of data points and
@@ -224,5 +234,5 @@
         """
 
-
+        #print 'Building interpolation matrix'
         
         #Convert point_coordinates to Numeric arrays, in case it was a list.
@@ -486,9 +496,9 @@
             #Build interpolator
             interpol = Interpolate(vertex_coordinates,
-                                     triangles,
-                                     #point_coordinates = \
-                                     #self.interpolation_points,
-                                     #alpha = 0,
-                                     verbose = verbose)
+                                   triangles,
+                                   #point_coordinates = \
+                                   #self.interpolation_points,
+                                   #alpha = 0,
+                                   verbose = verbose)
 
             if verbose: print 'Interpolate'
@@ -501,6 +511,6 @@
                     if len(quantities[name].shape) == 2:
                         result = interpol.interpolate(quantities[name][i,:],
-                                     point_coordinates = \
-                                     self.interpolation_points)
+                                                      point_coordinates = \
+                                                      self.interpolation_points)
                     else:
                        #Assume no time dependency 
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 3688)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 3689)
@@ -105,4 +105,60 @@
         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
                         [[1./3, 1./3, 1./3]])
+
+
+
+    def test_simple_interpolation_example(self):
+        
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        from abstract_2d_finite_volumes.quantity import Quantity
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+
+        #---------------
+        #Constant values
+        #---------------        
+        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
+                                    [4,4,4],[5,5,5]])
+
+
+        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+        # FIXME: This concat should roll into get_vertex_values
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        I = Interpolate(vertex_coordinates, triangles)
+        result = I.interpolate(vertex_values, interpolation_points)
+        assert allclose(result, answer)
+
+
+        #---------------
+        #Variable values
+        #---------------
+        quantity = Quantity(domain,[[0,1,2],[3,1,7],[2,1,2],[3,3,7],
+                                    [1,4,-9],[2,5,0]])
+        
+        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
+        vertex_coordinates = concatenate( (x[:, NewAxis], y[:, NewAxis]), axis=1 )
+        # FIXME: This concat should roll into get_vertex_values
+
+
+        # Get interpolated values at centroids
+        interpolation_points = domain.get_centroid_coordinates()
+        answer = quantity.get_values(location='centroids')
+
+        I = Interpolate(vertex_coordinates, triangles)
+        result = I.interpolate(vertex_values, interpolation_points)
+        assert allclose(result, answer)        
+        
 
     def test_quad_tree(self):
@@ -770,5 +826,5 @@
 
         
-    def test_interpolate_reuse(self):
+    def test_interpolate_reuse_if_None(self):
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -827,4 +883,51 @@
         except:
             pass
+        
+    def xxtest_interpolate_reuse_if_same(self):
+
+        # This on tests that repeated identical interpolation
+        # points makes use of precomputed matrix (Ole)
+        # This is not really a test and is disabled for now
+        
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0, 1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0, -2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
+
+
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0, 2.0],
+                        [1.0, 1.0],
+                        [2.0, 1.0],
+                        [0.0, 0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2, -0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0, 1.0]]
+
+        interp = Interpolate(vertices, triangles)
+        f = array([linear_function(vertices),2*linear_function(vertices) ])
+        f = transpose(f)
+        z = interp.interpolate(f, point_coords)
+        answer = [linear_function(point_coords),
+                  2*linear_function(point_coords) ]
+        answer = transpose(answer)
+
+        assert allclose(z, answer)
+        assert allclose(interp._A_can_be_reused, True)
+
+
+        z = interp.interpolate(f)    # None
+        assert allclose(z, answer)        
+        z = interp.interpolate(f, point_coords) # Repeated (not really a test)        
+        assert allclose(z, answer)
         
 
@@ -1469,5 +1572,5 @@
 
     suite = unittest.makeSuite(Test_Interpolate,'test')
-    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolate_sww2csv')
+    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_function_outside_point')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)
Index: /anuga_core/source/anuga/pmesh/mesh_interface.py
===================================================================
--- /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 3688)
+++ /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 3689)
@@ -28,4 +28,5 @@
                              mesh_geo_reference=None,
                              minimum_triangle_angle=28.0,
+                             use_cache=False,
                              verbose=True):
     """Create mesh from bounding polygons, and resolutions.
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 3688)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 3689)
@@ -70,4 +70,7 @@
 #$LastChangedRevision$
 #$LastChangedBy$
+
+from Numeric import zeros, ones, Float, array, sum, size
+from Numeric import compress, arange
 
 
@@ -217,4 +220,65 @@
 
 
+
+    def get_wet_elements(self, indices=None):
+        """Return indices for elements where h > minimum_allowed_height
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            indices = get_wet_elements()
+
+        Note, centroid values are used for this operation            
+        """
+
+        # Water depth below which it is considered to be 0 in the model
+        # FIXME (Ole): Allow this to be specified as a keyword argument as well
+        from anuga.config import minimum_allowed_height
+        
+
+        elevation = self.get_quantity('elevation').get_values(location='centroids', indices=indices)
+        stage = self.get_quantity('stage').get_values(location='centroids', indices=indices)                
+        depth = stage - elevation
+
+        # Select indices for which depth > 0
+        wet_indices = compress(depth > minimum_allowed_height, arange(len(depth)))
+        return wet_indices 
+
+
+    def get_maximum_inundation_elevation(self, indices=None):
+        """Return highest elevation where h > 0
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            q = get_maximum_inundation_elevation()
+
+        Note, centroid values are used for this operation            
+        """
+
+        wet_elements = self.get_wet_elements(indices)
+        return self.get_quantity('elevation').get_maximum_value(indices=wet_elements)
+
+
+    def get_maximum_inundation_location(self, indices=None):
+        """Return highest elevation where h > 0
+
+        Optional argument:
+            indices is the set of element ids that the operation applies to.
+
+        Usage:
+            q = get_maximum_inundation_elevation()
+
+        Note, centroid values are used for this operation            
+        """
+
+        wet_elements = self.get_wet_elements(indices)
+        return self.get_quantity('elevation').get_maximum_location(indices=wet_elements)    
+
+
+
+
     def initialise_visualiser(self,scale_z=1.0,rect=None):
         #Realtime visualisation
@@ -400,6 +464,4 @@
     """
 
-    from Numeric import zeros, Float
-
     assert len(q) == 3,\
            'Vector of conserved quantities must have length 3'\
@@ -451,5 +513,4 @@
     from anuga.config import g, epsilon
     from math import sqrt
-    from Numeric import array
 
     #Align momentums with x-axis
@@ -539,5 +600,4 @@
 
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -614,5 +674,4 @@
     """
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -643,5 +702,4 @@
 
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -794,6 +852,4 @@
     """
 
-    from Numeric import zeros, Float
-
     N = domain.number_of_elements
     beta_h = domain.beta_h
@@ -865,6 +921,4 @@
     Wrapper for c-extension
     """
-
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -1046,5 +1100,4 @@
         self.normals = domain.normals
 
-        from Numeric import zeros, Float
         self.conserved_quantities = zeros(3, Float)
 
@@ -1204,6 +1257,4 @@
     """
 
-    from Numeric import zeros, Float, array, sum
-
     xmom = domain.quantities['xmomentum'].explicit_update
     ymom = domain.quantities['ymomentum'].explicit_update
@@ -1380,7 +1431,4 @@
     2: a scalar
     """
-
-    from Numeric import ones, Float, array
-
 
     if callable(f):
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3688)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3689)
@@ -342,5 +342,7 @@
     else
       alpha = 1.0;  //Flat bed
-
+      
+      
+    //alpha = 1.0;  //Always deep FIXME: This actually looks good now
 
     //printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
Index: /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 3688)
+++ /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 3689)
@@ -5,5 +5,5 @@
 
 from anuga.config import g, epsilon
-from Numeric import allclose, array, zeros, ones, Float, take
+from Numeric import allclose, alltrue, array, zeros, ones, Float, take
 from anuga.utilities.numerical_tools import mean
 
@@ -668,6 +668,8 @@
 
 
-    def test_catching_negative_heights(self):
-
+    def xtest_catching_negative_heights(self):
+
+        #OBSOLETE
+        
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -701,4 +703,244 @@
 
 
+
+    def test_get_wet_elements(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+        val0 = 2.+2.0/3
+        val1 = 4.+4.0/3
+        val2 = 8.+2.0/3
+        val3 = 2.+8.0/3
+
+        zl=zr=5
+        domain.set_quantity('elevation', zl*ones( (4,3) ))
+        domain.set_quantity('stage', [[val0, val0-1, val0-2],
+                                      [val1, val1+1, val1],
+                                      [val2, val2-2, val2],
+                                      [val3-0.5, val3, val3]])
+
+
+
+        #print domain.get_quantity('elevation').get_values(location='centroids')
+        #print domain.get_quantity('stage').get_values(location='centroids')        
+        domain.check_integrity()
+
+        indices = domain.get_wet_elements()
+        assert allclose(indices, [1,2])
+
+        indices = domain.get_wet_elements(indices=[0,1,3])
+        assert allclose(indices, [1])
+        
+
+
+    def test_get_maximum_inundation_1(self):
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
+        domain.set_quantity('stage', 3)
+
+        domain.check_integrity()
+
+        indices = domain.get_wet_elements()
+        assert allclose(indices, [0])
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
+
+        x, y = domain.get_maximum_inundation_location()
+        assert allclose([x, y], domain.get_centroid_coordinates()[0])
+
+
+    def test_get_maximum_inundation_2(self):
+        """test_get_maximum_inundation_2(self)
+
+        Test multiple wet cells with same elevation
+        """
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
+        domain.set_quantity('stage', 4.1)
+
+        domain.check_integrity()
+
+        indices = domain.get_wet_elements()
+        assert allclose(indices, [0,1,2])
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, 4)        
+
+        x, y = domain.get_maximum_inundation_location()
+        assert allclose([x, y], domain.get_centroid_coordinates()[1])        
+        
+
+    def test_get_maximum_inundation_3(self):
+        """test_get_maximum_inundation_3(self)
+
+        Test real runup example
+        """
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+
+        initial_runup_height = -0.4
+        final_runup_height = -0.3
+
+
+        #--------------------------------------------------------------
+        # Setup computational domain
+        #--------------------------------------------------------------
+        N = 5
+        points, vertices, boundary = rectangular_cross(N, N) 
+        domain = Domain(points, vertices, boundary)            
+
+        #--------------------------------------------------------------
+        # Setup initial conditions
+        #--------------------------------------------------------------
+        def topography(x,y):
+            return -x/2                             # linear bed slope
+            
+
+        domain.set_quantity('elevation', topography)       # Use function for elevation
+        domain.set_quantity('friction', 0.)                # Zero friction 
+        domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
+
+
+        #--------------------------------------------------------------
+        # Setup boundary conditions
+        #--------------------------------------------------------------
+        Br = Reflective_boundary(domain)              # Reflective wall
+        Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
+                                 0,
+                                 0])
+
+        # All reflective to begin with (still water) 
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+        #--------------------------------------------------------------
+        # Test initial inundation height
+        #--------------------------------------------------------------
+
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+            get_values(location='centroids', indices=indices)
+        assert alltrue(z<initial_runup_height)
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy 
+
+        x, y = domain.get_maximum_inundation_location()
+
+        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
+        assert allclose(q, qref)
+
+
+        wet_elements = domain.get_wet_elements()
+        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
+                                                                     indices=wet_elements)
+        assert alltrue(wet_elevations<initial_runup_height)
+        assert allclose(wet_elevations, z)        
+
+
+        #print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
+        #print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
+        #print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
+
+        
+        #--------------------------------------------------------------
+        # Let triangles adjust
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
+            pass
+
+
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+            get_values(location='centroids', indices=indices)
+
+        assert alltrue(z<initial_runup_height)
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
+        
+        x, y = domain.get_maximum_inundation_location()
+        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
+        assert allclose(q, qref)        
+
+
+        #--------------------------------------------------------------
+        # Update boundary to allow inflow
+        #--------------------------------------------------------------
+        domain.modify_boundary({'right': Bd})
+
+        
+        #--------------------------------------------------------------
+        # Evolve system through time
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+            #domain.write_time()
+            pass
+    
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+            get_values(location='centroids', indices=indices)
+
+        assert alltrue(z<final_runup_height)
+
+        q = domain.get_maximum_inundation_elevation()
+        assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy 
+
+        x, y = domain.get_maximum_inundation_location()
+        qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
+        assert allclose(q, qref)        
+
+
+        wet_elements = domain.get_wet_elements()
+        wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
+                                                                     indices=wet_elements)
+        assert alltrue(wet_elevations<final_runup_height)
+        assert allclose(wet_elevations, z)        
+        
 
 
@@ -2757,5 +2999,5 @@
 
         #Initial condition
-        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
+        domain.set_quantity('stage', expression = 'elevation + 0.05')
         domain.check_integrity()
 
@@ -3354,5 +3596,7 @@
          self.failUnless(domain.geo_reference.xllcorner  == 140.0,
                           "bad geo_referece")
-    #************
+
+
+         #************
 
     
@@ -3398,10 +3642,12 @@
          self.failUnless(domain.geo_reference.xllcorner  == 140.0,
                           "bad geo_referece")
-    #************
+         #************
          os.remove(fileName)
 
         #-------------------------------------------------------------
+        
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Shallow_Water,'test')
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_3')
     runner = unittest.TextTestRunner(verbosity=1)    
     runner.run(suite)
