Changeset 7675


Ignore:
Timestamp:
Apr 6, 2010, 8:23:54 PM (15 years ago)
Author:
James Hudson
Message:

Ticket 113 is complete, and all tests pass.
A centroid list is built by Interpolation_function as it calculates the interpolation matrix, and this is passed out as extra quantities which are output into the gauge.csv file.

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

    r7673 r7675  
    190190                quake_time = time + quake_offset_time
    191191                points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
    192                 point_quantities = callable_sww(time,point_i) # __call__ is overridden
    193                                
     192                point_quantities = callable_sww(time, point_i) # __call__ is overridden
     193                           
    194194                for quantity in quantities:
    195195                    if quantity == NAN:
    … …  
    197197                                     % callable_sww.get_name)
    198198                    else:
     199                        #core quantities that are exported from the interpolator     
    199200                        if quantity == 'stage':
    200201                            points_list.append(point_quantities[0])
    … …  
    208209                        if quantity == 'ymomentum':
    209210                            points_list.append(point_quantities[3])
    210                            
     211
     212                        #derived quantities that are calculated from the core ones
    211213                        if quantity == 'depth':
    212214                            points_list.append(point_quantities[0]
    … …  
    237239                            points_list.append(calc_bearing(point_quantities[2],
    238240                                                            point_quantities[3]))
    239 
     241                        if quantity == 'xcentroid':
     242                            points_list.append(callable_sww.centroids[point_i][0])
     243
     244                        if quantity == 'ycentroid':
     245                            points_list.append(callable_sww.centroids[point_i][1])
     246                                                       
    240247                points_writer[point_i].writerow(points_list)
    241248

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

    r7673 r7675  
    451451
    452452        # create a csv file containing our gauge points
    453         points = [[2.0,1.0],[4.5,4.0]]
    454 
    455453        points_file = tempfile.mktemp(".csv")
    456454        file_id = open(points_file,"w")
    … …  
    510508
    511509
     510    def test_sww2csv_output_centroid_attribute(self):
     511
     512        def elevation_function(x, y):
     513            return -x
     514       
     515        """Check sww2csv timeseries at centroid, then output the centroid coordinates.
     516       
     517        Test the ability to get a timeseries at the centroid of a triangle, rather
     518                than the given gauge point, then output the results.
     519        """
     520       
     521        # Create rectangular mesh
     522        mesh_file = tempfile.mktemp(".tsh")   
     523        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
     524        m = Mesh()
     525        m.add_vertices(points)
     526        m.auto_segment()
     527        m.generate_mesh(verbose=False)
     528        m.export_mesh_file(mesh_file)
     529       
     530        # Create shallow water domain
     531        domain = Domain(mesh_file)
     532       
     533        domain.default_order=2
     534       
     535        # This test was made before tight_slope_limiters were introduced
     536        # Since were are testing interpolation values this is OK
     537        domain.tight_slope_limiters = 0
     538       
     539
     540        # Set some field values
     541        domain.set_quantity('elevation', elevation_function)
     542        domain.set_quantity('friction', 0.03)
     543
     544        ######################
     545        # Boundary conditions
     546        B = Transmissive_boundary(domain)
     547        domain.set_boundary( {'exterior': B})
     548
     549        # This call mangles the stage values.
     550        domain.distribute_to_vertices_and_edges()
     551        domain.set_quantity('stage', 1.0)
     552
     553
     554        domain.set_name('datatest' + str(time.time()))
     555        domain.smooth = True
     556        domain.reduction = mean
     557
     558
     559        sww = SWW_file(domain)
     560        sww.store_connectivity()
     561        sww.store_timestep()
     562        domain.set_quantity('stage', 10.0) # This is automatically limited
     563        # so it will not be less than the elevation
     564        domain.time = 2.
     565        sww.store_timestep()
     566
     567        # create a csv file containing our gauge points
     568        points_file = tempfile.mktemp(".csv")
     569        file_id = open(points_file,"w")
     570 
     571# These values are slightly off the centroids - will it find the centroids?
     572        file_id.write("name, easting, northing, elevation \n\
     573point1, 2.5, 4.25, 3.0\n")
     574
     575        file_id.close()
     576
     577        gauge_sww2csv(sww.filename,
     578                       points_file,
     579                       quantities=['stage', 'xcentroid', 'ycentroid'],
     580                       verbose=False,
     581                       use_cache=False,
     582                       output_centroids=True)
     583
     584        point1_answers_array = [[0.0,0.0,1.0,4.0,4.0], [2.0,2.0/3600.,10.0,4.0,4.0]]
     585        point1_filename = 'gauge_point1.csv'
     586        point1_handle = file(point1_filename)
     587        point1_reader = reader(point1_handle)
     588        point1_reader.next()
     589
     590        line=[]
     591        for i,row in enumerate(point1_reader):
     592            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),float(row[4])])
     593#            print 'assert line',line[i],'point1',point1_answers_array[i]
     594            assert num.allclose(line[i], point1_answers_array[i])
     595
     596        # clean up
     597        point1_handle.close()
     598        os.remove(mesh_file)           
     599        os.remove(sww.filename)
     600        os.remove(points_file)
     601        os.remove(point1_filename)
     602               
     603               
    512604
    513605#-------------------------------------------------------------

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r7673 r7675  
    502502                                   verbose=verbose,
    503503                                   gauge_neighbour_id=gauge_neighbour_id,
    504                                                                    output_centroids=output_centroids),
     504                                   output_centroids=output_centroids),
    505505            starttime)
    506506

  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r7673 r7675  
    218218
    219219
    220 
    221220    ##
    222221    # @brief Interpolate mesh data f to determine values, z, at points.
    … …  
    383382
    384383        # Unpack result
    385         self._A, self.inside_poly_indices, self.outside_poly_indices = X
     384        self._A, self.inside_poly_indices, self.outside_poly_indices, self.centroids = X
    386385
    387386        # Check that input dimensions are compatible
    … …  
    436435        #                         point.
    437436    # @param verbose True if this function is to be verbose.
    438     # @return Interpolation matrix A, plus lists of the points inside and outside the mesh.
     437    # @return Interpolation matrix A, plus lists of the points inside and outside the mesh
     438        #         and the list of centroids, if requested.
    439439    def _build_interpolation_matrix_A(self,
    440440                                      point_coordinates,
    … …  
    487487        n = len(inside_poly_indices)
    488488
     489        centroids = []
     490               
    489491        # Compute matrix elements for points inside the mesh
    490492        if verbose: log.critical('Building interpolation matrix from %d points'
    … …  
    516518                                    # If centroids are needed, weight all 3 vertices equally
    517519                    for j in js:
    518                         A[i, j] = 1.0/3.0                   
     520                        A[i, j] = 1.0/3.0
     521                    centroids.append(self.mesh.centroid_coordinates[k])                                         
    519522            else:
    520523                msg = 'Could not find triangle for point', x
    521524                raise Exception(msg)
    522         return A, inside_poly_indices, outside_poly_indices
     525        return A, inside_poly_indices, outside_poly_indices, centroids
    523526
    524527
    … …  
    841844        self.index = 0    # Initial time index
    842845        self.precomputed_values = {}
     846        self.centroids = []
    843847
    844848        # Precomputed spatial interpolation if requested
    … …  
    994998                                                      self.interpolation_points,
    995999                                                      verbose=False,
    996                                                                                                           output_centroids=output_centroids)
     1000                                                      output_centroids=output_centroids)
     1001                        self.centroids = interpol.centroids                                                                                                               
    9971002                    elif triangles is None and vertex_coordinates is not None:
    9981003                        result = interpolate_polyline(Q,
    … …  
    10031008
    10041009                    #assert len(result), len(interpolation_points)
    1005                     self.precomputed_values[name][i, :] = result
    1006 
     1010                    self.precomputed_values[name][i, :] = result                                                                       
     1011                                       
    10071012            # Report
    10081013            if verbose:
    1009                 log.critical(self.statistics())
     1014                log.critical(self.statistics())                 
    10101015        else:
    10111016            # Store quantitites as is

  • anuga_core/source/anuga/fit_interpolate/test_interpolate.py

    r7673 r7675  
    104104
    105105        interp = Interpolate(points, vertices)
    106         A, _, _ = interp._build_interpolation_matrix_A(data)
     106        A, _, _, _ = interp._build_interpolation_matrix_A(data)
    107107        assert num.allclose(A.todense(), [[1./3, 1./3, 1./3]])
    108108
    … …  
    285285                      0., 0. , 0., 0., 0., 0.]]
    286286
    287         A,_,_ = interp._build_interpolation_matrix_A(data)
     287        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    288288        assert num.allclose(A.todense(), answer)
    289289       
    … …  
    294294                      0., 0. , 0., 0., 0., 0.]]
    295295       
    296         A,_,_ = interp._build_interpolation_matrix_A(data)       
     296        A,_,_,_ = interp._build_interpolation_matrix_A(data)       
    297297        assert num.allclose(A.todense(), answer)
    298298
    … …  
    305305                      0., 0. , 0., 0., 0., 0.]]
    306306                     
    307         A,_,_ = interp._build_interpolation_matrix_A(data)       
     307        A,_,_,_ = interp._build_interpolation_matrix_A(data)       
    308308        assert num.allclose(A.todense(), answer)
    309309
    … …  
    326326                   [0., 0., 1.]]
    327327                   
    328         A,_,_ = interp._build_interpolation_matrix_A(data)
     328        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    329329        assert num.allclose(A.todense(), answer)
    330330
    … …  
    347347        interp = Interpolate(points, vertices)
    348348
    349         A,_,_ = interp._build_interpolation_matrix_A(data)
     349        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    350350        assert num.allclose(A.todense(), answer)
    351351
    … …  
    368368        interp = Interpolate(points, vertices)
    369369
    370         A,_,_ = interp._build_interpolation_matrix_A(data)
     370        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    371371        assert num.allclose(A.todense(), answer)
    372372
    … …  
    387387        #print "interp.get_A()", interp.get_A()
    388388       
    389         A,_,_ = interp._build_interpolation_matrix_A(data)
     389        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    390390        results = A.todense()
    391391        assert num.allclose(num.sum(results, axis=1), 1.0)
    … …  
    410410        answer = [third, third, third]
    411411               
    412         A,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
     412        A,_,_,_ = interp._build_interpolation_matrix_A(data, output_centroids=True)
    413413        results = A.todense()
    414414        assert num.allclose(results, answer)           
    … …  
    430430        interp = Interpolate(points, vertices)
    431431       
    432         A,_,_ = interp._build_interpolation_matrix_A(data)
     432        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    433433        results = A.todense()
    434434        assert num.allclose(num.sum(results, axis=1), [1,1,1,0])
    … …  
    461461
    462462
    463         A,_,_ = interp._build_interpolation_matrix_A(data)
     463        A,_,_,_ = interp._build_interpolation_matrix_A(data)
    464464        A = A.todense()
    465465        for i in range(A.shape[0]):
Note: See TracChangeset for help on using the changeset viewer.