Changeset 3689

Timestamp:

Oct 4, 2006, 11:00:56 AM (18 years ago)

Author:

ole

Message:

Added functionality for getting arbitrary interpolated values in Quantity as well as calculating inundation height and location. This work was done at SUT during the last week of September 2006.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (2 diffs)
• abstract_2d_finite_volumes/general_mesh.py (modified) (4 diffs)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/quantity.py (modified) (7 diffs)
• abstract_2d_finite_volumes/test_general_mesh.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (2 diffs)
• config.py (modified) (1 diff)
• damage_modelling/test_inundation_damage.py (modified) (1 diff)
• examples/netherlands.py (modified) (1 diff)
• examples/runup.py (modified) (6 diffs)
• examples/runup_convergence.py (modified) (5 diffs)
• fit_interpolate/general_fit_interpolate.py (modified) (2 diffs)
• fit_interpolate/interpolate.py (modified) (12 diffs)
• fit_interpolate/test_interpolate.py (modified) (4 diffs)
• pmesh/mesh_interface.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (12 diffs)
• shallow_water/shallow_water_ext.c (modified) (1 diff)
• shallow_water/test_shallow_water_domain.py (modified) (6 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -318 +318 @@
 
         See methods inside the quantity object for more options
+
+        FIXME: clean input args
         """
 
@@ -498 +500 @@
         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):

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -178 +178 @@
 
     def __repr__(self):
-        return 'Mesh: %d vertex coordinates, %d triangles'\
+        return 'Mesh: %d vertices, %d triangles'\
                %(self.coordinates.shape[0], len(self))
 
@@ -197 +197 @@
 
 
-    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 +212 @@
         (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 +236 @@
                 
         if obj is True:
-
             N = V.shape[0]
             return reshape(V, (3*N, 2))

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -175 +175 @@
 
     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))
 
 

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -15 +15 @@
 """
 
+from Numeric import array, zeros, Float, concatenate, NewAxis, argmax, allclose
+from anuga.utilities.numerical_tools import ensure_numeric
 
 class Quantity:
@@ -21 +23 @@
 
         from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
-        from Numeric import array, zeros, Float
 
         msg = 'First argument in Quantity.__init__ '
@@ -619 +620 @@
 
 
-        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 +752 @@
 
 
-    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 +892 @@
         """
         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 +930 @@
                 # 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 +1063 @@
             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()

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -18 +18 @@
     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 +83 @@
         assert unique_vertices == [0,2,4,5,6,7]
 
+
+        
+
 #-------------------------------------------------------------
 if __name__ == "__main__":

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -110 +110 @@
                                                [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 +1474 @@
         #      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):
         """

anuga_core/source/anuga/config.py

@@ -117 +117 @@
 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

anuga_core/source/anuga/damage_modelling/test_inundation_damage.py

@@ -142 +142 @@
     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)

anuga_core/source/anuga/examples/netherlands.py

@@ -166 +166 @@
 #
 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

anuga_core/source/anuga/examples/runup.py

@@ -9 +9 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
-
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 from anuga.shallow_water import Domain
@@ -15 +14 @@
 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 +20 @@
 # 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 +30 @@
 # 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 +42 @@
 # 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 +58 @@
 # 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()
     

anuga_core/source/anuga/examples/runup_convergence.py

@@ -34 +34 @@
 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 +52 @@
 
 # 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 +74 @@
                          interior_regions=[[interior_polygon,dx*dy/32]])
 domain = Domain(meshname, use_cache=True, verbose = True)
-
+domain.set_minimum_storable_height(0.01)
 
 
@@ -134 +134 @@
 
 
+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 +168 @@
     # 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 )
+
+
 
 #-----------------------------------------------------------------------------

anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

@@ -77 +77 @@
         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 +83 @@
         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)

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -76 +76 @@
               a mesh origin, since geospatial has its own mesh origin.
         """
+
+        # FIXME (Ole): Need an input check
         
         # Initialise variabels
@@ -82 +84 @@
         
         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 +116 @@
           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 +143 @@
                 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 +152 @@
                                            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 +165 @@
                     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 +180 @@
         return z
 
+
     def interpolate_block(self, f, point_coordinates = None, verbose=False):
         """
@@ -184 +194 @@
                 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 +216 @@
 
     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 +234 @@
         """
 
-
+        #print 'Building interpolation matrix'
         
         #Convert point_coordinates to Numeric arrays, in case it was a list.
@@ -486 +496 @@
             #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 +511 @@
                     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 

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -105 +105 @@
         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 +826 @@
 
         
-    def test_interpolate_reuse(self):
+    def test_interpolate_reuse_if_None(self):
         a = [-1.0, 0.0]
         b = [3.0, 4.0]
@@ -827 +883 @@
         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 +1572 @@
 
     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)

anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -28 +28 @@
                              mesh_geo_reference=None,
                              minimum_triangle_angle=28.0,
+                             use_cache=False,
                              verbose=True):
     """Create mesh from bounding polygons, and resolutions.

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -70 +70 @@
 #$LastChangedRevision$
 #$LastChangedBy$
+
+from Numeric import zeros, ones, Float, array, sum, size
+from Numeric import compress, arange
 
 
@@ -217 +220 @@
 
 
+
+    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 +464 @@
     """
 
-    from Numeric import zeros, Float
-
     assert len(q) == 3,\
            'Vector of conserved quantities must have length 3'\
@@ -451 +513 @@
     from anuga.config import g, epsilon
     from math import sqrt
-    from Numeric import array
 
     #Align momentums with x-axis
@@ -539 +600 @@
 
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -614 +674 @@
     """
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -643 +702 @@
 
     import sys
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -794 +852 @@
     """
 
-    from Numeric import zeros, Float
-
     N = domain.number_of_elements
     beta_h = domain.beta_h
@@ -865 +921 @@
     Wrapper for c-extension
     """
-
-    from Numeric import zeros, Float
 
     N = domain.number_of_elements
@@ -1046 +1100 @@
         self.normals = domain.normals
 
-        from Numeric import zeros, Float
         self.conserved_quantities = zeros(3, Float)
 
@@ -1204 +1257 @@
     """
 
-    from Numeric import zeros, Float, array, sum
-
     xmom = domain.quantities['xmomentum'].explicit_update
     ymom = domain.quantities['ymomentum'].explicit_update
@@ -1380 +1431 @@
     2: a scalar
     """
-
-    from Numeric import ones, Float, array
-
 
     if callable(f):

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -342 +342 @@
     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);

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -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 +668 @@
 
 
-    def test_catching_negative_heights(self):
-
+    def xtest_catching_negative_heights(self):
+
+        #OBSOLETE
+        
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -701 +703 @@
 
 
+
+    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 +2999 @@
 
         #Initial condition
-        domain.set_quantity('stage', Constant_height(x_slope, 0.05))
+        domain.set_quantity('stage', expression = 'elevation + 0.05')
         domain.check_integrity()
 
@@ -3354 +3596 @@
          self.failUnless(domain.geo_reference.xllcorner  == 140.0,
                           "bad geo_referece")
-    #************
+
+
+         #************
 
     
@@ -3398 +3642 @@
          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)