# Changeset 4820

Ignore:
Timestamp:
Nov 16, 2007, 2:23:00 PM (15 years ago)
Message:

Files:
5 edited
1 moved

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

 r4809 Boundary.__init__(self) #Get x,y vertex coordinates for all triangles # Get x,y vertex coordinates for all triangles V = domain.vertex_coordinates #Compute midpoint coordinates for all boundary elements #Only a subset may be invoked when boundary is evaluated but #we don't know which ones at this stage since this object can #be attached to #any tagged boundary later on. # Compute midpoint coordinates for all boundary elements # Only a subset may be invoked when boundary is evaluated but # we don't know which ones at this stage since this object can # be attached to # any tagged boundary later on. if verbose: print 'Find midpoint coordinates of entire boundary' #Make ordering unique #FIXME: should this happen in domain.py? # Make ordering unique #FIXME: should this happen in domain.py? boundary_keys.sort() #Record ordering #FIXME: should this also happen in domain.py? # Record ordering #FIXME: should this also happen in domain.py? self.boundary_indices = {} for i, (vol_id, edge_id) in enumerate(boundary_keys): x2, y2 = V[base_index+2, :] #Compute midpoints # Compute midpoints if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2]) if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2]) if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2]) #Convert to absolute UTM coordinates # Convert to absolute UTM coordinates m[0] += xllcorner m[1] += yllcorner #Register point and index # Register point and index self.midpoint_coordinates[i,:] = m #Register index of this boundary edge for use with evaluate # Register index of this boundary edge for use with evaluate self.boundary_indices[(vol_id, edge_id)] = i self.domain = domain #Test # Test # Here we'll flag indices outside the mesh as a warning # as suggested by Joaquim Luis in sourceforge posting # November 2007 # We won't make it an error as it is conceivable that # only part of mesh boundary is actually used with a given # file boundary sww file. if hasattr(self.F, 'indices_outside_mesh') and\ len(self.F.indices_outside_mesh) > 0: msg = 'WARNING: File_boundary has points outside the mesh ' msg += 'given in %s. ' %filename msg += 'See warning message issued by Interpolation_function ' msg += 'for details (should appear above somewhere if ' msg += 'verbose is True).\n' msg += 'This is perfectly OK as long as the points that are ' msg += 'outside aren\'t used on the actual boundary segment.' if verbose is True: print msg #raise Exception(msg) q = self.F(0, point_id=0) x,y=self.midpoint_coordinates[i,:] msg = 'NAN value found in file_boundary at ' msg += 'point id #%d: (%.2f, %.2f)' %(i, x, y) #print msg msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y) if hasattr(self.F, 'indices_outside_mesh') and\ len(self.F.indices_outside_mesh) > 0: # Check if NAN point is due it being outside # boundary defined in sww file. if i in self.F.indices_outside_mesh: msg += 'This point refers to one outside the ' msg += 'mesh defined by the file %s.\n'\ %self.F.filename msg += 'Make sure that the file covers ' msg += 'the boundary segment it is assigned to ' msg += 'in set_boundary.' else: msg += 'This point is inside the mesh defined ' msg += 'the file %s.\n' %self.F.filename msg += 'Check this file for NANs.' raise Exception, msg return res else: #raise 'Boundary call without point_id not implemented' #FIXME: What should the semantics be? # raise 'Boundary call without point_id not implemented' # FIXME: What should the semantics be? return self.F(t)
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

 r4787 #FIXME (OLE): Should check origin of domain against that of file #In fact, this is where origin should be converted to that of domain #Also, check that file covers domain fully. #Take into account: #- domain's georef #- sww file's georef #- interpolation points as absolute UTM coordinates # FIXME (OLE): Should check origin of domain against that of file # In fact, this is where origin should be converted to that of domain # Also, check that file covers domain fully. # Take into account: # - domain's georef # - sww file's georef # - interpolation points as absolute UTM coordinates f.starttime = starttime f.filename = filename if domain is not None: #FIXME: Check that model origin is the same as file's origin #(both in UTM coordinates) #If not - modify those from file to match domain #(origin should be passed in) #Take this code from e.g. dem2pts in data_manager.py #FIXME: Use geo_reference to read and write xllcorner... # FIXME: Check that model origin is the same as file's origin # (both in UTM coordinates) # If not - modify those from file to match domain # (origin should be passed in) # Take this code from e.g. dem2pts in data_manager.py # FIXME: Use geo_reference to read and write xllcorner... from Numeric import array, zeros, Float, alltrue, concatenate, reshape #Open NetCDF file # Open NetCDF file if verbose: print 'Reading', filename fid = NetCDFFile(filename, 'r') if quantity_names is None or len(quantity_names) < 1: #If no quantities are specified get quantities from file #x, y, time are assumed as the independent variables so #they are excluded as potentiol quantities # If no quantities are specified get quantities from file # x, y, time are assumed as the independent variables so # they are excluded as potentiol quantities quantity_names = [] for name in fid.variables.keys(): #Now assert that requested quantitites (and the independent ones) #are present in file # Now assert that requested quantitites (and the independent ones) # are present in file missing = [] for quantity in ['time'] + quantity_names: raise Exception, msg #Decide whether this data has a spatial dimension # Decide whether this data has a spatial dimension spatial = True for quantity in ['x', 'y']: raise msg #Get first timestep # Get first timestep try: starttime = fid.starttime[0] raise msg #Get variables #if verbose: print 'Get variables' # Get variables # if verbose: print 'Get variables' time = fid.variables['time'][:] # Get time independent stuff if spatial: #Get origin # Get origin xllcorner = fid.xllcorner[0] yllcorner = fid.yllcorner[0] if interpolation_points is not None: #Adjust for georef # Adjust for georef interpolation_points[:,0] -= xllcorner interpolation_points[:,1] -= yllcorner if domain_starttime is not None: #If domain_startime is *later* than starttime, #move time back - relative to domain's time # If domain_startime is *later* than starttime, # move time back - relative to domain's time if domain_starttime > starttime: time = time - domain_starttime + starttime #FIXME Use method in geo to reconcile #if spatial: #assert domain.geo_reference.xllcorner == xllcorner #assert domain.geo_reference.yllcorner == yllcorner #assert domain.geo_reference.zone == zone # FIXME Use method in geo to reconcile # if spatial: # assert domain.geo_reference.xllcorner == xllcorner # assert domain.geo_reference.yllcorner == yllcorner # assert domain.geo_reference.zone == zone if verbose: #from least_squares import Interpolation_function from anuga.fit_interpolate.interpolate import Interpolation_function
• ## anuga_core/source/anuga/fit_interpolate/interpolate.py

 r4779 # Check that all interpolation points fall within # mesh boundary as defined by triangles and vertex_coordinates. from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh from anuga.utilities.polygon import outside_polygon mesh = Mesh(vertex_coordinates, triangles) indices = outside_polygon(interpolation_points, mesh.get_boundary_polygon()) # Record result self.indices_outside_mesh = indices # Report if len(indices) > 0: msg = 'Interpolation points in Interpolation function fall ' msg += 'outside specified mesh. ' msg += 'Offending points:\n' for i in indices: msg += '%d: %s\n' %(i, interpolation_points[i]) # Joaquim Luis suggested this as an Exception, so # that the user can now what the problem is rather than # looking for NaN's. However, NANs are handy as they can # be ignored leaving good points for continued processing. if verbose: print msg #raise Exception(msg) m = len(self.interpolation_points) p = len(self.time) oldindex = self.index #Time index #Find current time slot # Find current time slot while t > self.time[self.index]: self.index += 1 while t < self.time[self.index]: self.index -= 1 if t == self.time[self.index]: #Protect against case where t == T[-1] (last time) # - also works in general when t == T[i] # Protect against case where t == T[-1] (last time) #  - also works in general when t == T[i] ratio = 0 else: #t is now between index and index+1 # t is now between index and index+1 ratio = (t - self.time[self.index])/\ (self.time[self.index+1] - self.time[self.index]) #Compute interpolated values # Compute interpolated values q = zeros(len(self.quantity_names), Float) #print "self.precomputed_values", self.precomputed_values # print "self.precomputed_values", self.precomputed_values for i, name in enumerate(self.quantity_names): Q = self.precomputed_values[name] if self.spatial is False: #If there is no spatial info # If there is no spatial info assert len(Q.shape) == 1 else: if x is not None and y is not None: #Interpolate to x, y # Interpolate to x, y raise 'x,y interpolation not yet implemented' else: #Use precomputed point # Use precomputed point Q0 = Q[self.index, point_id] if ratio > 0: Q1 = Q[self.index+1, point_id] #Linear temporal interpolation # Linear temporal interpolation if ratio > 0: if Q0 == NAN and Q1 == NAN: #Return vector of interpolated values #if len(q) == 1: #    return q[0] #else: #    return q #Return vector of interpolated values #FIXME: # Return vector of interpolated values # if len(q) == 1: #     return q[0] # else: #     return q # Return vector of interpolated values # FIXME: if self.spatial is True: return q else: #Replicate q according to x and y #This is e.g used for Wind_stress # Replicate q according to x and y # This is e.g used for Wind_stress if x is None or y is None: return q else: from Numeric import ones, Float #x is a vector - Create one constant column for each value # x is a vector - Create one constant column for each value N = len(x) assert len(y) == N, 'x and y must have same length'
• ## anuga_core/source/anuga/fit_interpolate/test_interpolate.py

 r4779 # interpolations in both time and space #Three timesteps # Three timesteps time = [1.0, 5.0, 6.0] #Setup mesh used to represent fitted function # Setup mesh used to represent fitted function a = [0.0, 0.0] b = [0.0, 2.0] #New datapoints where interpolated values are sought # New datapoints where interpolated values are sought interpolation_points = [[ 0.0, 0.0], [ 0.5, 0.5], [ 545354534, 4354354353]] # outside the mesh #One quantity # One quantity Q = zeros( (3,6), Float ) #Linear in time and space # Linear in time and space for i, t in enumerate(time): Q[i, :] = t*linear_function(points) #Check interpolation of one quantity using interpolaton points) # Check interpolation of one quantity using interpolaton points) I = Interpolation_function(time, Q, vertex_coordinates = points, interpolation_points = interpolation_points, verbose = False) assert alltrue(I.precomputed_values['Attribute'][:,4] != NAN)