# Changeset 3689

Ignore:
Timestamp:
Oct 4, 2006, 11:00:56 AM (17 years ago)
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:
18 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

 r3678 See methods inside the quantity object for more options FIXME: clean input args """ 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

 r3684 def __repr__(self): return 'Mesh: %d vertex coordinates, %d triangles'\ return 'Mesh: %d vertices, %d triangles'\ %(self.coordinates.shape[0], len(self)) 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 (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 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

 r3688 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

 r3678 """ from Numeric import array, zeros, Float, concatenate, NewAxis, argmax, allclose from anuga.utilities.numerical_tools import ensure_numeric class Quantity: from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh from Numeric import array, zeros, Float msg = 'First argument in Quantity.__init__ ' 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 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. """ 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']: # 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: 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

 r3560 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): assert unique_vertices == [0,2,4,5,6,7] #------------------------------------------------------------- if __name__ == "__main__":
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

 r3678 [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): #      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

 r3678 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

 r3567 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

 r3678 # 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

 r3563 # Import necessary modules #------------------------------------------------------------------------------ from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross from anuga.shallow_water import Domain 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 # 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 # 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 # 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 # 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

 r3630 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 # 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 interior_regions=[[interior_polygon,dx*dy/32]]) domain = Domain(meshname, use_cache=True, verbose = True) domain.set_minimum_storable_height(0.01) 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, # 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

 r3566 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. 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

 r3686 a mesh origin, since geospatial has its own mesh origin. """ # FIXME (Ole): Need an input check # Initialise variabels 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): 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) 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 \ verbose=verbose) else: #print 'BLOCKING' #Handle blocking self._A_can_be_reused = False start=0 start = 0 # creating a dummy array to concatenate to. 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) return z def interpolate_block(self, f, point_coordinates = None, verbose=False): """ 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) 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 """ #print 'Building interpolation matrix' #Convert point_coordinates to Numeric arrays, in case it was a list. #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' 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