Changeset 6226


Ignore:
Timestamp:
Jan 21, 2009, 5:28:57 PM (10 years ago)
Author:
rwilson
Message:

Did PEP8 pass, @brief.

Location:
anuga_core/source/anuga/abstract_2d_finite_volumes
Files:
2 edited

Legend:

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

    r6191 r6226  
    88"""
    99
     10import types
     11from time import time as walltime
     12
    1013from anuga.config import epsilon
    1114from anuga.config import beta_euler, beta_rk2
    12 
    1315from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
    1416from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     
    2224from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    2325     import Transmissive_boundary
    24 
    2526from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
    2627from anuga.abstract_2d_finite_volumes.region\
    2728     import Set_region as region_set_region
    28 
    2929from anuga.utilities.polygon import inside_polygon
    3030from anuga.abstract_2d_finite_volumes.util import get_textual_float
    31 
    3231from quantity import Quantity
    3332
    34 import types
    35 from time import time as walltime
    36 
    3733import Numeric as num
    3834
    3935
    4036##
    41 # @brief Basic Domain class
     37# @brief Generic Domain class
    4238class Domain:
    4339
     
    4945    # @param conserved_quantities List of names of quantities to be conserved.
    5046    # @param other_quantities List of names of other quantities.
    51     # @param tagged_elements
    52     # @param geo_reference
    53     # @param use_inscribed_circle
    54     # @param mesh_filename
    55     # @param use_cache
    56     # @param verbose
    57     # @param full_send_dict
    58     # @param ghost_recv_dict
    59     # @param processor
    60     # @param numproc
    61     # @param number_of_full_nodes
    62     # @param number_of_full_triangles
     47    # @param tagged_elements ??
     48    # @param geo_reference ??
     49    # @param use_inscribed_circle ??
     50    # @param mesh_filename ??
     51    # @param use_cache ??
     52    # @param verbose True if this method is to be verbose.
     53    # @param full_send_dict ??
     54    # @param ghost_recv_dict ??
     55    # @param processor ??
     56    # @param numproc ??
     57    # @param number_of_full_nodes ??
     58    # @param number_of_full_triangles ??
    6359    def __init__(self, source=None,
    6460                       triangles=None,
     
    119115                         number_of_full_triangles=number_of_full_triangles,
    120116                         verbose=verbose)
    121                          
     117
    122118        # Expose Mesh attributes (FIXME: Maybe turn into methods)
    123119        self.centroid_coordinates = self.mesh.centroid_coordinates
    124         self.vertex_coordinates = self.mesh.vertex_coordinates       
     120        self.vertex_coordinates = self.mesh.vertex_coordinates
    125121        self.boundary = self.mesh.boundary
    126122        self.neighbours = self.mesh.neighbours
    127         self.surrogate_neighbours = self.mesh.surrogate_neighbours       
     123        self.surrogate_neighbours = self.mesh.surrogate_neighbours
    128124        self.neighbour_edges = self.mesh.neighbour_edges
    129125        self.normals = self.mesh.normals
    130         self.edgelengths = self.mesh.edgelengths       
    131         self.radii = self.mesh.radii               
    132         self.areas = self.mesh.areas                       
    133                
    134         self.number_of_boundaries = self.mesh.number_of_boundaries       
     126        self.edgelengths = self.mesh.edgelengths
     127        self.radii = self.mesh.radii
     128        self.areas = self.mesh.areas
     129
     130        self.number_of_boundaries = self.mesh.number_of_boundaries
    135131        self.number_of_full_nodes = self.mesh.number_of_full_nodes
    136         self.number_of_full_triangles = self.mesh.number_of_full_triangles       
     132        self.number_of_full_triangles = self.mesh.number_of_full_triangles
    137133        self.number_of_triangles_per_node = self.mesh.number_of_triangles_per_node
    138134
    139135        self.vertex_value_indices = self.mesh.vertex_value_indices
    140         self.number_of_triangles = self.mesh.number_of_triangles       
     136        self.number_of_triangles = self.mesh.number_of_triangles
    141137
    142138        self.geo_reference = self.mesh.geo_reference
    143        
     139
    144140        if verbose: print 'Initialising Domain'
    145141
     
    276272
    277273        if verbose: print 'Domain: Done'
    278        
    279      
    280 
     274
     275    ######
    281276    # Expose underlying Mesh functionality
     277    ######
     278
    282279    def __len__(self):
    283280        return len(self.mesh)
     
    285282    def get_centroid_coordinates(self, *args, **kwargs):
    286283        return self.mesh.get_centroid_coordinates(*args, **kwargs)
    287        
     284
    288285    def get_radii(self, *args, **kwargs):
    289         return self.mesh.get_radii(*args, **kwargs)       
    290        
     286        return self.mesh.get_radii(*args, **kwargs)
     287
    291288    def get_areas(self, *args, **kwargs):
    292         return self.mesh.get_areas(*args, **kwargs)               
     289        return self.mesh.get_areas(*args, **kwargs)
    293290
    294291    def get_area(self, *args, **kwargs):
     
    296293
    297294    def get_vertex_coordinates(self, *args, **kwargs):
    298         return self.mesh.get_vertex_coordinates(*args, **kwargs)               
     295        return self.mesh.get_vertex_coordinates(*args, **kwargs)
    299296
    300297    def get_triangles(self, *args, **kwargs):
    301         return self.mesh.get_triangles(*args, **kwargs)               
    302        
     298        return self.mesh.get_triangles(*args, **kwargs)
     299
    303300    def get_nodes(self, *args, **kwargs):
    304301        return self.mesh.get_nodes(*args, **kwargs)
    305        
     302
    306303    def get_number_of_nodes(self, *args, **kwargs):
    307304        return self.mesh.get_number_of_nodes(*args, **kwargs)
    308        
     305
    309306    def get_normal(self, *args, **kwargs):
    310         return self.mesh.get_normal(*args, **kwargs)       
    311        
     307        return self.mesh.get_normal(*args, **kwargs)
     308
    312309    def get_intersecting_segments(self, *args, **kwargs):
    313310        return self.mesh.get_intersecting_segments(*args, **kwargs)
    314        
     311
    315312    def get_disconnected_triangles(self, *args, **kwargs):
    316313        return self.mesh.get_disconnected_triangles(*args, **kwargs)
    317        
     314
    318315    def get_boundary_tags(self, *args, **kwargs):
    319316        return self.mesh.get_boundary_tags(*args, **kwargs)
     
    321318    def get_boundary_polygon(self, *args, **kwargs):
    322319        return self.mesh.get_boundary_polygon(*args, **kwargs)
    323                
     320
    324321    def get_number_of_triangles_per_node(self, *args, **kwargs):
    325322        return self.mesh.get_number_of_triangles_per_node(*args, **kwargs)
    326        
     323
    327324    def get_triangles_and_vertices_per_node(self, *args, **kwargs):
    328325        return self.mesh.get_triangles_and_vertices_per_node(*args, **kwargs)
    329        
     326
    330327    def get_interpolation_object(self, *args, **kwargs):
    331         return self.mesh.get_interpolation_object(*args, **kwargs)       
    332        
     328        return self.mesh.get_interpolation_object(*args, **kwargs)
     329
    333330    def get_tagged_elements(self, *args, **kwargs):
    334         return self.mesh.get_tagged_elements(*args, **kwargs)               
    335        
     331        return self.mesh.get_tagged_elements(*args, **kwargs)
     332
    336333    def get_lone_vertices(self, *args, **kwargs):
    337         return self.mesh.get_lone_vertices(*args, **kwargs)       
    338        
     334        return self.mesh.get_lone_vertices(*args, **kwargs)
     335
    339336    def get_unique_vertices(self, *args, **kwargs):
    340         return self.mesh.get_unique_vertices(*args, **kwargs)               
     337        return self.mesh.get_unique_vertices(*args, **kwargs)
    341338
    342339    def get_georeference(self, *args, **kwargs):
    343340        return self.mesh.get_georeference(*args, **kwargs)
    344        
     341
    345342    def set_georeference(self, *args, **kwargs):
    346         self.mesh.set_georeference(*args, **kwargs)                   
    347        
     343        self.mesh.set_georeference(*args, **kwargs)
     344
    348345    def build_tagged_elements_dictionary(self, *args, **kwargs):
    349346        self.mesh.build_tagged_elements_dictionary(*args, **kwargs)
    350        
     347
    351348    def statistics(self, *args, **kwargs):
    352         return self.mesh.statistics(*args, **kwargs)       
    353                
    354        
    355        
     349        return self.mesh.statistics(*args, **kwargs)
     350
    356351    ##
    357352    # @brief Get conserved quantities for a volume.
     
    365360                                       vertex=None,
    366361                                       edge=None):
    367         """Get conserved quantities at volume vol_id
     362        """Get conserved quantities at volume vol_id.
    368363
    369364        If vertex is specified use it as index for vertex values
     
    397392    # @param time The new model time (seconds).
    398393    def set_time(self, time=0.0):
    399         """Set the model time (seconds)"""
     394        """Set the model time (seconds)."""
     395
    400396        # FIXME: this is setting the relative time
    401397        # Note that get_time and set_time are now not symmetric
     
    407403    # @return The absolute model time (seconds).
    408404    def get_time(self):
    409         """Get the absolute model time (seconds)"""
     405        """Get the absolute model time (seconds)."""
    410406
    411407        return self.time + self.starttime
     
    415411    # @param beta The new beta value.
    416412    def set_beta(self, beta):
    417         """Set default beta for limiting"""
     413        """Set default beta for limiting."""
    418414
    419415        self.beta = beta
    420416        for name in self.quantities:
    421             #print 'setting beta for quantity ',name
    422417            Q = self.quantities[name]
    423418            Q.set_beta(beta)
     
    427422    # @return The beta value used for limiting.
    428423    def get_beta(self):
    429         """Get default beta for limiting"""
     424        """Get default beta for limiting."""
    430425
    431426        return self.beta
     
    436431    # @note If 'n' is not 1 or 2, raise exception.
    437432    def set_default_order(self, n):
    438         """Set default (spatial) order to either 1 or 2"""
     433        """Set default (spatial) order to either 1 or 2."""
    439434
    440435        msg = 'Default order must be either 1 or 2. I got %s' % n
     
    443438        self.default_order = n
    444439        self._order_ = self.default_order
    445 
    446         if self.default_order == 1:
    447             pass
    448             #self.set_timestepping_method('euler')
    449             #self.set_all_limiters(beta_euler)
    450 
    451         if self.default_order == 2:
    452             pass
    453             #self.set_timestepping_method('rk2')
    454             #self.set_all_limiters(beta_rk2)
    455440
    456441    ##
     
    681666    ##
    682667    # @brief Set quantities based on a regional tag.
    683     # @param args 
    684     # @param kwargs 
     668    # @param args
     669    # @param kwargs
    685670    def set_region(self, *args, **kwargs):
    686671        """Set quantities based on a regional tag.
     
    12281213        self.starttime = float(time)
    12291214
    1230 #--------------------------
     1215################################################################################
    12311216# Main components of evolve
    1232 #--------------------------
     1217################################################################################
    12331218
    12341219    ##
     
    12751260               'evolving system, '
    12761261               'e.g. using the method set_boundary.\n'
    1277                'This system has the boundary tags %s ') \
    1278                    % self.get_boundary_tags()
     1262               'This system has the boundary tags %s '
     1263                   % self.get_boundary_tags())
    12791264        assert hasattr(self, 'boundary_objects'), msg
    12801265
     
    12871272
    12881273        if finaltime is not None and duration is not None:
    1289             # print 'F', finaltime, duration
    12901274            msg = 'Only one of finaltime and duration may be specified'
    12911275            raise Exception, msg
     
    12961280                self.finaltime = self.starttime + float(duration)
    12971281
    1298         N = len(self) # Number of triangles
    1299         self.yieldtime = 0.0 # Track time between 'yields'
     1282        N = len(self)            # Number of triangles
     1283        self.yieldtime = 0.0     # Track time between 'yields'
    13001284
    13011285        # Initialise interval of timestep sizes (for reporting only)
     
    13221306
    13231307        if skip_initial_step is False:
    1324             yield(self.time)  # Yield initial values
     1308            yield(self.time)      # Yield initial values
    13251309
    13261310        while True:
     
    13471331                if self.time > finaltime:
    13481332                    # FIXME (Ole, 30 April 2006): Do we need this check?
    1349                     # Probably not (Ole, 18 September 2008). Now changed to
    1350                     # Exception
     1333                    # Probably not (Ole, 18 September 2008).
     1334                    # Now changed to Exception.
    13511335                    msg = ('WARNING (domain.py): time overshot finaltime. '
    13521336                           'Contact Ole.Nielsen@ga.gov.au')
     
    14171401        self.backup_conserved_quantities()
    14181402
    1419         #--------------------------------------
     1403        ######
    14201404        # First euler step
    1421         #--------------------------------------
     1405        ######
    14221406
    14231407        # Compute fluxes across each element edge
     
    14421426        self.update_boundary()
    14431427
    1444         #------------------------------------
     1428        ######
    14451429        # Second Euler step
    1446         #------------------------------------
     1430        ######
    14471431
    14481432        # Compute fluxes across each element edge
     
    14521436        self.update_conserved_quantities()
    14531437
    1454         #------------------------------------
     1438        ######
    14551439        # Combine initial and final values
    14561440        # of conserved quantities and cleanup
    1457         #------------------------------------
     1441        ######
    14581442
    14591443        # Combine steps
     
    14841468        initial_time = self.time
    14851469
    1486         #--------------------------------------
     1470        ######
    14871471        # First euler step
    1488         #--------------------------------------
     1472        ######
    14891473
    14901474        # Compute fluxes across each element edge
     
    15091493        self.update_boundary()
    15101494
    1511         #------------------------------------
     1495        ######
    15121496        # Second Euler step
    1513         #------------------------------------
     1497        ######
    15141498
    15151499        # Compute fluxes across each element edge
     
    15191503        self.update_conserved_quantities()
    15201504
    1521         #------------------------------------
    1522         #Combine steps to obtain intermediate
    1523         #solution at time t^n + 0.5 h
    1524         #------------------------------------
     1505        ######
     1506        # Combine steps to obtain intermediate
     1507        # solution at time t^n + 0.5 h
     1508        ######
    15251509
    15261510        # Combine steps
     
    15391523        self.update_boundary()
    15401524
    1541         #------------------------------------
     1525        ######
    15421526        # Third Euler step
    1543         #------------------------------------
     1527        ######
    15441528
    15451529        # Compute fluxes across each element edge
     
    15491533        self.update_conserved_quantities()
    15501534
    1551         #------------------------------------
     1535        ######
    15521536        # Combine final and initial values
    15531537        # and cleanup
    1554         #------------------------------------
     1538        ######
    15551539
    15561540        # Combine steps
     
    16311615
    16321616    ##
    1633     # @brief 
    1634     # @param yieldstep 
    1635     # @param finaltime 
     1617    # @brief
     1618    # @param yieldstep
     1619    # @param finaltime
    16361620    def update_timestep(self, yieldstep, finaltime):
    16371621        from anuga.config import min_timestep, max_timestep
     
    17791763            elif self._order_ == 2:
    17801764                Q.extrapolate_second_order()
    1781                 #Q.limit()
    17821765            else:
    17831766                raise Exception, 'Unknown order'
    1784             #Q.interpolate_from_vertices_to_edges()
    17851767
    17861768    ##
    17871769    # @brief Calculate the norm of the centroid values of a specific quantity,
    17881770    #        using normfunc.
    1789     # @param quantity 
    1790     # @param normfunc 
     1771    # @param quantity
     1772    # @param normfunc
    17911773    def centroid_norm(self, quantity, normfunc):
    1792         """Calculate the norm of the centroid values
    1793         of a specific quantity, using normfunc.
     1774        """Calculate the norm of the centroid values of a specific quantity,
     1775        using normfunc.
    17941776
    17951777        normfunc should take a list to a float.
     
    18011783
    18021784
    1803 #------------------
     1785######
    18041786# Initialise module
    1805 #------------------
     1787######
    18061788
    18071789# Optimisation with psyco
    18081790from anuga.config import use_psyco
     1791
    18091792if use_psyco:
    18101793    try:
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r6221 r6226  
    1717from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
    1818from anuga.utilities.polygon import inside_polygon
    19 
    2019from anuga.geospatial_data.geospatial_data import Geospatial_data
    2120from anuga.fit_interpolate.fit import fit_to_mesh
    2221from anuga.config import points_file_block_line_size as default_block_line_size
    2322from anuga.config import epsilon
    24 from anuga.caching import cache
    25 
    26 
    2723
    2824import Numeric as num
    2925
    3026
     27##
     28# @brief Implement values at each triangular element.
    3129class Quantity:
    3230
    33     def __init__(self, domain, vertex_values=None):
    34    
     31    ##
     32    # @brief Construct values art each triangular element.
     33    # @param domain ??
     34    # @param vertex_values ??
     35    def __init__(self, domain,
     36                       vertex_values=None):
    3537        from anuga.abstract_2d_finite_volumes.domain import Domain
     38
    3639        msg = 'First argument in Quantity.__init__ '
    3740        msg += 'must be of class Domain (or a subclass thereof). '
    38         msg += 'I got %s.' %str(domain.__class__)
     41        msg += 'I got %s.' % str(domain.__class__)
    3942        assert isinstance(domain, Domain), msg
    4043
    4144        if vertex_values is None:
    42             N = len(domain) # number_of_elements
     45            N = len(domain)             # number_of_elements
    4346            self.vertex_values = num.zeros((N, 3), num.Float)
    4447        else:
     
    4649
    4750            N, V = self.vertex_values.shape
    48             assert V == 3,\
    49                    'Three vertex values per element must be specified'
    50 
    51 
    52             msg = 'Number of vertex values (%d) must be consistent with'\
    53                   %N
    54             msg += 'number of elements in specified domain (%d).'\
    55                    %len(domain)
    56 
     51            assert V == 3, 'Three vertex values per element must be specified'
     52
     53            msg = 'Number of vertex values (%d) must be consistent with' % N
     54            msg += 'number of elements in specified domain (%d).' % len(domain)
    5755            assert N == len(domain), msg
    5856
     
    6866
    6967        # Allocate space for Limiter Phi
    70         self.phi = num.zeros(N, num.Float)       
     68        self.phi = num.zeros(N, num.Float)
    7169
    7270        # Intialise centroid and edge_values
     
    8785        self.set_beta(1.0)
    8886
    89 
    90 
     87    #
    9188    # Methods for operator overloading
     89    #
     90
     91    ##
     92    # @brief Get length of object.
     93    # @return Return axis length of array.
    9294    def __len__(self):
    9395        return self.centroid_values.shape[0]
    9496
    95 
     97    ##
     98    # @brief Negate all values in this quantity.
    9699    def __neg__(self):
    97100        """Negate all values in this quantity giving meaning to the
     
    103106        return Q
    104107
    105 
     108    ##
     109    # @brief Add two quantities.
     110    # @param other The other quantity (to the right).
    106111    def __add__(self, other):
    107112        """Add to self anything that could populate a quantity
     
    119124        return result
    120125
     126    ##
     127    # @brief Add two quantities.
     128    # @param other The other quantity (to the left).
    121129    def __radd__(self, other):
    122130        """Handle cases like 7+Q, where Q is an instance of class Quantity
    123131        """
     132
    124133        return self + other
    125134
    126 
     135    ##
     136    # @brief Subtract two quantities.
     137    # @param other The other quantity (to the right).
    127138    def __sub__(self, other):
    128         return self + -other  #Invoke __neg__
    129 
     139        return self + -other            # Invoke self.__neg__()
     140
     141    ##
     142    # @brief Multiply two quantities.
     143    # @param other The other quantity (to the right).
    130144    def __mul__(self, other):
    131145        """Multiply self with anything that could populate a quantity
     
    138152        if isinstance(other, Quantity):
    139153            Q = other
    140         else:   
     154        else:
    141155            Q = Quantity(self.domain)
    142156            Q.set_values(other)
     
    151165        result.edge_values = self.edge_values * Q.edge_values
    152166        result.centroid_values = self.centroid_values * Q.centroid_values
    153        
     167
    154168        return result
    155169
     170    ##
     171    # @brief Multiply two quantities.
     172    # @param other The other quantity (to the left).
    156173    def __rmul__(self, other):
    157174        """Handle cases like 3*Q, where Q is an instance of class Quantity
    158175        """
     176
    159177        return self * other
    160178
     179    ##
     180    # @brief Divide two quantities.
     181    # @param other The other quantity (to the right).
    161182    def __div__(self, other):
    162183        """Divide self with anything that could populate a quantity
     
    172193        if isinstance(other, Quantity):
    173194            Q = other
    174         else:   
     195        else:
    175196            Q = Quantity(self.domain)
    176197            Q.set_values(other)
     
    188209        return result
    189210
     211    ##
     212    # @brief Divide two quantities.
     213    # @param other The other quantity (to the left).
    190214    def __rdiv__(self, other):
    191215        """Handle cases like 3/Q, where Q is an instance of class Quantity
    192216        """
     217
    193218        return self / other
    194219
     220    ##
     221    # @brief Raise a quaintity to some power.
     222    # @param other The power to raise quantity.
    195223    def __pow__(self, other):
    196224        """Raise quantity to (numerical) power
     
    201229        Example using __pow__:
    202230          Q = (Q1**2 + Q2**2)**0.5
    203 
    204231        """
    205232
    206233        if isinstance(other, Quantity):
    207234            Q = other
    208         else:   
     235        else:
    209236            Q = Quantity(self.domain)
    210237            Q.set_values(other)
     
    222249        return result
    223250
    224     #def __sqrt__(self, other):
    225     #    """Define in terms of x**0.5
    226     #    """
    227     #    pass
    228 
    229     def set_beta(self,beta):
    230         """Set default beta value for limiting
    231         """
     251    ##
     252    # @brief Set default beta value for limiting.
     253    # @param beta ??
     254    def set_beta(self, beta):
     255        """Set default beta value for limiting """
    232256
    233257        if beta < 0.0:
     
    235259        if beta > 2.0:
    236260            print 'WARNING: setting beta > 2.0'
    237            
     261
    238262        self.beta = beta
    239263
     264    ##
     265    # @brief Get the current beta value.
     266    # @return The current beta value.
    240267    def get_beta(self):
    241         """Get default beta value for limiting
    242         """
     268        """Get default beta value for limiting"""
    243269
    244270        return self.beta
    245271
     272    ##
     273    # @brief Compute interpolated values at edges and centroid.
     274    # @note vertex_values must be set before calling this.
    246275    def interpolate(self):
    247276        """Compute interpolated values at edges and centroid
    248277        Pre-condition: vertex_values have been set
    249278        """
    250        
     279
    251280        # FIXME (Ole): Maybe this function
    252281        # should move to the C-interface?
    253282        # However, it isn't called by validate_all.py, so it
    254283        # may not be that important to optimise it?
    255        
     284
    256285        N = self.vertex_values.shape[0]
    257286        for i in range(N):
     
    264293        self.interpolate_from_vertices_to_edges()
    265294
    266 
     295    ##
     296    # @brief ??
    267297    def interpolate_from_vertices_to_edges(self):
    268         # Call correct module function
    269         # (either from this module or C-extension)
     298        # Call correct module function (either from this module or C-extension)
    270299        interpolate_from_vertices_to_edges(self)
    271300
     301    ##
     302    # @brief ??
    272303    def interpolate_from_edges_to_vertices(self):
    273         # Call correct module function
    274         # (either from this module or C-extension)
     304        # Call correct module function (either from this module or C-extension)
    275305        interpolate_from_edges_to_vertices(self)
    276 
    277 
    278 
    279306
    280307    #---------------------------------------------
    281308    # Public interface for setting quantity values
    282309    #---------------------------------------------
    283     def set_values(self,
    284                    numeric=None,    # List, numeric array or constant
    285                    quantity=None,   # Another quantity
    286                    function=None,   # Callable object: f(x,y)
    287                    geospatial_data=None, # Arbitrary dataset
    288                    filename=None, attribute_name=None, # Input from file
    289                    alpha=None,
    290                    location='vertices',
    291                    polygon=None,
    292                    indices=None,
    293                    smooth=False,
    294                    verbose=False,
    295                    use_cache=False):
    296 
     310
     311    ##
     312    # @brief Set values for quantity based on different sources.
     313    # @param numeric A num array, list or constant value.
     314    # @param quantity Another Quantity.
     315    # @param function Any callable object that takes two 1d arrays.
     316    # @param geospatial_data Arbitrary instance of class Geospatial_data
     317    # @param filename Path to a points file.
     318    # @param attribute_name If specified any array using that name will be used.
     319    # @param alpha Smoothing parameter to be used with fit_interpolate.fit.
     320    # @param location Where to store values (vertices, edges, centroids).
     321    # @param polygon Restrict update to locations that fall inside polygon.
     322    # @param indices Restrict update to locations specified by this.
     323    # @param smooth If True, smooth vertex values.
     324    # @param verbose True if this method is to be verbose.
     325    # @param use_cache If True cache results for fit_interpolate.fit.
     326    # @note Exactly one of 'numeric', 'quantity', 'function', 'filename'
     327    #       must be present.
     328    def set_values(self, numeric=None,         # List, numeric array or constant
     329                         quantity=None,        # Another quantity
     330                         function=None,        # Callable object: f(x,y)
     331                         geospatial_data=None, # Arbitrary dataset
     332                         filename=None,
     333                         attribute_name=None,  # Input from file
     334                         alpha=None,
     335                         location='vertices',
     336                         polygon=None,
     337                         indices=None,
     338                         smooth=False,
     339                         verbose=False,
     340                         use_cache=False):
    297341        """Set values for quantity based on different sources.
    298342
     
    366410                 constant values so far.
    367411
    368         indices: Restrict update of quantity to locations that are 
     412        indices: Restrict update of quantity to locations that are
    369413                 identified by indices (e.g. node ids if location
    370414                 is 'unique vertices' or triangle ids otherwise).
    371        
     415
    372416        verbose: True means that output to stdout is generated
    373417
     
    391435        # perhaps the notion of location and indices simplified
    392436
    393         # FIXME (Ole): Need to compute indices based on polygon 
     437        # FIXME (Ole): Need to compute indices based on polygon
    394438        # (and location) and use existing code after that.
    395        
     439
    396440        # See ticket:275, ticket:250, ticeket:254 for refactoring plan
    397        
     441
    398442        if polygon is not None:
    399443            if indices is not None:
     
    405449            msg += 'a constant.'
    406450            if numeric is None:
    407                 raise Exception, msg           
     451                raise Exception, msg
    408452            else:
    409453                # Check that numeric is as constant
    410454                assert type(numeric) in [FloatType, IntType, LongType], msg
    411455
    412 
    413456            location = 'centroids'
    414 
    415457
    416458            points = self.domain.get_centroid_coordinates(absolute=True)
    417459            indices = inside_polygon(points, polygon)
    418            
    419             self.set_values_from_constant(numeric,
    420                                           location, indices, verbose)
    421 
     460
     461            self.set_values_from_constant(numeric, location, indices, verbose)
    422462
    423463            self.extrapolate_first_order()
    424464
    425465            if smooth:
    426                 self.smooth_vertex_values(use_cache=use_cache,
    427                                           verbose=verbose)
    428 
    429                
     466                self.smooth_vertex_values()
     467
    430468            return
    431        
    432        
    433 
    434 
    435 
    436469
    437470        # General input checks
    438471        L = [numeric, quantity, function, geospatial_data, filename]
    439         msg = 'Exactly one of the arguments '+\
    440               'numeric, quantity, function, geospatial_data, '+\
    441               'or filename must be present.'
     472        msg = ('Exactly one of the arguments numeric, quantity, function, '
     473               'geospatial_data, or filename must be present.')
    442474        assert L.count(None) == len(L)-1, msg
    443 
    444475
    445476        if location == 'edges':
    446477            msg = 'edges has been deprecated as valid location'
    447478            raise Exception, msg
    448            
     479
    449480        if location not in ['vertices', 'centroids', 'unique vertices']:
    450             msg = 'Invalid location: %s' %location
     481            msg = 'Invalid location: %s' % location
    451482            raise Exception, msg
    452 
    453483
    454484        msg = 'Indices must be a list or None'
    455485        assert type(indices) in [ListType, NoneType, num.ArrayType], msg
    456486
    457 
    458 
    459487        # Determine which 'set_values_from_...' to use
    460 
    461488        if numeric is not None:
    462489            if type(numeric) in [FloatType, IntType, LongType]:
    463                 self.set_values_from_constant(numeric,
    464                                               location, indices, verbose)
     490                self.set_values_from_constant(numeric, location,
     491                                              indices, verbose)
    465492            elif type(numeric) in [num.ArrayType, ListType]:
    466                 self.set_values_from_array(numeric,
    467                                            location, indices,
    468                                            use_cache=use_cache,
    469                                            verbose=verbose)
     493                self.set_values_from_array(numeric, location, indices, verbose)
    470494            elif callable(numeric):
    471                 self.set_values_from_function(numeric,
    472                                               location, indices,
    473                                               use_cache=use_cache,
    474                                               verbose=verbose)
     495                self.set_values_from_function(numeric, location,
     496                                              indices, verbose)
    475497            elif isinstance(numeric, Quantity):
    476                 self.set_values_from_quantity(numeric,
    477                                               location, indices,
    478                                               verbose=verbose)
     498                self.set_values_from_quantity(numeric, location,
     499                                              indices, verbose)
    479500            elif isinstance(numeric, Geospatial_data):
    480                 self.set_values_from_geospatial_data(numeric,
    481                                                      alpha,
     501                self.set_values_from_geospatial_data(numeric, alpha, location,
     502                                                     indices, verbose=verbose,
     503                                                     use_cache=use_cache)
     504            else:
     505                msg = 'Illegal type for argument numeric: %s' % str(numeric)
     506                raise msg
     507        elif quantity is not None:
     508            self.set_values_from_quantity(quantity, location, indices, verbose)
     509        elif function is not None:
     510            msg = 'Argument function must be callable'
     511            assert callable(function), msg
     512            self.set_values_from_function(function, location, indices, verbose)
     513        elif geospatial_data is not None:
     514                self.set_values_from_geospatial_data(geospatial_data, alpha,
    482515                                                     location, indices,
    483516                                                     verbose=verbose,
    484517                                                     use_cache=use_cache)
    485             else:
    486                 msg = 'Illegal type for argument numeric: %s' %str(numeric)
    487                 raise msg
    488 
    489         elif quantity is not None:
    490             self.set_values_from_quantity(quantity,
    491                                           location, indices, verbose)
    492         elif function is not None:
    493             msg = 'Argument function must be callable'
    494             assert callable(function), msg
    495             self.set_values_from_function(function,
    496                                           location,
    497                                           indices,
    498                                           use_cache=use_cache,
    499                                           verbose=verbose)
    500         elif geospatial_data is not None:
    501                 self.set_values_from_geospatial_data(geospatial_data,
    502                                                      alpha,
    503                                                      location, indices,
    504                                                      verbose=verbose,
    505                                                      use_cache=use_cache)
    506            
    507518        elif filename is not None:
    508519            if hasattr(self.domain, 'points_file_block_line_size'):
     
    510521            else:
    511522                max_read_lines = default_block_line_size
    512             self.set_values_from_file(filename, attribute_name, alpha,
    513                                       location, indices,
    514                                       verbose=verbose,
     523            self.set_values_from_file(filename, attribute_name, alpha, location,
     524                                      indices, verbose=verbose,
    515525                                      max_read_lines=max_read_lines,
    516526                                      use_cache=use_cache)
    517527        else:
    518             raise Exception, 'This can\'t happen :-)'
    519 
    520 
     528            raise Exception, "This can't happen :-)"
    521529
    522530        # Update all locations in triangles
     
    531539
    532540
    533     #-------------------------------------------------------------       
     541    #---------------------------------------------------------------------------
    534542    # Specific internal functions for setting values based on type
    535     #-------------------------------------------------------------           
    536    
    537     def set_values_from_constant(self, X,
    538                                  location, indices, verbose):
    539         """Set quantity values from specified constant X
    540         """
     543    #---------------------------------------------------------------------------
     544
     545    ##
     546    # @brief Set quantity values from specified constant.
     547    # @param X The constant to set quantity values to.
     548    # @param location
     549    # @param indices
     550    # @param verbose
     551    def set_values_from_constant(self, X, location, indices, verbose):
     552        """Set quantity values from specified constant X"""
    541553
    542554        # FIXME (Ole): Somehow indices refer to centroids
    543555        # rather than vertices as default. See unit test
    544556        # test_set_vertex_values_using_general_interface_with_subset(self):
    545        
    546557
    547558        if location == 'centroids':
     
    552563                for i in indices:
    553564                    self.centroid_values[i] = X
    554 
    555         #elif location == 'edges':
    556         #    if indices is None:
    557         #        self.edge_values[:] = X
    558         #    else:
    559         #        # Brute force
    560         #        for i in indices:
    561         #            self.edge_values[i] = X
    562 
    563565        elif location == 'unique vertices':
    564566            if indices is None:
    565567                self.edge_values[:] = X  #FIXME (Ole): Shouldn't this be vertex_values?
    566568            else:
    567 
    568569                # Go through list of unique vertices
    569570                for unique_vert_id in indices:
    570 
    571                     triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
    572                    
     571                    triangles = \
     572                        self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
     573
    573574                    # In case there are unused points
    574575                    if len(triangles) == 0:
    575576                        continue
    576                    
     577
    577578                    # Go through all triangle, vertex pairs
    578579                    # and set corresponding vertex value
     
    590591                    self.vertex_values[i_vertex] = X
    591592
    592 
    593 
    594 
     593    ##
     594    # @brief Set values for a quantity.
     595    # @param values Array of values.
     596    # @param location Where values are to be stored.
     597    # @param indices Limit update to these indices.
     598    # @param verbose True if this method is to be verbose.
    595599    def set_values_from_array(self, values,
    596                               location='vertices',
    597                               indices=None,
    598                               use_cache=False,
    599                               verbose=False):
     600                                    location='vertices',
     601                                    indices=None,
     602                                    verbose=False):
    600603        """Set values for quantity
    601604
     
    628631        if indices is not None:
    629632            indices = num.array(indices, num.Int)
    630             msg = 'Number of values must match number of indices:'
    631             msg += ' You specified %d values and %d indices'\
    632                    %(values.shape[0], indices.shape[0])
     633            msg = ('Number of values must match number of indices: You '
     634                   'specified %d values and %d indices'
     635                   % (values.shape[0], indices.shape[0]))
    633636            assert values.shape[0] == indices.shape[0], msg
    634637
     
    650653                for i in range(len(indices)):
    651654                    self.centroid_values[indices[i]] = values[i]
    652 
    653655        elif location == 'unique vertices':
    654             assert len(values.shape) == 1 or num.allclose(values.shape[1:], 1),\
    655                    'Values array must be 1d'
    656 
    657             self.set_vertex_values(values.flat,
    658                                    indices=indices,
    659                                    use_cache=use_cache,
    660                                    verbose=verbose)
    661            
     656            assert (len(values.shape) == 1 or num.allclose(values.shape[1:], 1),
     657                    'Values array must be 1d')
     658
     659            self.set_vertex_values(values.flat, indices=indices)
    662660        else:
    663661            # Location vertices
    664662            if len(values.shape) == 1:
    665                 # This is the common case arising from fitted
    666                 # values (e.g. from pts file).
    667                 self.set_vertex_values(values,
    668                                        indices=indices,
    669                                        use_cache=use_cache,
    670                                        verbose=verbose)
    671 
     663                self.set_vertex_values(values, indices=indices)
    672664            elif len(values.shape) == 2:
    673665                # Vertex values are given as a triplet for each triangle
    674 
    675666                msg = 'Array must be N x 3'
    676667                assert values.shape[1] == 3, msg
     
    683674            else:
    684675                msg = 'Values array must be 1d or 2d'
    685                 raise msg
    686            
    687 
    688     def set_values_from_quantity(self, q,
    689                                  location, indices, verbose):
     676                raise Exception, msg
     677
     678    ##
     679    # @brief Set quantity values from a specified quantity instance.
     680    # @param q The quantity instance to take values from.
     681    # @param location IGNORED, 'vertices' ALWAYS USED!
     682    # @param indices ??
     683    # @param verbose True if this method is to be verbose.
     684    def set_values_from_quantity(self, q, location, indices, verbose):
    690685        """Set quantity values from specified quantity instance q
    691686
     
    702697
    703698        self.set_values(A, location='vertices',
    704                         indices=indices,
    705                         verbose=verbose)
    706 
    707 
     699                        indices=indices, verbose=verbose)
     700
     701    ##
     702    # @brief Set quantity values from a specified quantity instance.
     703    # @param f Callable that takes two 1d array -> 1d array.
     704    # @param location Where values are to be stored.
     705    # @param indices ??
     706    # @param verbose True if this method is to be verbose.
    708707    def set_values_from_function(self, f,
    709                                  location='vertices',
    710                                  indices=None,
    711                                  use_cache=False,
    712                                  verbose=False):
     708                                       location='vertices',
     709                                       indices=None,
     710                                       verbose=False):
    713711        """Set values for quantity using specified function
    714712
    715713        Input
    716        
    717714        f: x, y -> z Function where x, y and z are arrays
    718715        location: Where values are to be stored.
     
    720717                  unique vertices
    721718                  Default is "vertices"
    722         indices: 
    723 
    724                  
     719        indices:
    725720        """
    726721
    727722        # FIXME: Should check that function returns something sensible and
    728         # raise a meaningfull exception if it returns None for example
     723        # raise a meaningful exception if it returns None for example
    729724
    730725        # FIXME: Should supply absolute coordinates
    731 
    732726
    733727        # Compute the function values and call set_values again
     
    735729            if indices is None:
    736730                indices = range(len(self))
    737                
     731
    738732            V = num.take(self.domain.get_centroid_coordinates(), indices)
    739 
    740             x = V[:,0]; y = V[:,1]
    741             if use_cache is True:
    742                 res = cache(f, (x, y),
    743                             verbose=verbose)
    744             else:
    745                 res = f(x, y)
    746 
    747             self.set_values(res,
    748                             location=location,
    749                             indices=indices)
    750            
     733            self.set_values(f(V[:,0], V[:,1]),
     734                            location=location, indices=indices)
    751735        elif location == 'vertices':
    752             # This is the default branch taken by set_quantity
    753            
    754736            M = self.domain.number_of_triangles
    755737            V = self.domain.get_vertex_coordinates()
    756738
    757             x = V[:,0]; y = V[:,1]
    758             if use_cache is True:
    759                 #print 'Caching function'
    760                 values = cache(f, (x, y),
    761                                verbose=verbose)               
    762             else:
    763                 if verbose is True:
    764                     print 'Evaluating function in set_values'
    765                 values = f(x, y)
    766 
    767             #print 'value', min(x), max(x), min(y), max(y), max(values), len(values)
    768                
     739            x = V[:,0];
     740            y = V[:,1];
     741            values = f(x, y)
    769742
    770743            # FIXME (Ole): This code should replace all the
     
    773746            # If that could be resolved this one will be
    774747            # more robust and simple.
    775            
    776             #values = num.reshape(values, (M,3))
    777             #self.set_values(values,
    778             #                location='vertices',
    779             #                indices=indices)
    780 
    781748
    782749            # This should be removed
    783750            if is_scalar(values):
    784751                # Function returned a constant value
    785                 self.set_values_from_constant(values,
    786                                               location, indices, verbose)
     752                self.set_values_from_constant(values, location,
     753                                              indices, verbose)
    787754                return
    788755
    789             # This should be removed           
     756            # This should be removed
    790757            if indices is None:
    791758                for j in range(3):
    792                     self.vertex_values[:,j] = values[j::3]                 
    793             else:   
     759                    self.vertex_values[:, j] = values[j::3]
     760            else:
    794761                # Brute force
    795762                for i in indices:
    796763                    for j in range(3):
    797                         self.vertex_values[i,j] = values[3*i+j]
    798 
    799 
     764                        self.vertex_values[i, j] = values[3*i + j]
    800765        else:
    801             raise 'Not implemented: %s' %location
    802 
    803 
    804 
    805     def set_values_from_geospatial_data(self, geospatial_data, alpha,
    806                                         location, indices,
    807                                         verbose=False,
    808                                         use_cache=False):
    809         """ Set values based on geo referenced geospatial data object.
    810         """
     766            raise Exception, 'Not implemented: %s' % location
     767
     768    ##
     769    # @brief Set values based on geo referenced geospatial data object.
     770    # @param geospatial_data ??
     771    # @param alpha ??
     772    # @param location ??
     773    # @param indices ??
     774    # @param verbose ??
     775    # @param use_cache ??
     776    def set_values_from_geospatial_data(self, geospatial_data,
     777                                              alpha,
     778                                              location,
     779                                              indices,
     780                                              verbose=False,
     781                                              use_cache=False):
     782        """Set values based on geo referenced geospatial data object."""
     783
     784        from anuga.coordinate_transforms.geo_reference import Geo_reference
    811785
    812786        points = geospatial_data.get_data_points(absolute=False)
     
    814788        data_georef = geospatial_data.get_geo_reference()
    815789
    816 
    817         from anuga.coordinate_transforms.geo_reference import Geo_reference
    818 
    819 
    820790        points = ensure_numeric(points, num.Float)
    821791        values = ensure_numeric(values, num.Float)
    822792
    823793        if location != 'vertices':
    824             msg = 'set_values_from_points is only defined for '+\
    825                   'location=\'vertices\''
    826             raise ms
    827 
    828         #coordinates = self.domain.get_nodes()
    829         #triangles = self.domain.get_triangles()
    830 
     794            msg = ("set_values_from_points is only defined for "
     795                   "location='vertices'")
     796            raise Exception, msg
    831797
    832798        # Take care of georeferencing
     
    834800            data_georef = Geo_reference()
    835801
    836 
    837802        mesh_georef = self.domain.geo_reference
    838803
    839 
    840804        # Call fit_interpolate.fit function
    841         # args = (coordinates, triangles, points, values)
    842805        args = (points, )
    843806        kwargs = {'vertex_coordinates': None,
     
    850813                  'verbose': verbose}
    851814
    852         vertex_attributes = apply(fit_to_mesh,
    853                                   args, kwargs)       
     815        vertex_attributes = apply(fit_to_mesh, args, kwargs)
    854816
    855817        # Call underlying method using array values
    856         self.set_values_from_array(vertex_attributes,
    857                                    location, indices,
    858                                    use_cache=use_cache,
    859                                    verbose=verbose)
    860 
    861 
    862 
    863     def set_values_from_points(self, points, values, alpha,
    864                                location, indices,
    865                                data_georef=None,
    866                                verbose=False,
    867                                use_cache=False):
    868         """
    869         Set quantity values from arbitray data points using
    870         fit_interpolate.fit
    871         """
     818        self.set_values_from_array(vertex_attributes, location,
     819                                   indices, verbose)
     820
     821    ##
     822    # @brief Set quantity values from arbitray data points.
     823    # @param points ??
     824    # @param values ??
     825    # @param alpha ??
     826    # @param location ??
     827    # @param indices ??
     828    # @param data_georef ??
     829    # @param verbose True if this method is to be verbose.
     830    # @param use_cache ??
     831    def set_values_from_points(self, points,
     832                                     values,
     833                                     alpha,
     834                                     location,
     835                                     indices,
     836                                     data_georef=None,
     837                                     verbose=False,
     838                                     use_cache=False):
     839        """Set quantity values from arbitray data points using fit_interpolate.fit"""
    872840
    873841        raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
    874        
    875 
    876     def set_values_from_file(self, filename, attribute_name, alpha,
    877                              location, indices,
    878                              verbose=False,
    879                              use_cache=False,
    880                              max_read_lines=None):
    881         """Set quantity based on arbitrary points in a points file
    882         using attribute_name selects name of attribute
    883         present in file.
     842
     843    ##
     844    # @brief Set quantity based on arbitrary points in a points file.
     845    # @param filename Path to the points file.
     846    # @param attribute_name
     847    # @param alpha
     848    # @param location
     849    # @param indices
     850    # @param verbose True if this method is to be verbose.
     851    # @param use_cache
     852    # @param max_read_lines
     853    def set_values_from_file(self, filename,
     854                                   attribute_name,
     855                                   alpha,
     856                                   location,
     857                                   indices,
     858                                   verbose=False,
     859                                   use_cache=False,
     860                                   max_read_lines=None):
     861        """Set quantity based on arbitrary points in a points file using
     862        attribute_name selects name of attribute present in file.
    884863        If attribute_name is not specified, use first available attribute
    885         as defined in geospatial_data.
     864        as defined in geospatial_data.
    886865        """
    887866
    888867        from types import StringType
     868
    889869        msg = 'Filename must be a text string'
    890870        assert type(filename) == StringType, msg
    891871
    892 
    893872        if location != 'vertices':
    894             msg = 'set_values_from_file is only defined for '+\
    895                   'location=\'vertices\''
    896             raise msg
    897 
    898         if True: 
     873            msg = "set_values_from_file is only defined for location='vertices'"
     874            raise Exception, msg
     875
     876        if True:
    899877            # Use mesh as defined by domain
    900             # This used to cause problems for caching
    901             # due to quantities changing, but
    902             # it now works using the appropriate Mesh object.
     878            # This used to cause problems for caching due to quantities
     879            # changing, but it now works using the appropriate Mesh object.
    903880            # This should close ticket:242
    904881            vertex_attributes = fit_to_mesh(filename,
    905                                             mesh=self.domain.mesh, 
     882                                            mesh=self.domain.mesh,
    906883                                            alpha=alpha,
    907884                                            attribute_name=attribute_name,
     
    911888        else:
    912889            # This variant will cause Mesh object to be recreated
    913             # in fit_to_mesh thus doubling up on the neighbour structure 
     890            # in fit_to_mesh thus doubling up on the neighbour structure
    914891            # FIXME(Ole): This is now obsolete 19 Jan 2009.
    915892            nodes = self.domain.get_nodes(absolute=True)
    916             triangles = self.domain.get_triangles()     
     893            triangles = self.domain.get_triangles()
    917894            vertex_attributes = fit_to_mesh(filename,
    918                                             nodes, triangles, 
     895                                            nodes, triangles,
    919896                                            mesh=None,
    920897                                            alpha=alpha,
     
    923900                                            verbose=verbose,
    924901                                            max_read_lines=max_read_lines)
    925                                            
     902
    926903        # Call underlying method using array values
    927         if verbose:
    928             print 'Applying fitted data to domain'
    929         self.set_values_from_array(vertex_attributes,
    930                                    location, indices,
    931                                    use_cache=use_cache,
    932                                    verbose=verbose)
    933 
    934    
    935    
    936     #-----------------------------------------------------   
     904        self.set_values_from_array(vertex_attributes, location,
     905                                   indices, verbose)
     906
     907    ##
     908    # @brief Get index for maximum or minimum value of quantity.
     909    # @param mode Either 'max' or 'min'.
     910    # @param indices Set of IDs of elements to work on.
    937911    def get_extremum_index(self, mode=None, indices=None):
    938912        """Return index for maximum or minimum value of quantity (on centroids)
     
    958932        if mode is None or mode == 'max':
    959933            i = num.argmax(V)
    960         elif mode == 'min':   
     934        elif mode == 'min':
    961935            i = num.argmin(V)
    962 
    963            
     936        else:
     937            raise ValueError, 'Bad mode value, got: %s' % str(mode)
     938
    964939        if indices is None:
    965940            return i
     
    967942            return indices[i]
    968943
    969 
     944    ##
     945    # @brief Get index for maximum value of quantity.
     946    # @param indices Set of IDs of elements to work on.
    970947    def get_maximum_index(self, indices=None):
    971         """See get extreme index for details
    972         """
    973 
    974         return self.get_extremum_index(mode='max',
    975                                        indices=indices)
    976 
    977 
    978        
     948        """See get extreme index for details"""
     949
     950        return self.get_extremum_index(mode='max', indices=indices)
     951
     952    ##
     953    # @brief Return maximum value of quantity (on centroids).
     954    # @param indices Set of IDs of elements to work on.
    979955    def get_maximum_value(self, indices=None):
    980956        """Return maximum value of quantity (on centroids)
     
    987963
    988964        Note, we do not seek the maximum at vertices as each vertex can
    989         have multiple values - one for each triangle sharing it           
    990         """
    991 
     965        have multiple values - one for each triangle sharing it
     966        """
    992967
    993968        i = self.get_maximum_index(indices)
    994         V = self.get_values(location='centroids') #, indices=indices)
    995        
     969        V = self.get_values(location='centroids')      #, indices=indices)
     970
    996971        return V[i]
    997        
    998 
     972
     973    ##
     974    # @brief Get location of maximum value of quantity (on centroids).
     975    # @param indices Set of IDs of elements to work on.
    999976    def get_maximum_location(self, indices=None):
    1000977        """Return location of maximum value of quantity (on centroids)
     
    1005982        Usage:
    1006983            x, y = get_maximum_location()
    1007 
    1008984
    1009985        Notes:
     
    1012988
    1013989            If there are multiple cells with same maximum value, the
    1014             first cell encountered in the triangle array is returned.       
     990            first cell encountered in the triangle array is returned.
    1015991        """
    1016992
     
    1020996        return x, y
    1021997
    1022 
     998    ##
     999    # @brief  Get index for minimum value of quantity.
     1000    # @param indices Set of IDs of elements to work on.
    10231001    def get_minimum_index(self, indices=None):
    1024         """See get extreme index for details
    1025         """       
    1026 
    1027         return self.get_extremum_index(mode='min',
    1028                                        indices=indices)
    1029 
    1030 
     1002        """See get extreme index for details"""
     1003
     1004        return self.get_extremum_index(mode='min', indices=indices)
     1005
     1006    ##
     1007    # @brief Return minimum value of quantity (on centroids).
     1008    # @param indices Set of IDs of elements to work on.
    10311009    def get_minimum_value(self, indices=None):
    10321010        """Return minimum value of quantity (on centroids)
     
    10381016            v = get_minimum_value()
    10391017
    1040         See get_maximum_value for more details.   
    1041         """
    1042 
     1018        See get_maximum_value for more details.
     1019        """
    10431020
    10441021        i = self.get_minimum_index(indices)
    10451022        V = self.get_values(location='centroids')
    1046        
     1023
    10471024        return V[i]
    1048        
    1049 
     1025
     1026
     1027    ##
     1028    # @brief Get location of minimum value of quantity (on centroids).
     1029    # @param indices Set of IDs of elements to work on.
    10501030    def get_minimum_location(self, indices=None):
    10511031        """Return location of minimum value of quantity (on centroids)
     
    10561036        Usage:
    10571037            x, y = get_minimum_location()
    1058 
    10591038
    10601039        Notes:
     
    10631042
    10641043            If there are multiple cells with same maximum value, the
    1065             first cell encountered in the triangle array is returned.       
     1044            first cell encountered in the triangle array is returned.
    10661045        """
    10671046
     
    10711050        return x, y
    10721051
    1073 
    1074 
     1052    ##
     1053    # @brief Get values at interpolation points.
     1054    # @param interpolation_points List of UTM coords or geospatial data object.
     1055    # @param use_cache ??
     1056    # @param verbose True if this method is to be verbose.
    10751057    def get_interpolated_values(self, interpolation_points,
    1076                                 use_cache=False,
    1077                                 verbose=False):
    1078         """ Get values at interpolation points
    1079        
    1080         The argument interpolation points must be given as either a 
     1058                                      use_cache=False,
     1059                                      verbose=False):
     1060        """Get values at interpolation points
     1061
     1062        The argument interpolation points must be given as either a
    10811063        list of absolute UTM coordinates or a geospatial data object.
    10821064        """
    1083        
    10841065
    10851066        # FIXME (Ole): Points might be converted to coordinates relative to mesh origin
    1086         # This could all be refactored using the 
    1087         # 'change_points_geo_ref' method of Class geo_reference. 
     1067        # This could all be refactored using the
     1068        # 'change_points_geo_ref' method of Class geo_reference.
    10881069        # The purpose is to make interpolation points relative
    10891070        # to the mesh origin.
    10901071        #
    10911072        # Speed is also a consideration here.
    1092        
    1093        
    1094         # Ensure that interpolation points is either a list of
     1073
     1074        # Ensure that interpolation points is either a list of
    10951075        # points, Nx2 array, or geospatial and convert to Numeric array
    1096         if isinstance(interpolation_points, Geospatial_data):       
     1076        if isinstance(interpolation_points, Geospatial_data):
    10971077            # Ensure interpolation points are in absolute UTM coordinates
    1098             interpolation_points = interpolation_points.get_data_points(absolute=True)
    1099                
     1078            interpolation_points = \
     1079                interpolation_points.get_data_points(absolute=True)
     1080
    11001081        # Reconcile interpolation points with georeference of domain
    1101         interpolation_points = self.domain.geo_reference.get_relative(interpolation_points)
     1082        interpolation_points = \
     1083            self.domain.geo_reference.get_relative(interpolation_points)
    11021084        interpolation_points = ensure_numeric(interpolation_points)
    11031085
    1104        
     1086
    11051087        # Get internal representation (disconnected) of vertex values
    11061088        vertex_values, triangles = self.get_vertex_values(xy=False,
    1107                                                           smooth=False)               
    1108    
     1089                                                          smooth=False)
     1090
    11091091        # Get possibly precomputed interpolation object
    11101092        I = self.domain.get_interpolation_object()
    11111093
    1112         # Call interpolate method with interpolation points               
     1094        # Call interpolate method with interpolation points
    11131095        result = I.interpolate_block(vertex_values, interpolation_points,
    1114                                      use_cache=use_cache,
    1115                                      verbose=verbose)
    1116                                
     1096                                     use_cache=use_cache, verbose=verbose)
     1097
    11171098        return result
    1118        
    1119        
    1120 
    1121 
    1122     def get_values(self,
    1123                    interpolation_points=None,
    1124                    location='vertices',
    1125                    indices=None,
    1126                    use_cache=False,
    1127                    verbose=False):
    1128         """get values for quantity
    1129 
    1130         Extract values for quantity as a Numeric array.       
    1131        
     1099
     1100    ##
     1101    # @brief Get values as an array.
     1102    # @param interpolation_points List of coords to get values at.
     1103    # @param location Where to store results.
     1104    # @param indices Set of IDs of elements to work on.
     1105    # @param use_cache
     1106    # @param verbose True if this method is to be verbose.
     1107    def get_values(self, interpolation_points=None,
     1108                         location='vertices',
     1109                         indices=None,
     1110                         use_cache=False,
     1111                         verbose=False):
     1112        """Get values for quantity
     1113
     1114        Extract values for quantity as a Numeric array.
     1115
    11321116        Inputs:
    11331117           interpolation_points: List of x, y coordinates where value is
    1134                                  sought (using interpolation). If points 
    1135                                  are given, values of location and indices 
     1118                                 sought (using interpolation). If points
     1119                                 are given, values of location and indices
    11361120                                 are ignored. Assume either absolute UTM
    11371121                                 coordinates or geospatial data object.
    1138        
     1122
    11391123           location: Where values are to be stored.
    11401124                     Permissible options are: vertices, edges, centroids
     
    11481132        values will be a list or a Numerical array of length N, N being
    11491133        the number of elements.
    1150        
     1134
    11511135        In case of location == 'vertices' or 'edges' the dimension of
    11521136        returned values will be of dimension Nx3
     
    11541138        In case of location == 'unique vertices' the average value at
    11551139        each vertex will be returned and the dimension of returned values
    1156         will be a 1d array of length "number of vertices" 
    1157        
     1140        will be a 1d array of length "number of vertices"
     1141
    11581142        Indices is the set of element ids that the operation applies to.
    11591143
     
    11611145        internal ordering.
    11621146        """
    1163        
    11641147
    11651148        # FIXME (Ole): I reckon we should have the option of passing a
    11661149        #              polygon into get_values. The question becomes how
    11671150        #              resulting values should be ordered.
    1168        
     1151
    11691152        if verbose is True:
    1170             print 'Getting values from %s' %location
     1153            print 'Getting values from %s' % location
    11711154
    11721155        if interpolation_points is not None:
     
    11741157                                                use_cache=use_cache,
    11751158                                                verbose=verbose)
    1176        
    1177        
     1159
    11781160        # FIXME (Ole): Consider deprecating 'edges' - but not if it is used
    1179         # elsewhere in ANUGA. 
     1161        # elsewhere in ANUGA.
    11801162        # Edges have already been deprecated in set_values, see changeset:5521,
    11811163        # but *might* be useful in get_values. Any thoughts anyone?
    1182        
    1183         if location not in ['vertices', 'centroids', 'edges',
    1184                             'unique vertices']:
    1185             msg = 'Invalid location: %s' %location
     1164
     1165        if location not in ['vertices', 'centroids',
     1166                            'edges', 'unique vertices']:
     1167            msg = 'Invalid location: %s' % location
    11861168            raise msg
    11871169
    11881170        import types
     1171
    11891172        assert type(indices) in [types.ListType, types.NoneType,
    1190                                  num.ArrayType],\
    1191                                  'Indices must be a list or None'
     1173                                 num.ArrayType], \
     1174                   'Indices must be a list or None'
    11921175
    11931176        if location == 'centroids':
    11941177            if (indices ==  None):
    11951178                indices = range(len(self))
    1196             return num.take(self.centroid_values,indices)
     1179            return num.take(self.centroid_values, indices)
    11971180        elif location == 'edges':
    11981181            if (indices ==  None):
    11991182                indices = range(len(self))
    1200             return num.take(self.edge_values,indices)
     1183            return num.take(self.edge_values, indices)
    12011184        elif location == 'unique vertices':
    12021185            if (indices ==  None):
     
    12071190            for unique_vert_id in indices:
    12081191                triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
    1209                    
     1192
    12101193                # In case there are unused points
    12111194                if len(triangles) == 0:
    12121195                    msg = 'Unique vertex not associated with triangles'
    1213                     raise msg
     1196                    raise Exception, msg
    12141197
    12151198                # Go through all triangle, vertex pairs
    12161199                # Average the values
    1217                
    12181200                # FIXME (Ole): Should we merge this with get_vertex_values
    12191201                sum = 0
    12201202                for triangle_id, vertex_id in triangles:
    12211203                    sum += self.vertex_values[triangle_id, vertex_id]
    1222                 vert_values.append(sum/len(triangles))
     1204                vert_values.append(sum / len(triangles))
    12231205            return num.array(vert_values, num.Float)
    12241206        else:
     
    12271209            return num.take(self.vertex_values, indices)
    12281210
    1229 
    1230 
    1231     def set_vertex_values(self, A,
    1232                           indices=None,
    1233                           use_cache=False,
    1234                           verbose=False):
     1211    ##
     1212    # @brief Set vertex values for all unique vertices based on array.
     1213    # @param A Array to set values with.
     1214    # @param indices Set of IDs of elements to work on.
     1215    def set_vertex_values(self, A, indices=None):
    12351216        """Set vertex values for all unique vertices based on input array A
    1236         which has one entry per unique vertex, i.e.
    1237         one value for each row in array self.domain.nodes.
     1217        which has one entry per unique vertex, i.e. one value for each row in
     1218        array self.domain.nodes.
    12381219
    12391220        indices is the list of vertex_id's that will be set.
     
    12421223        """
    12431224
    1244 
    1245         # Assert that A can be converted to a Numeric array of appropriate dim
     1225        # Check that A can be converted to array and is of appropriate dim
    12461226        A = ensure_numeric(A, num.Float)
    1247 
    1248         # print 'SHAPE A', A.shape
    12491227        assert len(A.shape) == 1
    12501228
     
    12561234            vertex_list = indices
    12571235
    1258 
    1259         #FIXME(Ole): This function ought to be faster.
    1260         # We need to get the triangles_and_vertices list
    1261         # from domain in one hit, then cache the computation of the
    1262         # Nx3 array of vertex values that can then be assigned using
    1263         # set_values_from_array.
    1264         #
    1265         # Alternatively, some C code would be handy
    1266         #
    1267         self._set_vertex_values(vertex_list, A)
    1268            
    1269 
    1270     def _set_vertex_values(self, vertex_list, A):
    1271         """Go through list of unique vertices
    1272         This is the common case e.g. when values
    1273         are obtained from a pts file through fitting
    1274         """
    1275        
     1236        # Go through list of unique vertices
    12761237        for i_index, unique_vert_id in enumerate(vertex_list):
    1277 
    12781238            triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
    1279                    
     1239
    12801240            # In case there are unused points
    1281             if len(triangles) == 0: continue
     1241            if len(triangles) == 0:
     1242                continue
    12821243
    12831244            # Go through all triangle, vertex pairs
     
    12891250        self.interpolate()
    12901251
    1291 
    1292     def smooth_vertex_values(self,
    1293                              use_cache=False,
    1294                              verbose=False):
    1295         """ Smooths vertex values.
    1296         """
    1297 
    1298         A,V = self.get_vertex_values(xy=False, smooth=True)
    1299         self.set_vertex_values(A,
    1300                                use_cache=use_cache,
    1301                                verbose=verbose)                               
    1302                                
    1303 
    1304 
     1252    ##
     1253    # @brief Smooth vertex values.
     1254    def smooth_vertex_values(self):
     1255        """Smooths vertex values."""
     1256
     1257        A, V = self.get_vertex_values(xy=False, smooth=True)
     1258        self.set_vertex_values(A)
     1259
     1260    ############################################################################
    13051261    # Methods for outputting model results
    1306     def get_vertex_values(self,
    1307                           xy=True,
    1308                           smooth=None,
    1309                           precision=None):
     1262    ############################################################################
     1263
     1264    ##
     1265    # @brief Get vertex values like an OBJ format i.e. one value per node.
     1266    # @param xy True if we return X and Y as well as A and V.
     1267    # @param smooth True if vertex values are to be smoothed.
     1268    # @param precision The type of the result values (default float).
     1269    def get_vertex_values(self, xy=True,
     1270                                smooth=None,
     1271                                precision=None):
    13101272        """Return vertex values like an OBJ format i.e. one value per node.
    13111273
     
    13161278        Mx3, where M is the number of triangles. Each row has three indices
    13171279        defining the triangle and they correspond to elements in the arrays
    1318         X, Y and A.
    1319 
    1320         if smooth is True, vertex values corresponding to one common
    1321         coordinate set will be smoothed by taking the average of vertex values for each node.
    1322         In this case vertex coordinates will be
    1323         de-duplicated corresponding to the original nodes as obtained from
    1324         the method general_mesh.get_nodes()
    1325 
    1326         If no smoothings is required, vertex coordinates and values will
    1327         be aggregated as a concatenation of values at
    1328         vertices 0, vertices 1 and vertices 2. This corresponds to
    1329         the node coordinates obtained from the method
    1330         general_mesh.get_vertex_coordinates()
    1331 
     1280        X, Y and A.
     1281
     1282        If smooth is True, vertex values corresponding to one common coordinate
     1283        set will be smoothed by taking the average of vertex values for each
     1284        node.  In this case vertex coordinates will be de-duplicated
     1285        corresponding to the original nodes as obtained from the method
     1286        general_mesh.get_nodes()
     1287
     1288        If no smoothings is required, vertex coordinates and values will be
     1289        aggregated as a concatenation of values at vertices 0, vertices 1 and
     1290        vertices 2.  This corresponds to the node coordinates obtained from the
     1291        method general_mesh.get_vertex_coordinates()
    13321292
    13331293        Calling convention
    13341294        if xy is True:
    1335            X,Y,A,V = get_vertex_values
     1295           X, Y, A, V = get_vertex_values
    13361296        else:
    1337            A,V = get_vertex_values
    1338 
    1339         """
    1340 
     1297           A, V = get_vertex_values
     1298        """
    13411299
    13421300        if smooth is None:
     
    13491307        if precision is None:
    13501308            precision = num.Float
    1351            
    13521309
    13531310        if smooth is True:
    1354             # Ensure continuous vertex values by averaging
    1355             # values at each node
    1356            
     1311            # Ensure continuous vertex values by averaging values at each node
    13571312            V = self.domain.get_triangles()
    13581313            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    13591314            A = num.zeros(N, num.Float)
    1360             points = self.domain.get_nodes()           
    1361            
     1315            points = self.domain.get_nodes()
     1316
    13621317            if 1:
    13631318                # Fast C version
     
    13671322                                      A)
    13681323                A = A.astype(precision)
    1369             else:   
    1370 
     1324            else:
    13711325                # Slow Python version
    1372                
    13731326                current_node = 0
    13741327                k = 0 # Track triangles touching on node
     
    13761329                for index in self.domain.vertex_value_indices:
    13771330                    if current_node == N:
    1378                         msg = 'Current node exceeding number of nodes (%d) ' %(N)
     1331                        msg = 'Current node exceeding number of nodes (%d) ' % N
    13791332                        raise msg
    1380                    
    1381 
    1382                    
     1333
    13831334                    k += 1
    1384                    
     1335
    13851336                    volume_id = index / 3
    13861337                    vertex_id = index % 3
    1387                  
    1388                     #assert V[volume_id, vertex_id] == current_node
    1389                
     1338
    13901339                    v = self.vertex_values[volume_id, vertex_id]
    13911340                    total += v
    13921341
    1393                     #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
    13941342                    if self.domain.number_of_triangles_per_node[current_node] == k:
    13951343                        A[current_node] = total/k
    1396                
    1397                    
     1344
    13981345                        # Move on to next node
    13991346                        total = 0.0
    14001347                        k = 0
    14011348                        current_node += 1
    1402 
    1403 
    1404 
    14051349        else:
    1406             # Return disconnected internal vertex values 
     1350            # Return disconnected internal vertex values
    14071351            V = self.domain.get_disconnected_triangles()
    14081352            points = self.domain.get_vertex_coordinates()
    14091353            A = self.vertex_values.flat.astype(precision)
    14101354
    1411 
    1412         # Return   
     1355        # Return
    14131356        if xy is True:
    14141357            X = points[:,0].astype(precision)
    14151358            Y = points[:,1].astype(precision)
    1416            
     1359
    14171360            return X, Y, A, V
    14181361        else:
    1419             return A, V           
    1420 
    1421 
    1422 
     1362            return A, V
     1363
     1364    ##
     1365    # @brief Extrapolate conserved quantities from centroid.
    14231366    def extrapolate_first_order(self):
    1424         """Extrapolate conserved quantities from centroid to
    1425         vertices and edges for each volume using
    1426         first order scheme.
     1367        """Extrapolate conserved quantities from centroid to vertices and edges
     1368        for each volume using first order scheme.
    14271369        """
    14281370
     
    14381380        self.y_gradient *= 0.0
    14391381
    1440 
     1382    ##
     1383    # @brief Compute the integral of quantity across entire domain.
     1384    # @return The integral.
    14411385    def get_integral(self):
    1442         """Compute the integral of quantity across entire domain
    1443         """
     1386        """Compute the integral of quantity across entire domain."""
     1387
    14441388        areas = self.domain.get_areas()
    14451389        integral = 0
     
    14511395        return integral
    14521396
     1397    ##
     1398    # @brief get the gradients.
    14531399    def get_gradients(self):
    1454         """Provide gradients. Use compute_gradients first
    1455         """
     1400        """Provide gradients. Use compute_gradients first."""
    14561401
    14571402        return self.x_gradient, self.y_gradient
    14581403
    1459 
     1404    ##
     1405    # @brief ??
     1406    # @param timestep ??
    14601407    def update(self, timestep):
    14611408        # Call correct module function
     
    14631410        return update(self, timestep)
    14641411
     1412    ##
     1413    # @brief ??
    14651414    def compute_gradients(self):
    14661415        # Call correct module function
     
    14681417        return compute_gradients(self)
    14691418
     1419    ##
     1420    # @brief ??
    14701421    def limit(self):
    14711422        # Call correct module depending on whether
     
    14731424        limit_old(self)
    14741425
     1426    ##
     1427    # @brief ??
    14751428    def limit_vertices_by_all_neighbours(self):
    14761429        # Call correct module function
     
    14781431        limit_vertices_by_all_neighbours(self)
    14791432
     1433    ##
     1434    # @brief ??
    14801435    def limit_edges_by_all_neighbours(self):
    14811436        # Call correct module function
     
    14831438        limit_edges_by_all_neighbours(self)
    14841439
     1440    ##
     1441    # @brief ??
    14851442    def limit_edges_by_neighbour(self):
    14861443        # Call correct module function
    14871444        # (either from this module or C-extension)
    1488         limit_edges_by_neighbour(self)               
    1489 
     1445        limit_edges_by_neighbour(self)
     1446
     1447    ##
     1448    # @brief ??
    14901449    def extrapolate_second_order(self):
    14911450        # Call correct module function
     
    14931452        compute_gradients(self)
    14941453        extrapolate_from_gradient(self)
    1495        
     1454
     1455    ##
     1456    # @brief ??
    14961457    def extrapolate_second_order_and_limit_by_edge(self):
    14971458        # Call correct module function
     
    14991460        extrapolate_second_order_and_limit_by_edge(self)
    15001461
     1462    ##
     1463    # @brief ??
    15011464    def extrapolate_second_order_and_limit_by_vertex(self):
    15021465        # Call correct module function
     
    15041467        extrapolate_second_order_and_limit_by_vertex(self)
    15051468
     1469    ##
     1470    # @brief ??
     1471    # @param bound ??
    15061472    def bound_vertices_below_by_constant(self, bound):
    15071473        # Call correct module function
     
    15091475        bound_vertices_below_by_constant(self, bound)
    15101476
     1477    ##
     1478    # @brief ??
     1479    # @param quantity ??
    15111480    def bound_vertices_below_by_quantity(self, quantity):
    15121481        # Call correct module function
     
    15151484        # check consistency
    15161485        assert self.domain == quantity.domain
    1517         bound_vertices_below_by_quantity(self, quantity)                       
    1518 
     1486        bound_vertices_below_by_quantity(self, quantity)
     1487
     1488    ##
     1489    # @brief ??
    15191490    def backup_centroid_values(self):
    15201491        # Call correct module function
     
    15221493        backup_centroid_values(self)
    15231494
    1524     def saxpy_centroid_values(self,a,b):
     1495    ##
     1496    # @brief ??
     1497    # @param a ??
     1498    # @param b ??
     1499    def saxpy_centroid_values(self, a, b):
    15251500        # Call correct module function
    15261501        # (either from this module or C-extension)
    1527         saxpy_centroid_values(self,a,b)
    1528    
    1529 #Conserved_quantity = Quantity
    1530 
     1502        saxpy_centroid_values(self, a, b)
     1503
     1504
     1505##
     1506# @brief OBSOLETE!
    15311507class Conserved_quantity(Quantity):
    1532     """Class conserved quantity being removed, use Quantity
    1533 
    1534     """
     1508    """Class conserved quantity being removed, use Quantity."""
    15351509
    15361510    def __init__(self, domain, vertex_values=None):
    1537         #Quantity.__init__(self, domain, vertex_values)
    1538 
    15391511        msg = 'ERROR: Use Quantity instead of Conserved_quantity'
    1540 
    15411512        raise Exception, msg
    15421513
    15431514
     1515######
     1516# Prepare the C extensions.
     1517######
    15441518
    15451519from anuga.utilities import compile
    1546 if compile.can_use_C_extension('quantity_ext.c'):   
    1547     # Underlying C implementations can be accessed
     1520
     1521if compile.can_use_C_extension('quantity_ext.c'):
     1522    # Underlying C implementations can be accessed
    15481523
    15491524    from quantity_ext import \
     
    15641539         interpolate_from_vertices_to_edges,\
    15651540         interpolate_from_edges_to_vertices,\
    1566          update   
     1541         update
    15671542else:
    1568     msg = 'C implementations could not be accessed by %s.\n ' %__file__
     1543    msg = 'C implementations could not be accessed by %s.\n ' % __file__
    15691544    msg += 'Make sure compile_all.py has been run as described in '
    15701545    msg += 'the ANUGA installation guide.'
    15711546    raise Exception, msg
    1572 
    1573 
    1574 
Note: See TracChangeset for help on using the changeset viewer.