Changeset 4779


Ignore:
Timestamp:
Nov 2, 2007, 5:18:29 PM (17 years ago)
Author:
duncan
Message:

reduce memory use in quantity.set_value. fit_to_mesh can now use an existing mesh instance, which it does in quantity.set_value.

Location:
anuga_core/source/anuga
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r4599 r4779  
    259259               
    260260        return V
    261        
     261   
     262    def get_node(self, i,
     263                 absolute=False):
     264        """Return node coordinates.
     265
     266        The nodes are ordered in an Nx2 array where N is the number of nodes.
     267        This is the same format they were provided in the constructor
     268        i.e. without any duplication.
     269
     270        Boolean keyword argument absolute determines whether coordinates
     271        are to be made absolute by taking georeference into account
     272        Default is False as many parts of ANUGA expects relative coordinates.
     273        (To see which, switch to default absolute=True and run tests).       
     274        """
     275
     276       
     277        V = self.nodes[i,:]
     278        if absolute is True:
     279            if not self.geo_reference.is_absolute():
     280                V[0] += self.geo_reference.xllcorner
     281                V[1] += self.geo_reference.yllcorner
     282               
     283        return V
     284   
    262285       
    263286

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

    r4769 r4779  
    726726
    727727        # Call fit_interpolate.fit function
    728         args = (coordinates, triangles, points, values)
    729         kwargs = {'data_origin': data_georef.get_origin(),
     728        #args = (coordinates, triangles, points, values)
     729        args = (points, )
     730        kwargs = {'vertex_coordinates': coordinates,
     731                  'triangles': triangles,
     732                  'mesh': None,
     733                  'point_attributes': values,
     734                  'data_origin': data_georef.get_origin(),
    730735                  'mesh_origin': mesh_georef.get_origin(),
    731736                  'alpha': alpha,
    … …  
    776781            raise msg
    777782
    778         coordinates = self.domain.get_nodes(absolute=True)
    779         triangles = self.domain.triangles      #FIXME
    780            
    781         vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
     783        #coordinates = self.domain.get_nodes(absolute=True)
     784        #triangles = self.domain.triangles      #FIXME
     785        #vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
     786        vertex_attributes = fit_to_mesh(filename,
     787                                        mesh = self.domain,
    782788                                        alpha=alpha,
    783789                                        attribute_name=attribute_name,

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

    r4599 r4779  
    99from Numeric import allclose, array, ones, Float
    1010from general_mesh import General_mesh
     11from anuga.coordinate_transforms.geo_reference import Geo_reference
    1112
    1213
    … …  
    216217        assert unique_vertices == [0,2,4,5,6,7]
    217218
    218 
    219        
     219    def test_get_node(self):
     220        """test_get_triangles_and_vertices_per_node -
     221
     222        Test that tuples of triangle, vertex can be extracted
     223        from inverted triangles structure
     224
     225        """
     226       
     227        x0 = 314036.58727982
     228        y0 = 6224951.2960092
     229        geo = Geo_reference(56, x0, y0)
     230       
     231        a = [0.0, 0.0]
     232        b = [0.0, 2.0]
     233        c = [2.0, 0.0]
     234        d = [0.0, 4.0]
     235        e = [2.0, 2.0]
     236        f = [4.0, 0.0]
     237
     238        nodes = array([a, b, c, d, e, f])
     239
     240        nodes_absolute = geo.get_absolute(nodes)
     241       
     242        #bac, bce, ecf, dbe, daf, dae
     243        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
     244
     245        domain = General_mesh(nodes, triangles,
     246                       geo_reference = geo)
     247        node = domain.get_node(2)       
     248        self.assertEqual(c, node)
     249       
     250        node = domain.get_node(2, absolute=True)     
     251        self.assertEqual(nodes_absolute[2], node)
     252       
     253
    220254
    221255#-------------------------------------------------------------
    222256if __name__ == "__main__":
    223     suite = unittest.makeSuite(Test_General_Mesh,'test')   
     257    suite = unittest.makeSuite(Test_General_Mesh,'test')
     258    #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')   
    224259    runner = unittest.TextTestRunner()
    225260    runner.run(suite)

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

    r4712 r4779  
    541541
    542542        #Now try by setting the same values directly
    543         vertex_attributes = fit_to_mesh(quantity.domain.get_nodes(),
     543        vertex_attributes = fit_to_mesh(data_points,
     544                                        quantity.domain.get_nodes(),
    544545                                        quantity.domain.triangles, #FIXME
    545                                         data_points,
    546                                         z,
     546                                        point_attributes=z,
    547547                                        alpha = 0,
    548548                                        verbose=False)
    … …  
    584584
    585585        #Reference
    586         ref = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
     586        ref = fit_to_mesh(data_points, vertex_coordinates, triangles,
     587                          point_attributes=z,
    587588                          data_origin = data_georef.get_origin(),
    588589                          mesh_origin = mesh_georef.get_origin(),
    … …  
    10481049        x0 = 314036.58727982
    10491050        y0 = 6224951.2960092
     1051        #x0 = 0.0
     1052        #y0 = 0.0
    10501053
    10511054        a = [0.0, 0.0]

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

    r4663 r4779  
    2727import types
    2828
    29 from Numeric import zeros, Float, ArrayType,take
    30 
     29from Numeric import zeros, Float, ArrayType,take, Int
     30
     31from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
    3132from anuga.caching import cache           
    3233from anuga.geospatial_data.geospatial_data import Geospatial_data, \
    … …  
    4950   
    5051    def __init__(self,
    51                  vertex_coordinates,
    52                  triangles,
     52                 vertex_coordinates=None,
     53                 triangles=None,
     54                 mesh=None,
    5355                 mesh_origin=None,
    5456                 alpha = None,
    … …  
    98100                 vertex_coordinates,
    99101                 triangles,
     102                 mesh,
    100103                 mesh_origin,
    101104                 verbose,
    … …  
    434437############################################################################
    435438
    436 def fit_to_mesh(vertex_coordinates,
    437                 triangles,
    438                 point_coordinates, # this can also be a points file name
     439def fit_to_mesh(point_coordinates, # this can also be a points file name
     440                vertex_coordinates=None,
     441                triangles=None,
     442                mesh=None,
    439443                point_attributes=None,
    440444                alpha=DEFAULT_ALPHA,
    … …  
    450454    """
    451455   
    452     args = (vertex_coordinates, triangles, point_coordinates, )
    453     kwargs = {'point_attributes': point_attributes,
     456    args = (point_coordinates, )
     457    kwargs = {'vertex_coordinates': vertex_coordinates,
     458              'triangles': triangles,
     459              'mesh': mesh,
     460              'point_attributes': point_attributes,
    454461              'alpha': alpha,
    455462              'verbose': verbose,
    … …  
    480487                     args, kwargs)
    481488
    482 def _fit_to_mesh(vertex_coordinates,
    483                  triangles,
    484                  point_coordinates, # this can also be a points file name
     489def _fit_to_mesh(point_coordinates, # this can also be a points file name
     490                 vertex_coordinates=None,
     491                 triangles=None,
     492                 mesh=None,
    485493                 point_attributes=None,
    486494                 alpha=DEFAULT_ALPHA,
    … …  
    513521          alpha: Smoothing parameter.
    514522
    515           acceptable overshoot: controls the allowed factor by which fitted values
     523          acceptable overshoot: NOT IMPLEMENTED
     524          controls the allowed factor by which
     525          fitted values
    516526          may exceed the value of input data. The lower limit is defined
    517527          as min(z) - acceptable_overshoot*delta z and upper limit
    518528          as max(z) + acceptable_overshoot*delta z
     529         
    519530
    520531          mesh_origin: A geo_reference object or 3-tuples consisting of
    … …  
    531542    # Duncan and Ole think that this isn't worth caching.
    532543    # Caching happens at the higher level anyway.
    533     interp = Fit(vertex_coordinates,
    534                  triangles,
     544   
     545    if mesh is None:
     546        # Fixme (DSG) Throw errors if triangles or vertex_coordinates
     547        # are None
     548           
     549        #Convert input to Numeric arrays
     550        triangles = ensure_numeric(triangles, Int)
     551        vertex_coordinates = ensure_absolute(vertex_coordinates,
     552                                             geo_reference = mesh_origin)
     553
     554        if verbose: print 'FitInterpolate: Building mesh'       
     555        mesh = Mesh(vertex_coordinates, triangles)
     556        mesh.check_integrity()
     557   
     558    interp = Fit(mesh=mesh,
    535559                 verbose=verbose,
    536                  mesh_origin=mesh_origin,
    537560                 alpha=alpha)
    538561
    … …  
    624647    if verbose: print "points file loaded"
    625648    if verbose: print "fitting to mesh"
    626     f = fit_to_mesh(vertex_coordinates,
     649    f = fit_to_mesh(point_coordinates,
     650                    vertex_coordinates,
    627651                    triangles,
    628                     point_coordinates,
     652                    None,
    629653                    point_attributes,
    630654                    alpha = alpha,

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

    r4739 r4779  
    4444   
    4545    def __init__(self,
    46                  vertex_coordinates,
    47                  triangles,
     46                 vertex_coordinates=None,
     47                 triangles=None,
     48                 mesh=None,
    4849                 mesh_origin=None,
    4950                 verbose=False,
    … …  
    5455        function values at vertices to function values at data points
    5556
     57        Pass in a mesh instance or vertex_coordinates and triangles
     58        and optionally mesh_origin
     59       
    5660        Inputs:
    5761
    … …  
    6468          triangles: List of 3-tuples (or a Numeric array) of
    6569              integers representing indices of all vertices in the mesh.
     70
     71        mesh: A mesh instance describing the mesh.
    6672
    6773          mesh_origin: A geo_reference object or 3-tuples consisting of
    … …  
    7985        if max_vertices_per_cell == None:
    8086            max_vertices_per_cell = MAX_VERTICES_PER_CELL
    81        
    82         #Convert input to Numeric arrays
    83         triangles = ensure_numeric(triangles, Int)
    84         vertex_coordinates = ensure_absolute(vertex_coordinates,
    85                                              geo_reference = mesh_origin)
    8687
    87         #Don't pass geo_reference to mesh.  It doesn't work. (Still??)
    88        
    89         if verbose: print 'FitInterpolate: Building mesh'       
    90         self.mesh = Mesh(vertex_coordinates, triangles)
    91         self.mesh.check_integrity()
     88        if mesh is None:
     89            # Fixme (DSG) Throw errors if triangles or vertex_coordinates
     90            # are None
     91           
     92            #Convert input to Numeric arrays
     93            triangles = ensure_numeric(triangles, Int)
     94            vertex_coordinates = ensure_absolute(vertex_coordinates,
     95                                                 geo_reference = mesh_origin)
     96
     97            if verbose: print 'FitInterpolate: Building mesh'       
     98            self.mesh = Mesh(vertex_coordinates, triangles)
     99            self.mesh.check_integrity()
     100        else:
     101            self.mesh = mesh
    92102       
    93103        if verbose: print 'FitInterpolate: Building quad tree'

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

    r4663 r4779  
    8484       
    8585        FitInterpolate.__init__(self,
    86                                 vertex_coordinates,
    87                                 triangles,
    88                                 mesh_origin,
    89                                 verbose,
    90                                 max_vertices_per_cell)
     86                                vertex_coordinates=vertex_coordinates,
     87                                triangles=triangles,
     88                                mesh_origin=mesh_origin,
     89                                verbose=verbose,
     90                                max_vertices_per_cell=max_vertices_per_cell)
    9191
    9292    # FIXME: What is a good start_blocking_len value?

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

    r4659 r4779  
    270270        points = [a, b, c]
    271271        elements = [[0,2,1]]
    272         f = fit_to_mesh(points, elements, txt_file,
     272        f = fit_to_mesh(txt_file, points, elements,
    273273                        alpha=0.0, max_read_lines=2)
    274274        answer = linear_function(points)
    … …  
    312312        points = geo.get_data_points(absolute=True)
    313313        atts = geo.get_attributes()
    314         f = fit_to_mesh(vertices, triangles,points,atts,
     314        f = fit_to_mesh(points,atts, vertices, triangles,
    315315                                alpha=0.0, max_read_lines=2,
    316316                        use_cache=True, verbose=True)
    … …  
    352352        fileName_pts = tempfile.mktemp(".pts")
    353353        geo.export_points_file(fileName_pts)
    354         f = fit_to_mesh(vertices, triangles,fileName_pts,
     354        f = fit_to_mesh(fileName_pts, vertices, triangles,
     355                                alpha=0.0, max_read_lines=2)
     356        answer = linear_function(vertices)
     357        #print "f\n",f
     358        #print "answer\n",answer
     359        assert allclose(f, answer)
     360        os.remove(fileName)
     361        os.remove(fileName_pts)
     362       
     363    def test_fit_to_mesh_pts_passing_mesh_in(self):
     364        a = [-1.0, 0.0]
     365        b = [3.0, 4.0]
     366        c = [4.0,1.0]
     367        d = [-3.0, 2.0] #3
     368        e = [-1.0,-2.0]
     369        f = [1.0, -2.0] #5
     370
     371        vertices = [a, b, c, d,e,f]
     372        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
     373
     374
     375        fileName = tempfile.mktemp(".txt")
     376        file = open(fileName,"w")
     377        file.write(" x, y, elevation \n\
     378-2.0, 2.0, 0.\n\
     379-1.0, 1.0, 0.\n\
     3800.0, 2.0 , 2.\n\
     3811.0, 1.0 , 2.\n\
     382 2.0,  1.0 ,3. \n\
     383 0.0,  0.0 , 0.\n\
     384 1.0,  0.0 , 1.\n\
     385 0.0,  -1.0, -1.\n\
     386 -0.2, -0.5, -0.7\n\
     387 -0.9, -1.5, -2.4\n\
     388 0.5,  -1.9, -1.4\n\
     389 3.0,  1.0 , 4.\n")
     390        file.close()
     391        geo = Geospatial_data(fileName)
     392        fileName_pts = tempfile.mktemp(".pts")
     393        geo.export_points_file(fileName_pts)
     394        f = fit_to_mesh(fileName_pts, vertices, triangles,
    355395                                alpha=0.0, max_read_lines=2)
    356396        answer = linear_function(vertices)
    … …  
    391431        file.close()
    392432       
    393         f = fit_to_mesh(vertices, triangles,fileName,
     433        f = fit_to_mesh(fileName, vertices, triangles,
    394434                                alpha=0.0, max_read_lines=2)
    395435                        #use_cache=True, verbose=True)
    … …  
    432472        file.close()
    433473       
    434         f = fit_to_mesh(vertices, triangles,fileName,
     474        f = fit_to_mesh(fileName, vertices, triangles,
    435475                        alpha=0.0,
    436476                        attribute_name='elevation', max_read_lines=2)
    … …  
    686726        z = [z1, z2, z3]
    687727
    688         f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
     728        f = fit_to_mesh(data_coords, vertex_coordinates=points,
     729                        triangles=triangles,
     730                        point_attributes=z, alpha=0.0)
    689731        answer = [[0, 0], [5., 10.], [5., 10.]]
    690732
    … …  
    720762       
    721763        #Fit
    722         zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
     764        zz = fit_to_mesh(data_points, vertex_coordinates=vertex_coordinates,
     765                         triangles=triangles,
     766                         point_attributes=z,
    723767                         data_origin = data_geo.get_origin(),
    724768                         mesh_origin = mesh_geo.get_origin(),
    … …  
    10521096if __name__ == "__main__":
    10531097    suite = unittest.makeSuite(Test_Fit,'test')
    1054     #suite = unittest.makeSuite(Test_Fit,'cache_test_fit_to_mesh_pts')
     1098    #suite = unittest.makeSuite(Test_Fit,'test_fit_to_mesh_pts_passing_mesh_in')
    10551099    #suite = unittest.makeSuite(Test_Fit,'')
    10561100    runner = unittest.TextTestRunner(verbosity=1)

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

    r4739 r4779  
    1414
    1515from Scientific.IO.NetCDF import NetCDFFile
    16 from Numeric import allclose, array, transpose, zeros, Float, sometrue, alltrue, take, where
     16from Numeric import allclose, array, transpose, zeros, Float, sometrue, \
     17     alltrue, take, where
    1718
    1819

  • anuga_core/source/anuga/utilities/quad.py

    r4739 r4779  
    142142        if len(args) == 2:
    143143            point_id = int(args[1])
    144             x, y = self.mesh.get_nodes()[point_id]
     144            x, y = self.mesh.get_node(point_id, absolute=True)
    145145
    146146            #print point_id, x, y
    … …  
    183183                    self.store(point)
    184184                else:
    185                     x = self.mesh.coordinates[point][0]
    186                     y = self.mesh.coordinates[point][1]
    187                     print "(" + str(x) + "," + str(y) + ")"
     185                    # Have to take into account of georef.
     186                    #x = self.mesh.coordinates[point][0]
     187                    #y = self.mesh.coordinates[point][1]
     188                    node = self.mesh.get_node(point, absolute=True)
     189                    print "node", node
     190                    print "(" + str(node[0]) + "," + str(node[1]) + ")"
    188191                    raise 'point not in region: %s' %str(point)
    189192               
    … …  
    229232                    if not triangles.has_key(k):
    230233                        # print 'k',k
    231                         tri = self.mesh.get_vertex_coordinates(k)
     234                        tri = self.mesh.get_vertex_coordinates(k,
     235                                                               absolute=True)
    232236                        n0 = self.mesh.get_normal(k, 0)
    233237                        n1 = self.mesh.get_normal(k, 1)
    … …  
    441445    #print mesh.coordinates
    442446
    443     nodes = mesh.get_nodes()
     447   
     448    nodes = mesh.get_nodes(absolute=True)
    444449    xmin = min(nodes[:,0])
    445450    xmax = max(nodes[:,0])
    … …  
    447452    ymax = max(nodes[:,1])
    448453
     454    # Don't know why this didn't work
     455    #xmin, xmax, ymin, ymax = mesh.get_extent(absolute=True)
    449456   
    450457    # Ensure boundary points are fully contained in region
Note: See TracChangeset for help on using the changeset viewer.