Changeset 4820

Timestamp:

Nov 16, 2007, 2:23:00 PM (17 years ago)

Author:

ole

Message:

Addressed ticket:217

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (5 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/util.py (modified) (12 diffs)
• anuga_core/source/anuga/fit_interpolate/interpolate.py (modified) (5 diffs)
• anuga_core/source/anuga/fit_interpolate/test_interpolate.py (modified) (4 diffs)
• anuga_validation/automated_validation_tests/pixel_registration/README.txt (modified) (1 diff)
• anuga_validation/automated_validation_tests/pixel_registration/no_validate_pixel_registration.py (moved) (moved from anuga_validation/automated_validation_tests/pixel_registration/validate_pixel_registration.py)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -198 +198 @@
         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'
@@ -216 +216 @@
         
 
-        #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):
@@ -228 +228 @@
             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
 
@@ -252 +252 @@
         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)
 
@@ -279 +300 @@
                 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

@@ -75 +75 @@
 
 
-    #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
 
 
@@ -136 +136 @@
 
     f.starttime = starttime
+    f.filename = filename
     
     if domain is not None:
@@ -213 +214 @@
     
     
-    #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...
         
 
@@ -226 +227 @@
     from Numeric import array, zeros, Float, alltrue, concatenate, reshape
 
-    #Open NetCDF file
+    # Open NetCDF file
     if verbose: print 'Reading', filename
     fid = NetCDFFile(filename, 'r')
@@ -234 +235 @@
     
     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():
@@ -253 +254 @@
 
 
-    #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:
@@ -266 +267 @@
         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']:
@@ -280 +281 @@
         raise msg        
 
-    #Get first timestep
+    # Get first timestep
     try:
         starttime = fid.starttime[0]
@@ -287 +288 @@
         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]
@@ -307 +308 @@
 
         if interpolation_points is not None:
-            #Adjust for georef
+            # Adjust for georef
             interpolation_points[:,0] -= xllcorner
             interpolation_points[:,1] -= yllcorner        
@@ -316 +317 @@
     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:
@@ -345 +346 @@
     
 
-    #from least_squares import Interpolation_function
     from anuga.fit_interpolate.interpolate import Interpolation_function
 

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -521 +521 @@
 
 
+            # 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)
@@ -626 +659 @@
         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
 
@@ -654 +687 @@
             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:
@@ -673 +706 @@
 
 
-        #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
@@ -696 +729 @@
                 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

@@ -1299 +1299 @@
         # 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]
@@ -1315 +1315 @@
 
 
-        #New datapoints where interpolated values are sought
+        # New datapoints where interpolated values are sought
         interpolation_points = [[ 0.0, 0.0],
                                 [ 0.5, 0.5],
@@ -1323 +1323 @@
                                 [ 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,
@@ -1336 +1337 @@
                                    interpolation_points = interpolation_points,
                                    verbose = False)
-
+        
         
         assert alltrue(I.precomputed_values['Attribute'][:,4] != NAN)

anuga_validation/automated_validation_tests/pixel_registration/README.txt

@@ -8 +8 @@
 the boundary condition.
 
-After Joaquim's fix (November 2007), this passed.
+--
+It turned out that this test didn't reveal the registration problem after all
+(see sourceforge postings Nov 2007), but rather that ANUGA didn't 
+give a proper error if the domain exceeds the boundary defined in Field_boundary.
 
-The file layer01.png shows the polygons constituting this example. 
+We don't want to raise an Exception, but the code now has appropriate warnings.
+This directory is therefore not really part of the validation suite after all.
+