Changeset 4820


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

Addressed ticket:217

Files:
5 edited
1 moved

Legend:

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

    r4809 r4820  
    198198        Boundary.__init__(self)
    199199
    200         #Get x,y vertex coordinates for all triangles
     200        # Get x,y vertex coordinates for all triangles
    201201        V = domain.vertex_coordinates
    202202
    203         #Compute midpoint coordinates for all boundary elements
    204         #Only a subset may be invoked when boundary is evaluated but
    205         #we don't know which ones at this stage since this object can
    206         #be attached to
    207         #any tagged boundary later on.
     203        # Compute midpoint coordinates for all boundary elements
     204        # Only a subset may be invoked when boundary is evaluated but
     205        # we don't know which ones at this stage since this object can
     206        # be attached to
     207        # any tagged boundary later on.
    208208
    209209        if verbose: print 'Find midpoint coordinates of entire boundary'
     
    216216       
    217217
    218         #Make ordering unique #FIXME: should this happen in domain.py?
     218        # Make ordering unique #FIXME: should this happen in domain.py?
    219219        boundary_keys.sort()
    220220
    221         #Record ordering #FIXME: should this also happen in domain.py?
     221        # Record ordering #FIXME: should this also happen in domain.py?
    222222        self.boundary_indices = {}
    223223        for i, (vol_id, edge_id) in enumerate(boundary_keys):
     
    228228            x2, y2 = V[base_index+2, :]
    229229           
    230             #Compute midpoints
     230            # Compute midpoints
    231231            if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
    232232            if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
    233233            if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
    234234
    235             #Convert to absolute UTM coordinates
     235            # Convert to absolute UTM coordinates
    236236            m[0] += xllcorner
    237237            m[1] += yllcorner
    238238           
    239             #Register point and index
     239            # Register point and index
    240240            self.midpoint_coordinates[i,:] = m
    241241
    242             #Register index of this boundary edge for use with evaluate
     242            # Register index of this boundary edge for use with evaluate
    243243            self.boundary_indices[(vol_id, edge_id)] = i
    244244
     
    252252        self.domain = domain
    253253
    254         #Test
     254        # Test
     255
     256        # Here we'll flag indices outside the mesh as a warning
     257        # as suggested by Joaquim Luis in sourceforge posting
     258        # November 2007
     259        # We won't make it an error as it is conceivable that
     260        # only part of mesh boundary is actually used with a given
     261        # file boundary sww file.
     262        if hasattr(self.F, 'indices_outside_mesh') and\
     263               len(self.F.indices_outside_mesh) > 0:
     264            msg = 'WARNING: File_boundary has points outside the mesh '
     265            msg += 'given in %s. ' %filename
     266            msg += 'See warning message issued by Interpolation_function '
     267            msg += 'for details (should appear above somewhere if '
     268            msg += 'verbose is True).\n'
     269            msg += 'This is perfectly OK as long as the points that are '
     270            msg += 'outside aren\'t used on the actual boundary segment.'
     271            if verbose is True:           
     272                print msg
     273            #raise Exception(msg)
     274
     275       
    255276        q = self.F(0, point_id=0)
    256277
     
    279300                x,y=self.midpoint_coordinates[i,:]
    280301                msg = 'NAN value found in file_boundary at '
    281                 msg += 'point id #%d: (%.2f, %.2f)' %(i, x, y)
    282                 #print msg
     302                msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
     303
     304                if hasattr(self.F, 'indices_outside_mesh') and\
     305                       len(self.F.indices_outside_mesh) > 0:
     306                    # Check if NAN point is due it being outside
     307                    # boundary defined in sww file.
     308
     309                    if i in self.F.indices_outside_mesh:
     310                        msg += 'This point refers to one outside the '
     311                        msg += 'mesh defined by the file %s.\n'\
     312                               %self.F.filename
     313                        msg += 'Make sure that the file covers '
     314                        msg += 'the boundary segment it is assigned to '
     315                        msg += 'in set_boundary.'
     316                    else:
     317                        msg += 'This point is inside the mesh defined '
     318                        msg += 'the file %s.\n' %self.F.filename
     319                        msg += 'Check this file for NANs.'
    283320                raise Exception, msg
    284321           
    285322            return res
    286323        else:
    287             #raise 'Boundary call without point_id not implemented'
    288             #FIXME: What should the semantics be?
     324            # raise 'Boundary call without point_id not implemented'
     325            # FIXME: What should the semantics be?
    289326            return self.F(t)
    290327
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r4787 r4820  
    7575
    7676
    77     #FIXME (OLE): Should check origin of domain against that of file
    78     #In fact, this is where origin should be converted to that of domain
    79     #Also, check that file covers domain fully.
    80 
    81     #Take into account:
    82     #- domain's georef
    83     #- sww file's georef
    84     #- interpolation points as absolute UTM coordinates
     77    # FIXME (OLE): Should check origin of domain against that of file
     78    # In fact, this is where origin should be converted to that of domain
     79    # Also, check that file covers domain fully.
     80
     81    # Take into account:
     82    # - domain's georef
     83    # - sww file's georef
     84    # - interpolation points as absolute UTM coordinates
    8585
    8686
     
    136136
    137137    f.starttime = starttime
     138    f.filename = filename
    138139   
    139140    if domain is not None:
     
    213214   
    214215   
    215     #FIXME: Check that model origin is the same as file's origin
    216     #(both in UTM coordinates)
    217     #If not - modify those from file to match domain
    218     #(origin should be passed in)
    219     #Take this code from e.g. dem2pts in data_manager.py
    220     #FIXME: Use geo_reference to read and write xllcorner...
     216    # FIXME: Check that model origin is the same as file's origin
     217    # (both in UTM coordinates)
     218    # If not - modify those from file to match domain
     219    # (origin should be passed in)
     220    # Take this code from e.g. dem2pts in data_manager.py
     221    # FIXME: Use geo_reference to read and write xllcorner...
    221222       
    222223
     
    226227    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
    227228
    228     #Open NetCDF file
     229    # Open NetCDF file
    229230    if verbose: print 'Reading', filename
    230231    fid = NetCDFFile(filename, 'r')
     
    234235   
    235236    if quantity_names is None or len(quantity_names) < 1:
    236         #If no quantities are specified get quantities from file
    237         #x, y, time are assumed as the independent variables so
    238         #they are excluded as potentiol quantities
     237        # If no quantities are specified get quantities from file
     238        # x, y, time are assumed as the independent variables so
     239        # they are excluded as potentiol quantities
    239240        quantity_names = []
    240241        for name in fid.variables.keys():
     
    253254
    254255
    255     #Now assert that requested quantitites (and the independent ones)
    256     #are present in file
     256    # Now assert that requested quantitites (and the independent ones)
     257    # are present in file
    257258    missing = []
    258259    for quantity in ['time'] + quantity_names:
     
    266267        raise Exception, msg
    267268
    268     #Decide whether this data has a spatial dimension
     269    # Decide whether this data has a spatial dimension
    269270    spatial = True
    270271    for quantity in ['x', 'y']:
     
    280281        raise msg       
    281282
    282     #Get first timestep
     283    # Get first timestep
    283284    try:
    284285        starttime = fid.starttime[0]
     
    287288        raise msg
    288289
    289     #Get variables
    290     #if verbose: print 'Get variables'   
     290    # Get variables
     291    # if verbose: print 'Get variables'   
    291292    time = fid.variables['time'][:]   
    292293
    293294    # Get time independent stuff
    294295    if spatial:
    295         #Get origin
     296        # Get origin
    296297        xllcorner = fid.xllcorner[0]
    297298        yllcorner = fid.yllcorner[0]
     
    307308
    308309        if interpolation_points is not None:
    309             #Adjust for georef
     310            # Adjust for georef
    310311            interpolation_points[:,0] -= xllcorner
    311312            interpolation_points[:,1] -= yllcorner       
     
    316317    if domain_starttime is not None:
    317318
    318         #If domain_startime is *later* than starttime,
    319         #move time back - relative to domain's time
     319        # If domain_startime is *later* than starttime,
     320        # move time back - relative to domain's time
    320321        if domain_starttime > starttime:
    321322            time = time - domain_starttime + starttime
    322323
    323         #FIXME Use method in geo to reconcile
    324         #if spatial:
    325         #assert domain.geo_reference.xllcorner == xllcorner
    326         #assert domain.geo_reference.yllcorner == yllcorner
    327         #assert domain.geo_reference.zone == zone       
     324        # FIXME Use method in geo to reconcile
     325        # if spatial:
     326        # assert domain.geo_reference.xllcorner == xllcorner
     327        # assert domain.geo_reference.yllcorner == yllcorner
     328        # assert domain.geo_reference.zone == zone       
    328329       
    329330    if verbose:
     
    345346   
    346347
    347     #from least_squares import Interpolation_function
    348348    from anuga.fit_interpolate.interpolate import Interpolation_function
    349349
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r4779 r4820  
    521521
    522522
     523            # Check that all interpolation points fall within
     524            # mesh boundary as defined by triangles and vertex_coordinates.
     525            from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     526            from anuga.utilities.polygon import outside_polygon           
     527           
     528            mesh = Mesh(vertex_coordinates, triangles)
     529           
     530            indices = outside_polygon(interpolation_points,
     531                                      mesh.get_boundary_polygon())
     532
     533
     534            # Record result
     535            self.indices_outside_mesh = indices
     536
     537            # Report
     538            if len(indices) > 0:
     539                msg = 'Interpolation points in Interpolation function fall '
     540                msg += 'outside specified mesh. '
     541                msg += 'Offending points:\n'
     542                for i in indices:
     543                    msg += '%d: %s\n' %(i, interpolation_points[i])
     544
     545                # Joaquim Luis suggested this as an Exception, so
     546                # that the user can now what the problem is rather than
     547                # looking for NaN's. However, NANs are handy as they can
     548                # be ignored leaving good points for continued processing.
     549                if verbose:
     550                    print msg
     551                #raise Exception(msg)
     552
     553
     554           
     555
    523556            m = len(self.interpolation_points)
    524557            p = len(self.time)
     
    626659        oldindex = self.index #Time index
    627660
    628         #Find current time slot
     661        # Find current time slot
    629662        while t > self.time[self.index]: self.index += 1
    630663        while t < self.time[self.index]: self.index -= 1
    631664
    632665        if t == self.time[self.index]:
    633             #Protect against case where t == T[-1] (last time)
    634             # - also works in general when t == T[i]
     666            # Protect against case where t == T[-1] (last time)
     667            #  - also works in general when t == T[i]
    635668            ratio = 0
    636669        else:
    637             #t is now between index and index+1
     670            # t is now between index and index+1
    638671            ratio = (t - self.time[self.index])/\
    639672                    (self.time[self.index+1] - self.time[self.index])
    640673
    641         #Compute interpolated values
     674        # Compute interpolated values
    642675        q = zeros(len(self.quantity_names), Float)
    643         #print "self.precomputed_values", self.precomputed_values
     676        # print "self.precomputed_values", self.precomputed_values
    644677        for i, name in enumerate(self.quantity_names):
    645678            Q = self.precomputed_values[name]
    646679
    647680            if self.spatial is False:
    648                 #If there is no spatial info               
     681                # If there is no spatial info               
    649682                assert len(Q.shape) == 1
    650683
     
    654687            else:
    655688                if x is not None and y is not None:
    656                     #Interpolate to x, y
     689                    # Interpolate to x, y
    657690                   
    658691                    raise 'x,y interpolation not yet implemented'
    659692                else:
    660                     #Use precomputed point
     693                    # Use precomputed point
    661694                    Q0 = Q[self.index, point_id]
    662695                    if ratio > 0:
    663696                        Q1 = Q[self.index+1, point_id]
    664697
    665             #Linear temporal interpolation   
     698            # Linear temporal interpolation   
    666699            if ratio > 0:
    667700                if Q0 == NAN and Q1 == NAN:
     
    673706
    674707
    675         #Return vector of interpolated values
    676         #if len(q) == 1:
    677         #    return q[0]
    678         #else:
    679         #    return q
    680 
    681 
    682         #Return vector of interpolated values
    683         #FIXME:
     708        # Return vector of interpolated values
     709        # if len(q) == 1:
     710        #     return q[0]
     711        # else:
     712        #     return q
     713
     714
     715        # Return vector of interpolated values
     716        # FIXME:
    684717        if self.spatial is True:
    685718            return q
    686719        else:
    687             #Replicate q according to x and y
    688             #This is e.g used for Wind_stress
     720            # Replicate q according to x and y
     721            # This is e.g used for Wind_stress
    689722            if x is None or y is None:
    690723                return q
     
    696729                else:
    697730                    from Numeric import ones, Float
    698                     #x is a vector - Create one constant column for each value
     731                    # x is a vector - Create one constant column for each value
    699732                    N = len(x)
    700733                    assert len(y) == N, 'x and y must have same length'
  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r4779 r4820  
    12991299        # interpolations in both time and space
    13001300   
    1301         #Three timesteps
     1301        # Three timesteps
    13021302        time = [1.0, 5.0, 6.0]   
    13031303
    1304         #Setup mesh used to represent fitted function
     1304        # Setup mesh used to represent fitted function
    13051305        a = [0.0, 0.0]
    13061306        b = [0.0, 2.0]
     
    13151315
    13161316
    1317         #New datapoints where interpolated values are sought
     1317        # New datapoints where interpolated values are sought
    13181318        interpolation_points = [[ 0.0, 0.0],
    13191319                                [ 0.5, 0.5],
     
    13231323                                [ 545354534, 4354354353]] # outside the mesh
    13241324
    1325         #One quantity
     1325        # One quantity
    13261326        Q = zeros( (3,6), Float )
    13271327
    1328         #Linear in time and space
     1328        # Linear in time and space
    13291329        for i, t in enumerate(time):
    13301330            Q[i, :] = t*linear_function(points)
    13311331
    1332         #Check interpolation of one quantity using interpolaton points)
     1332        # Check interpolation of one quantity using interpolaton points)
     1333
    13331334        I = Interpolation_function(time, Q,
    13341335                                   vertex_coordinates = points,
     
    13361337                                   interpolation_points = interpolation_points,
    13371338                                   verbose = False)
    1338 
     1339       
    13391340       
    13401341        assert alltrue(I.precomputed_values['Attribute'][:,4] != NAN)
  • anuga_validation/automated_validation_tests/pixel_registration/README.txt

    r4810 r4820  
    88the boundary condition.
    99
    10 After Joaquim's fix (November 2007), this passed.
     10--
     11It turned out that this test didn't reveal the registration problem after all
     12(see sourceforge postings Nov 2007), but rather that ANUGA didn't
     13give a proper error if the domain exceeds the boundary defined in Field_boundary.
    1114
    12 The file layer01.png shows the polygons constituting this example.
     15We don't want to raise an Exception, but the code now has appropriate warnings.
     16This directory is therefore not really part of the validation suite after all.
     17 
Note: See TracChangeset for help on using the changeset viewer.