Changeset 7685

Timestamp:

Apr 14, 2010, 6:53:12 PM (13 years ago)

Author:

hudson

Message:

Ticket 344 completed and unit tests pass.
Stripped out a heap of tabs and replaced them with 4 spaces.
Removed some redundant imports and deprecated functions.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/gauge.py (modified) (6 diffs)
• abstract_2d_finite_volumes/test_gauge.py (modified) (6 diffs)
• abstract_2d_finite_volumes/util.py (modified) (6 diffs)
• fit_interpolate/general_fit_interpolate.py (modified) (1 diff)
• fit_interpolate/interpolate.py (modified) (14 diffs)
• fit_interpolate/test_interpolate.py (modified) (3 diffs)
• shallow_water/data_manager.py (modified) (2 diffs)
• shallow_water/data_manager_joaquims_patch.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

@@ -92 +92 @@
     import string
     from anuga.shallow_water.data_manager import get_all_swwfiles
-    from anuga.abstract_2d_finite_volumes.util import file_function     
+    from anuga.abstract_2d_finite_volumes.util import file_function    
 
     assert type(gauge_file) == type(''), 'Gauge filename must be a string'
@@ -246 +246 @@
                         if quantity == 'ycentroid':
                             points_list.append(callable_sww.centroids[point_i][1])
-                                                        
+                            
                 points_writer[point_i].writerow(points_list)
 
@@ -281 +281 @@
                    use_cache=False,
                    verbose=False,
-                                   output_centroids=False):
+                   output_centroids=False):
     """ Read sww file and plot the time series for the
     prescribed quantities at defined gauge locations and
@@ -399 +399 @@
                         use_cache,
                         verbose,
-                                                output_centroids = output_centroids)
+                        output_centroids = output_centroids)
     return k
 
@@ -435 +435 @@
                     use_cache = False,
                     verbose = False,
-                                        output_centroids = False):   
+                    output_centroids = False):   
         
     # FIXME(Ole): Shouldn't print statements here be governed by verbose?
@@ -499 +499 @@
                           verbose = verbose,
                           use_cache = use_cache,
-                                                  output_centroids = output_centroids)
+                          output_centroids = output_centroids)
 
         # determine which gauges are contained in sww file

anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

@@ -389 +389 @@
         os.remove(point1_filename)
         os.remove(point2_filename)
-                
+        
 
     def test_sww2csv_centroid(self):
@@ -399 +399 @@
         
         Test the ability to get a timeseries at the centroid of a triangle, rather
-                than the given gauge point.
+        than the given gauge point.
         """
         
@@ -453 +453 @@
         points_file = tempfile.mktemp(".csv")
         file_id = open(points_file,"w")
-# These values are where the centroids should be                
+# These values are where the centroids should be        
 #        file_id.write("name, easting, northing, elevation \n\
 #point1, 2.0, 2.0, 3.0\n\
@@ -501 +501 @@
         point1_handle.close()
         point2_handle.close()
-        os.remove(mesh_file)            
+        os.remove(mesh_file)        
         os.remove(sww.filename)
         os.remove(points_file)
@@ -516 +516 @@
         
         Test the ability to get a timeseries at the centroid of a triangle, rather
-                than the given gauge point, then output the results.
+        than the given gauge point, then output the results.
         """
         
@@ -596 +596 @@
         # clean up
         point1_handle.close()
-        os.remove(mesh_file)            
+        os.remove(mesh_file)        
         os.remove(sww.filename)
         os.remove(points_file)
         os.remove(point1_filename)
-                
-                
+        
+        
 
 #-------------------------------------------------------------

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -158 +158 @@
     """Temporary Interface to new location"""
 
-    import anuga.utilities.numerical_tools as NT        
+    import anuga.utilities.numerical_tools as NT    
     
     msg = 'angle has moved from util.py.  '
@@ -187 +187 @@
     warn(msg, DeprecationWarning) 
 
-    return utilities.polygon.point_on_line(*args, **kwargs)     
+    return utilities.polygon.point_on_line(*args, **kwargs)    
 
     
@@ -211 +211 @@
     return utilities.polygon.outside_polygon(*args, **kwargs)    
 
-
-# @note TEMP
-def separate_points_by_polygon(*args, **kwargs):
-    """Temporary Interface to new location"""
-
-    log.critical('separate_points_by_polygon has moved from util.py.')
-    log.critical('Please use "from anuga.utilities.polygon import '
-                 'separate_points_by_polygon"')
-
-    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)    
-
-
+    
 # @note TEMP
 def read_polygon(*args, **kwargs):
@@ -947 +936 @@
         # verts=num.take(verts,X,axis=0)  to Remove the loners from verts
         # but I think it would use more memory
-        new_i = lone_start      # point at first loner - 'shuffle down' target
+        new_i = lone_start    # point at first loner - 'shuffle down' target
         for i in range(lone_start, N):
-            if loners[i] >= N:  # [i] is a loner, leave alone
+            if loners[i] >= N:    # [i] is a loner, leave alone
                 pass
-            else:               # a non-loner, move down
+            else:        # a non-loner, move down
                 loners[i] = new_i
                 verts[new_i] = verts[i]
@@ -1575 +1564 @@
     """
     from gauge import sww2csv_gauges as sww2csv
-        
+    
     return sww2csv(sww_file, gauge_file, out_name, quantities, verbose, use_cache)
-        
+    
 def sww2timeseries(swwfiles,
                    gauge_filename,
@@ -1610 +1599 @@
                    use_cache,
                    verbose,
-                   output_centroids)                               
-        
+                   output_centroids)                   
+    
 ##
 # @brief Get a wave height at a certain depth given wave height at another depth.

anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

@@ -27 +27 @@
 from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError
 from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.utilities.polygon import in_and_outside_polygon
 from anuga.utilities.quad import build_quadtree
 

anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -155 +155 @@
                                  use_cache=use_cache,
                                  verbose=verbose,
-                                                                 output_centroids=output_centroids)
+                                 output_centroids=output_centroids)
 
     return result
@@ -230 +230 @@
                     start_blocking_len=500000,
                     verbose=False,
-                                        output_centroids=False):
+                    output_centroids=False):
         """Interpolate mesh data f to determine values, z, at points.
 
@@ -324 +324 @@
     # @param use_cache True if caching should be used to accelerate the calculations
     # @param verbose True if this function is verbose.
-    # @return ??
+    # @return interpolated f
     def interpolate_block(self, f, point_coordinates,
                           use_cache=False, verbose=False, output_centroids=False):
@@ -407 +407 @@
     ##
     # @brief Get interpolated data at given points.
-        #        Applies a transform to all points to calculate the
-        #        interpolated values. Points outside the mesh are returned as NaN.
-        # @note self._A matrix must be valid
+    #        Applies a transform to all points to calculate the
+    #        interpolated values. Points outside the mesh are returned as NaN.
+    # @note self._A matrix must be valid
     # @param f Array of arbitrary data
     # @param verbose True if this function is to be verbose.
@@ -431 +431 @@
     # @brief Build NxM interpolation matrix.
     # @param point_coordinates Points to sample at
-        # @param output_centroids set to True to always sample from the centre
-        #                         of the intersected triangle, instead of the intersection
-        #                         point.
+    # @param output_centroids set to True to always sample from the centre
+    #                         of the intersected triangle, instead of the intersection
+    #                         point.
     # @param verbose True if this function is to be verbose.
     # @return Interpolation matrix A, plus lists of the points inside and outside the mesh
-        #         and the list of centroids, if requested.
+    #         and the list of centroids, if requested.
     def _build_interpolation_matrix_A(self,
                                       point_coordinates,
@@ -464 +464 @@
         if verbose: log.critical('Getting indices inside mesh boundary')
 
-        inside_poly_indices, outside_poly_indices = \
+        # Quick test against boundary, but will not deal with holes in the mesh 
+        inside_boundary_indices, outside_poly_indices = \
             in_and_outside_polygon(point_coordinates,
                                    self.mesh.get_boundary_polygon(),
@@ -474 +475 @@
 
         # Since you can block, throw a warning, not an error.
-        if verbose and 0 == len(inside_poly_indices):
+        if verbose and 0 == len(inside_boundary_indices):
             log.critical('WARNING: No points within the mesh!')
 
@@ -485 +486 @@
         A = Sparse(n,m)
 
-        n = len(inside_poly_indices)
+        n = len(inside_boundary_indices)
 
         centroids = []
-                
+        
+        inside_poly_indices = []
+        
         # Compute matrix elements for points inside the mesh
         if verbose: log.critical('Building interpolation matrix from %d points'
                                  % n)
 
-        for d, i in enumerate(inside_poly_indices):
+        for d, i in enumerate(inside_boundary_indices):
             # For each data_coordinate point
             if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d'
@@ -501 +504 @@
             element_found, sigma0, sigma1, sigma2, k = \
                            search_tree_of_vertices(self.root, x)
+                       
+        # Update interpolation matrix A if necessary
+            if element_found is True:
+
+                if verbose:
+                    print 'Point is within mesh:', d, i            
             
-            # Update interpolation matrix A if necessary
-            if element_found is True:
+                inside_poly_indices.append(i)
+                
                 # Assign values to matrix A
                 j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0
@@ -509 +518 @@
                 j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2
                 js = [j0, j1, j2]
-                                
+                
                 if output_centroids is False:
                     # Weight each vertex according to its distance from x
@@ -516 +525 @@
                         A[i, j] = sigmas[j]
                 else:
-                                    # If centroids are needed, weight all 3 vertices equally
+                    # If centroids are needed, weight all 3 vertices equally
                     for j in js:
                         A[i, j] = 1.0/3.0
-                    centroids.append(self.mesh.centroid_coordinates[k])                                         
+                    centroids.append(self.mesh.centroid_coordinates[k])                        
             else:
-                 # Mesh has a hole - move this point to the outside list
-#                 num.delete(inside_poly_indices,[i], axis=0)
-                 num.append(outside_poly_indices, [i], axis=0)
-#                msg = 'Could not find triangle for point', x
-#                raise Exception(msg)
+                if verbose:
+                    log.critical('Mesh has a hole - moving this point to outside list')
+
+                # This is a numpy arrays, so we need to do a slow transfer
+                outside_poly_indices = num.append(outside_poly_indices, [i], axis=0)
 
         return A, inside_poly_indices, outside_poly_indices, centroids
@@ -775 +784 @@
                  verbose=False,
                  gauge_neighbour_id=None,
-                                 output_centroids=False):
+                 output_centroids=False):
         """Initialise object and build spatial interpolation if required
 
@@ -1003 +1012 @@
                                                       verbose=False,
                                                       output_centroids=output_centroids)
-                        self.centroids = interpol.centroids                                                                                                               
+                        self.centroids = interpol.centroids                                                          
                     elif triangles is None and vertex_coordinates is not None:
                         result = interpolate_polyline(Q,
@@ -1012 +1021 @@
 
                     #assert len(result), len(interpolation_points)
-                    self.precomputed_values[name][i, :] = result                                                                        
-                                        
+                    self.precomputed_values[name][i, :] = result                                    
+                    
             # Report
             if verbose:
-                log.critical(self.statistics())                 
+                log.critical(self.statistics())            
         else:
             # Store quantitites as is

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -107 +107 @@
         assert num.allclose(A.todense(), [[1./3, 1./3, 1./3]])
 
-
+    def test_datapoint_in_hole(self):
+        # create 3 right-angled triangles arranged in a bigger triangle       
+        a = [0.0, 0.0] #0
+        b = [0.0, 2.0] #1
+        c = [2.0,0.0]  #2
+        d = [0.0,4.0]  #3
+        e = [2.0,2.0]  #4
+        f = [4.0,0.0]  #5
+        points = [a, b, c, d, e, f]
+        vertices = [ [1,0,2], [3,1,4], [4,2,5] ]   #bac dbe ecf
+
+        point_in_hole = [1.5, 1.5] # a point in the hole
+        data = [ [20, 20], [0.3, 0.3], point_in_hole, [2.5, 0.3], [30, 30] ] #some points inside and outside the hole
+
+        # any function for the vertices, we don't care about the result
+        f = num.array([linear_function(points), 2*linear_function(points)])
+        f = num.transpose(f)
+
+        interp = Interpolate(points, vertices)
+        interp.interpolate(f, data)
+        
+        assert not set(interp.inside_poly_indices).intersection(set(interp.outside_poly_indices)), \
+               'Some points are in both lists!'    
+        assert len(interp.inside_poly_indices) == 2
+        assert len(interp.outside_poly_indices) == 3
+        
+        interp.outside_poly_indices.sort()
+        assert interp.outside_poly_indices[1] == 2, \
+               'third outside point should be inside the hole!'
 
     def test_simple_interpolation_example(self):
@@ -391 +419 @@
         assert num.allclose(num.sum(results, axis=1), 1.0)
 
-                
+        
     def test_arbitrary_datapoints_return_centroids(self):
         #Try arbitrary datapoints, making sure they return
@@ -409 +437 @@
         third = [1.0/3.0, 1.0/3.0, 1.0/3.0]
         answer = [third, third, third]
-                
+        
         A,_,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
         results = A.todense()
-        assert num.allclose(results, answer)            
-                
-                
+        assert num.allclose(results, answer)        
+        
+        
     def test_arbitrary_datapoints_some_outside(self):
         #Try arbitrary datapoints one outside the triangle.

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -2258 +2258 @@
     import types
 
-    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
-         separate_points_by_polygon
+    from anuga.utilities.polygon import inside_polygon, outside_polygon
     from anuga.abstract_2d_finite_volumes.util import \
          apply_expression_to_dictionary
@@ -2690 +2689 @@
 
     import sys
-    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
-             separate_points_by_polygon
+    from anuga.utilities.polygon import inside_polygon, outside_polygon
     from anuga.abstract_2d_finite_volumes.util import \
              apply_expression_to_dictionary

anuga_core/source/anuga/shallow_water/data_manager_joaquims_patch.py

@@ -1847 +1847 @@
     from Numeric import array2string
 
-    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
-         separate_points_by_polygon
+    from anuga.utilities.polygon import inside_polygon, outside_polygon
     from anuga.abstract_2d_finite_volumes.util import \
          apply_expression_to_dictionary
@@ -2232 +2231 @@
     from Numeric import array2string
 
-    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
+    from anuga.utilities.polygon import inside_polygon, outside_polygon
     from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary