source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py @ 7317

Last change on this file since 7317 was 7317, checked in by rwilson, 15 years ago

Replaced 'print' statements with log.critical() calls.

File size: 58.6 KB
RevLine 
[5897]1"""Class Quantity - Implements values at each triangular element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 3 array of values at each vertex for each element.
10                  Default None
11
12   If vertex_values are None Create array of zeros compatible with domain.
13   Otherwise check that it is compatible with dimenions of domain.
14   Otherwise raise an exception
15"""
16
[7276]17from types import FloatType, IntType, LongType, NoneType
18
[5897]19from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
20from anuga.utilities.polygon import inside_polygon
21from anuga.geospatial_data.geospatial_data import Geospatial_data
22from anuga.fit_interpolate.fit import fit_to_mesh
23from anuga.config import points_file_block_line_size as default_block_line_size
24from anuga.config import epsilon
[6228]25from anuga.caching import cache
[7317]26import anuga.utilities.log as log
[5897]27
[7276]28import anuga.utilities.numerical_tools as aunt
[6145]29
[7276]30import numpy as num
[6145]31
[7276]32
[6226]33##
34# @brief Implement values at each triangular element.
[5897]35class Quantity:
36
[6226]37    ##
38    # @brief Construct values art each triangular element.
39    # @param domain ??
40    # @param vertex_values ??
[6228]41    def __init__(self, domain, vertex_values=None):
[6191]42        from anuga.abstract_2d_finite_volumes.domain import Domain
[6226]43
[6228]44        msg = ('First argument in Quantity.__init__() must be of class Domain '
45               '(or a subclass thereof). I got %s.' % str(domain.__class__))
[6195]46        assert isinstance(domain, Domain), msg
[5897]47
48        if vertex_values is None:
[6226]49            N = len(domain)             # number_of_elements
[7276]50            self.vertex_values = num.zeros((N, 3), num.float)
[5897]51        else:
[7276]52            self.vertex_values = num.array(vertex_values, num.float)
[5897]53
54            N, V = self.vertex_values.shape
[6226]55            assert V == 3, 'Three vertex values per element must be specified'
[5897]56
[6226]57            msg = 'Number of vertex values (%d) must be consistent with' % N
58            msg += 'number of elements in specified domain (%d).' % len(domain)
[5897]59            assert N == len(domain), msg
60
61        self.domain = domain
62
63        # Allocate space for other quantities
[7276]64        self.centroid_values = num.zeros(N, num.float)
65        self.edge_values = num.zeros((N, 3), num.float)
[5897]66
67        # Allocate space for Gradient
[7276]68        self.x_gradient = num.zeros(N, num.float)
69        self.y_gradient = num.zeros(N, num.float)
[5897]70
71        # Allocate space for Limiter Phi
[7276]72        self.phi = num.zeros(N, num.float)
[5897]73
74        # Intialise centroid and edge_values
75        self.interpolate()
76
77        # Allocate space for boundary values
78        L = len(domain.boundary)
[7276]79        self.boundary_values = num.zeros(L, num.float)
[5897]80
81        # Allocate space for updates of conserved quantities by
82        # flux calculations and forcing functions
83
84        # Allocate space for update fields
[7276]85        self.explicit_update = num.zeros(N, num.float )
86        self.semi_implicit_update = num.zeros(N, num.float )
87        self.centroid_backup_values = num.zeros(N, num.float)
[5897]88
89        self.set_beta(1.0)
90
[6228]91    ############################################################################
[6226]92    # Methods for operator overloading
[6228]93    ############################################################################
[5897]94
95    def __len__(self):
96        return self.centroid_values.shape[0]
97
98    def __neg__(self):
99        """Negate all values in this quantity giving meaning to the
100        expression -Q where Q is an instance of class Quantity
101        """
102
103        Q = Quantity(self.domain)
104        Q.set_values(-self.vertex_values)
105        return Q
106
107    def __add__(self, other):
108        """Add to self anything that could populate a quantity
109
110        E.g other can be a constant, an array, a function, another quantity
111        (except for a filename or points, attributes (for now))
112        - see set_values for details
113        """
114
115        Q = Quantity(self.domain)
116        Q.set_values(other)
117
118        result = Quantity(self.domain)
119        result.set_values(self.vertex_values + Q.vertex_values)
120        return result
121
122    def __radd__(self, other):
123        """Handle cases like 7+Q, where Q is an instance of class Quantity
124        """
[6226]125
[5897]126        return self + other
127
128    def __sub__(self, other):
[6226]129        return self + -other            # Invoke self.__neg__()
[5897]130
131    def __mul__(self, other):
132        """Multiply self with anything that could populate a quantity
133
134        E.g other can be a constant, an array, a function, another quantity
135        (except for a filename or points, attributes (for now))
136        - see set_values for details
137        """
138
139        if isinstance(other, Quantity):
140            Q = other
[6226]141        else:
[5897]142            Q = Quantity(self.domain)
143            Q.set_values(other)
144
145        result = Quantity(self.domain)
146
147        # The product of vertex_values, edge_values and centroid_values
148        # are calculated and assigned directly without using
149        # set_values (which calls interpolate). Otherwise
150        # edge and centroid values wouldn't be products from q1 and q2
151        result.vertex_values = self.vertex_values * Q.vertex_values
152        result.edge_values = self.edge_values * Q.edge_values
153        result.centroid_values = self.centroid_values * Q.centroid_values
[6226]154
[5897]155        return result
156
157    def __rmul__(self, other):
158        """Handle cases like 3*Q, where Q is an instance of class Quantity
159        """
[6226]160
[5897]161        return self * other
162
163    def __div__(self, other):
164        """Divide self with anything that could populate a quantity
165
166        E.g other can be a constant, an array, a function, another quantity
167        (except for a filename or points, attributes (for now))
168        - see set_values for details
169
170        Zero division is dealt with by adding an epsilon to the divisore
171        FIXME (Ole): Replace this with native INF once we migrate to NumPy
172        """
173
174        if isinstance(other, Quantity):
175            Q = other
[6226]176        else:
[5897]177            Q = Quantity(self.domain)
178            Q.set_values(other)
179
180        result = Quantity(self.domain)
181
182        # The quotient of vertex_values, edge_values and centroid_values
183        # are calculated and assigned directly without using
184        # set_values (which calls interpolate). Otherwise
185        # edge and centroid values wouldn't be quotient of q1 and q2
186        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
187        result.edge_values = self.edge_values/(Q.edge_values + epsilon)
188        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
189
190        return result
191
192    def __rdiv__(self, other):
193        """Handle cases like 3/Q, where Q is an instance of class Quantity
194        """
[6226]195
[5897]196        return self / other
197
198    def __pow__(self, other):
199        """Raise quantity to (numerical) power
200
201        As with __mul__ vertex values are processed entry by entry
202        while centroid and edge values are re-interpolated.
203
204        Example using __pow__:
205          Q = (Q1**2 + Q2**2)**0.5
206        """
207
208        if isinstance(other, Quantity):
209            Q = other
[6226]210        else:
[5897]211            Q = Quantity(self.domain)
212            Q.set_values(other)
213
214        result = Quantity(self.domain)
215
216        # The power of vertex_values, edge_values and centroid_values
217        # are calculated and assigned directly without using
218        # set_values (which calls interpolate). Otherwise
219        # edge and centroid values wouldn't be correct
220        result.vertex_values = self.vertex_values ** other
221        result.edge_values = self.edge_values ** other
222        result.centroid_values = self.centroid_values ** other
223
224        return result
225
[6228]226    ############################################################################
227    # Setters/Getters
228    ############################################################################
229
[6226]230    ##
231    # @brief Set default beta value for limiting.
232    # @param beta ??
233    def set_beta(self, beta):
234        """Set default beta value for limiting """
[5897]235
236        if beta < 0.0:
[7317]237            log.critical('WARNING: setting beta < 0.0')
[5897]238        if beta > 2.0:
[7317]239            log.critical('WARNING: setting beta > 2.0')
[6226]240
[5897]241        self.beta = beta
242
[6226]243    ##
244    # @brief Get the current beta value.
245    # @return The current beta value.
[5897]246    def get_beta(self):
[6226]247        """Get default beta value for limiting"""
[5897]248
249        return self.beta
250
[6226]251    ##
252    # @brief Compute interpolated values at edges and centroid.
253    # @note vertex_values must be set before calling this.
[5897]254    def interpolate(self):
255        """Compute interpolated values at edges and centroid
256        Pre-condition: vertex_values have been set
257        """
[6226]258
[5897]259        # FIXME (Ole): Maybe this function
260        # should move to the C-interface?
261        # However, it isn't called by validate_all.py, so it
262        # may not be that important to optimise it?
[6226]263
[5897]264        N = self.vertex_values.shape[0]
265        for i in range(N):
266            v0 = self.vertex_values[i, 0]
267            v1 = self.vertex_values[i, 1]
268            v2 = self.vertex_values[i, 2]
269
270            self.centroid_values[i] = (v0 + v1 + v2)/3
271
272        self.interpolate_from_vertices_to_edges()
273
[6226]274    ##
275    # @brief ??
[5897]276    def interpolate_from_vertices_to_edges(self):
[6226]277        # Call correct module function (either from this module or C-extension)
[5897]278        interpolate_from_vertices_to_edges(self)
279
[6226]280    ##
281    # @brief ??
[5897]282    def interpolate_from_edges_to_vertices(self):
[6226]283        # Call correct module function (either from this module or C-extension)
[5897]284        interpolate_from_edges_to_vertices(self)
285
286    #---------------------------------------------
287    # Public interface for setting quantity values
288    #---------------------------------------------
289
[6226]290    ##
291    # @brief Set values for quantity based on different sources.
292    # @param numeric A num array, list or constant value.
293    # @param quantity Another Quantity.
294    # @param function Any callable object that takes two 1d arrays.
295    # @param geospatial_data Arbitrary instance of class Geospatial_data
296    # @param filename Path to a points file.
297    # @param attribute_name If specified any array using that name will be used.
298    # @param alpha Smoothing parameter to be used with fit_interpolate.fit.
299    # @param location Where to store values (vertices, edges, centroids).
300    # @param polygon Restrict update to locations that fall inside polygon.
301    # @param indices Restrict update to locations specified by this.
302    # @param smooth If True, smooth vertex values.
303    # @param verbose True if this method is to be verbose.
304    # @param use_cache If True cache results for fit_interpolate.fit.
305    # @note Exactly one of 'numeric', 'quantity', 'function', 'filename'
306    #       must be present.
307    def set_values(self, numeric=None,         # List, numeric array or constant
308                         quantity=None,        # Another quantity
309                         function=None,        # Callable object: f(x,y)
310                         geospatial_data=None, # Arbitrary dataset
311                         filename=None,
312                         attribute_name=None,  # Input from file
313                         alpha=None,
314                         location='vertices',
315                         polygon=None,
316                         indices=None,
317                         smooth=False,
318                         verbose=False,
319                         use_cache=False):
[5897]320        """Set values for quantity based on different sources.
321
322        numeric:
[7276]323          Compatible list, numeric array (see below) or constant.
[5897]324          If callable it will treated as a function (see below)
325          If instance of another Quantity it will be treated as such.
326          If geo_spatial object it will be treated as such
327
328        quantity:
329          Another quantity (compatible quantity, e.g. obtained as a
330          linear combination of quantities)
331
332        function:
333          Any callable object that takes two 1d arrays x and y
334          each of length N and returns an array also of length N.
335          The function will be evaluated at points determined by
336          location and indices in the underlying mesh.
337
338        geospatial_data:
339          Arbitrary geo spatial dataset in the form of the class
340          Geospatial_data. Mesh points are populated using
341          fit_interpolate.fit fitting
342
343        filename:
344          Name of a points file containing data points and attributes for
345          use with fit_interpolate.fit.
346
347        attribute_name:
348          If specified, any array matching that name
349          will be used. from file or geospatial_data.
350          Otherwise a default will be used.
351
352        alpha:
353          Smoothing parameter to be used with fit_interpolate.fit.
354          See module fit_interpolate.fit for further details about alpha.
355          Alpha will only be used with points, values or filename.
356          Otherwise it will be ignored.
357
358
359        location: Where values are to be stored.
360                  Permissible options are: vertices, edges, centroids
361                  Default is 'vertices'
362
363                  In case of location == 'centroids' the dimension values must
[7276]364                  be a list of a numerical array of length N,
[5897]365                  N being the number of elements.
366                  Otherwise it must be of dimension Nx3
367
368
369                  The values will be stored in elements following their
370                  internal ordering.
371
372                  If location is 'unique vertices' indices refers the set
373                  of node ids that the operation applies to.
374                  If location is not 'unique vertices' indices refers the
375                  set of triangle ids that the operation applies to.
376
377
378                  If selected location is vertices, values for
379                  centroid and edges will be assigned interpolated
380                  values.  In any other case, only values for the
381                  specified locations will be assigned and the others
382                  will be left undefined.
383
384
385        polygon: Restrict update of quantity to locations that fall
386                 inside polygon. Polygon works by selecting indices
387                 and calling set_values recursively.
388                 Polygon mode has only been implemented for
389                 constant values so far.
390
[6226]391        indices: Restrict update of quantity to locations that are
[5897]392                 identified by indices (e.g. node ids if location
393                 is 'unique vertices' or triangle ids otherwise).
[6226]394
[5897]395        verbose: True means that output to stdout is generated
396
397        use_cache: True means that caching of intermediate results is
398                   attempted for fit_interpolate.fit.
399
400
401
402
403        Exactly one of the arguments
404          numeric, quantity, function, filename
405        must be present.
406        """
407
408        from anuga.geospatial_data.geospatial_data import Geospatial_data
409
410        # Treat special case: Polygon situation
411        # Location will be ignored and set to 'centroids'
412        # FIXME (Ole): This needs to be generalised and
413        # perhaps the notion of location and indices simplified
414
[6226]415        # FIXME (Ole): Need to compute indices based on polygon
[5990]416        # (and location) and use existing code after that.
[6226]417
[5990]418        # See ticket:275, ticket:250, ticeket:254 for refactoring plan
[6226]419
[5897]420        if polygon is not None:
421            if indices is not None:
422                msg = 'Only one of polygon and indices can be specified'
423                raise Exception, msg
424
425            msg = 'With polygon selected, set_quantity must provide '
426            msg += 'the keyword numeric and it must (currently) be '
427            msg += 'a constant.'
428            if numeric is None:
[6226]429                raise Exception, msg
[5897]430            else:
431                # Check that numeric is as constant
432                assert type(numeric) in [FloatType, IntType, LongType], msg
433
434            location = 'centroids'
435
436            points = self.domain.get_centroid_coordinates(absolute=True)
437            indices = inside_polygon(points, polygon)
438
[6226]439            self.set_values_from_constant(numeric, location, indices, verbose)
[5897]440
441            self.extrapolate_first_order()
442
443            if smooth:
[6228]444                self.smooth_vertex_values(use_cache=use_cache,
445                                          verbose=verbose)
[5897]446
447            return
448
449        # General input checks
450        L = [numeric, quantity, function, geospatial_data, filename]
[6226]451        msg = ('Exactly one of the arguments numeric, quantity, function, '
452               'geospatial_data, or filename must be present.')
[5897]453        assert L.count(None) == len(L)-1, msg
454
455        if location == 'edges':
456            msg = 'edges has been deprecated as valid location'
457            raise Exception, msg
[6226]458
[5897]459        if location not in ['vertices', 'centroids', 'unique vertices']:
[6226]460            msg = 'Invalid location: %s' % location
[5897]461            raise Exception, msg
462
[7276]463        msg = 'Indices must be a list, array or None'
464        assert isinstance(indices, (NoneType, list, num.ndarray)), msg
[5897]465
466        # Determine which 'set_values_from_...' to use
467        if numeric is not None:
[7276]468            if isinstance(numeric, (list, num.ndarray)):
[6228]469                self.set_values_from_array(numeric, location, indices,
470                                           use_cache=use_cache, verbose=verbose)
[5897]471            elif callable(numeric):
[6228]472                self.set_values_from_function(numeric, location, indices,
473                                              use_cache=use_cache,
474                                              verbose=verbose)
[5897]475            elif isinstance(numeric, Quantity):
[6228]476                self.set_values_from_quantity(numeric, location, indices,
477                                              verbose=verbose)
[5897]478            elif isinstance(numeric, Geospatial_data):
[6226]479                self.set_values_from_geospatial_data(numeric, alpha, location,
480                                                     indices, verbose=verbose,
[5897]481                                                     use_cache=use_cache)
[7276]482            else:   # see if it's coercible to a float (float, int or long, etc)
483                try:
484                    numeric = float(numeric)
485                except ValueError:
486                    msg = ("Illegal type for variable 'numeric': %s"
487                           % type(numeric))
488                    raise Exception, msg
489                self.set_values_from_constant(numeric, location,
490                                              indices, verbose)
[5897]491        elif quantity is not None:
[6226]492            self.set_values_from_quantity(quantity, location, indices, verbose)
[5897]493        elif function is not None:
494            msg = 'Argument function must be callable'
495            assert callable(function), msg
[6228]496            self.set_values_from_function(function, location, indices,
497                                          use_cache=use_cache, verbose=verbose)
[5897]498        elif geospatial_data is not None:
[6226]499                self.set_values_from_geospatial_data(geospatial_data, alpha,
[5897]500                                                     location, indices,
501                                                     verbose=verbose,
502                                                     use_cache=use_cache)
503        elif filename is not None:
504            if hasattr(self.domain, 'points_file_block_line_size'):
505                max_read_lines = self.domain.points_file_block_line_size
506            else:
507                max_read_lines = default_block_line_size
[6226]508            self.set_values_from_file(filename, attribute_name, alpha, location,
509                                      indices, verbose=verbose,
[5897]510                                      max_read_lines=max_read_lines,
511                                      use_cache=use_cache)
512        else:
[6226]513            raise Exception, "This can't happen :-)"
[5897]514
515        # Update all locations in triangles
516        if location == 'vertices' or location == 'unique vertices':
517            # Intialise centroid and edge_values
518            self.interpolate()
519
520        if location == 'centroids':
521            # Extrapolate 1st order - to capture notion of area being specified
522            self.extrapolate_first_order()
523
[6228]524    ############################################################################
[5897]525    # Specific internal functions for setting values based on type
[6228]526    ############################################################################
[5897]527
[6226]528    ##
529    # @brief Set quantity values from specified constant.
530    # @param X The constant to set quantity values to.
531    # @param location
532    # @param indices
533    # @param verbose
534    def set_values_from_constant(self, X, location, indices, verbose):
535        """Set quantity values from specified constant X"""
536
[5897]537        # FIXME (Ole): Somehow indices refer to centroids
538        # rather than vertices as default. See unit test
539        # test_set_vertex_values_using_general_interface_with_subset(self):
540
541        if location == 'centroids':
542            if indices is None:
543                self.centroid_values[:] = X
544            else:
545                # Brute force
546                for i in indices:
547                    self.centroid_values[i] = X
548        elif location == 'unique vertices':
549            if indices is None:
550                self.edge_values[:] = X  #FIXME (Ole): Shouldn't this be vertex_values?
551            else:
552                # Go through list of unique vertices
553                for unique_vert_id in indices:
[6226]554                    triangles = \
555                        self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
[5897]556
557                    # In case there are unused points
558                    if len(triangles) == 0:
559                        continue
[6226]560
[5897]561                    # Go through all triangle, vertex pairs
562                    # and set corresponding vertex value
563                    for triangle_id, vertex_id in triangles:
564                        self.vertex_values[triangle_id, vertex_id] = X
565
566                    # Intialise centroid and edge_values
567                    self.interpolate()
568        else:
569            if indices is None:
570                self.vertex_values[:] = X
571            else:
572                # Brute force
573                for i_vertex in indices:
574                    self.vertex_values[i_vertex] = X
575
[6226]576    ##
577    # @brief Set values for a quantity.
578    # @param values Array of values.
579    # @param location Where values are to be stored.
580    # @param indices Limit update to these indices.
[6228]581    # @param use_cache ??
[6226]582    # @param verbose True if this method is to be verbose.
[5897]583    def set_values_from_array(self, values,
[6226]584                                    location='vertices',
585                                    indices=None,
[6228]586                                    use_cache=False,
[6226]587                                    verbose=False):
[5897]588        """Set values for quantity
589
[7276]590        values: numeric array
[5897]591        location: Where values are to be stored.
592        Permissible options are: vertices, centroid, unique vertices
593        Default is 'vertices'
594
595        indices - if this action is carried out on a subset of
596        elements or unique vertices
597        The element/unique vertex indices are specified here.
598
599        In case of location == 'centroid' the dimension values must
[7276]600        be a list of a numerical array of length N, N being the number
[5897]601        of elements.
602
603        Otherwise it must be of dimension Nx3
604
605        The values will be stored in elements following their
606        internal ordering.
607
608        If selected location is vertices, values for centroid and edges
609        will be assigned interpolated values.
610        In any other case, only values for the specified locations
611        will be assigned and the others will be left undefined.
612        """
613
[7276]614        values = num.array(values, num.float)
[5897]615
616        if indices is not None:
[7276]617            indices = num.array(indices, num.int)
[6226]618            msg = ('Number of values must match number of indices: You '
619                   'specified %d values and %d indices'
620                   % (values.shape[0], indices.shape[0]))
[5897]621            assert values.shape[0] == indices.shape[0], msg
622
623        N = self.centroid_values.shape[0]
624
625        if location == 'centroids':
626            assert len(values.shape) == 1, 'Values array must be 1d'
627
628            if indices is None:
629                msg = 'Number of values must match number of elements'
630                assert values.shape[0] == N, msg
631
632                self.centroid_values = values
633            else:
634                msg = 'Number of values must match number of indices'
635                assert values.shape[0] == indices.shape[0], msg
636
637                # Brute force
638                for i in range(len(indices)):
639                    self.centroid_values[indices[i]] = values[i]
640        elif location == 'unique vertices':
[6226]641            assert (len(values.shape) == 1 or num.allclose(values.shape[1:], 1),
642                    'Values array must be 1d')
[5897]643
[7276]644            self.set_vertex_values(values.flatten(), indices=indices,
[6228]645                                   use_cache=use_cache, verbose=verbose)
[5897]646        else:
647            # Location vertices
648            if len(values.shape) == 1:
[6228]649                # This is the common case arising from fitted
650                # values (e.g. from pts file).
651                self.set_vertex_values(values, indices=indices,
652                                       use_cache=use_cache, verbose=verbose)
[5897]653            elif len(values.shape) == 2:
654                # Vertex values are given as a triplet for each triangle
655                msg = 'Array must be N x 3'
656                assert values.shape[1] == 3, msg
657
658                if indices is None:
659                    self.vertex_values = values
660                else:
661                    for element_index, value in map(None, indices, values):
662                        self.vertex_values[element_index] = value
663            else:
664                msg = 'Values array must be 1d or 2d'
[6226]665                raise Exception, msg
[5897]666
[6226]667    ##
668    # @brief Set quantity values from a specified quantity instance.
669    # @param q The quantity instance to take values from.
670    # @param location IGNORED, 'vertices' ALWAYS USED!
671    # @param indices ??
672    # @param verbose True if this method is to be verbose.
673    def set_values_from_quantity(self, q, location, indices, verbose):
[5897]674        """Set quantity values from specified quantity instance q
675
676        Location is ignored - vertices will always be used here.
677        """
678
679
680        A = q.vertex_values
681
682        msg = 'Quantities are defined on different meshes. '+\
683              'This might be a case for implementing interpolation '+\
684              'between different meshes.'
[6145]685        assert num.allclose(A.shape, self.vertex_values.shape), msg
[5897]686
687        self.set_values(A, location='vertices',
[6226]688                        indices=indices, verbose=verbose)
[5897]689
[6226]690    ##
691    # @brief Set quantity values from a specified quantity instance.
692    # @param f Callable that takes two 1d array -> 1d array.
693    # @param location Where values are to be stored.
694    # @param indices ??
[6228]695    # @param use_cache ??
[6226]696    # @param verbose True if this method is to be verbose.
[7276]697    def set_values_from_function(self,
698                                 f,
699                                 location='vertices',
700                                 indices=None,
701                                 use_cache=False,
702                                 verbose=False):
[5897]703        """Set values for quantity using specified function
704
705        Input
706        f: x, y -> z Function where x, y and z are arrays
707        location: Where values are to be stored.
708                  Permissible options are: vertices, centroid,
709                  unique vertices
710                  Default is "vertices"
[6226]711        indices:
[5897]712        """
713
714        # FIXME: Should check that function returns something sensible and
[6226]715        # raise a meaningful exception if it returns None for example
[5897]716
717        # FIXME: Should supply absolute coordinates
718
719        # Compute the function values and call set_values again
720        if location == 'centroids':
721            if indices is None:
722                indices = range(len(self))
[6226]723
[7276]724            V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0)
[6228]725            x = V[:,0]; y = V[:,1]
726            if use_cache is True:
727                res = cache(f, (x, y), verbose=verbose)
728            else:
729                res = f(x, y)
730
731            self.set_values(res, location=location, indices=indices)
[5897]732        elif location == 'vertices':
[6228]733            # This is the default branch taken by set_quantity
[5897]734            M = self.domain.number_of_triangles
735            V = self.domain.get_vertex_coordinates()
736
[6226]737            x = V[:,0];
[6228]738            y = V[:,1]
739            if use_cache is True:
740                values = cache(f, (x, y), verbose=verbose)               
741            else:
742                if verbose is True:
[7317]743                    log.critical('Evaluating function in set_values')
[6228]744                values = f(x, y)
[5897]745
746            # FIXME (Ole): This code should replace all the
747            # rest of this function and it would work, except
748            # one unit test in test_region fails.
749            # If that could be resolved this one will be
750            # more robust and simple.
751
752            # This should be removed
753            if is_scalar(values):
754                # Function returned a constant value
[6226]755                self.set_values_from_constant(values, location,
756                                              indices, verbose)
[5897]757                return
758
[6226]759            # This should be removed
[5897]760            if indices is None:
761                for j in range(3):
[6226]762                    self.vertex_values[:, j] = values[j::3]
763            else:
[5897]764                # Brute force
765                for i in indices:
766                    for j in range(3):
[6226]767                        self.vertex_values[i, j] = values[3*i + j]
[5897]768        else:
[6226]769            raise Exception, 'Not implemented: %s' % location
[5897]770
[6226]771    ##
772    # @brief Set values based on geo referenced geospatial data object.
773    # @param geospatial_data ??
774    # @param alpha ??
775    # @param location ??
776    # @param indices ??
777    # @param verbose ??
778    # @param use_cache ??
[7276]779    def set_values_from_geospatial_data(self,
780                                        geospatial_data,
781                                        alpha,
782                                        location,
783                                        indices,
784                                        verbose=False,
785                                        use_cache=False):
[6226]786        """Set values based on geo referenced geospatial data object."""
[5897]787
[6226]788        from anuga.coordinate_transforms.geo_reference import Geo_reference
[5897]789
790        points = geospatial_data.get_data_points(absolute=False)
791        values = geospatial_data.get_attributes()
792        data_georef = geospatial_data.get_geo_reference()
793
[6228]794        from anuga.coordinate_transforms.geo_reference import Geo_reference
795
[7276]796        points = ensure_numeric(points, num.float)
797        values = ensure_numeric(values, num.float)
[5897]798
799        if location != 'vertices':
[6226]800            msg = ("set_values_from_points is only defined for "
801                   "location='vertices'")
802            raise Exception, msg
[5897]803
804        # Take care of georeferencing
805        if data_georef is None:
806            data_georef = Geo_reference()
807
808        mesh_georef = self.domain.geo_reference
809
810        # Call fit_interpolate.fit function
811        args = (points, )
[6197]812        kwargs = {'vertex_coordinates': None,
813                  'triangles': None,
814                  'mesh': self.domain.mesh,
[5897]815                  'point_attributes': values,
816                  'data_origin': data_georef.get_origin(),
817                  'mesh_origin': mesh_georef.get_origin(),
818                  'alpha': alpha,
819                  'verbose': verbose}
820
[6226]821        vertex_attributes = apply(fit_to_mesh, args, kwargs)
[5897]822
823        # Call underlying method using array values
[6228]824        self.set_values_from_array(vertex_attributes, location, indices,
825                                   use_cache=use_cache, verbose=verbose)
[5897]826
[6226]827    ##
828    # @brief Set quantity values from arbitray data points.
829    # @param points ??
830    # @param values ??
831    # @param alpha ??
832    # @param location ??
833    # @param indices ??
834    # @param data_georef ??
835    # @param verbose True if this method is to be verbose.
836    # @param use_cache ??
[7276]837    def set_values_from_points(self,
838                               points,
839                               values,
840                               alpha,
841                               location,
842                               indices,
843                               data_georef=None,
844                               verbose=False,
845                               use_cache=False):
[6226]846        """Set quantity values from arbitray data points using fit_interpolate.fit"""
[5897]847
848        raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
849
[6226]850    ##
851    # @brief Set quantity based on arbitrary points in a points file.
852    # @param filename Path to the points file.
853    # @param attribute_name
854    # @param alpha
855    # @param location
856    # @param indices
857    # @param verbose True if this method is to be verbose.
858    # @param use_cache
859    # @param max_read_lines
[7276]860    def set_values_from_file(self,
861                             filename,
862                             attribute_name,
863                             alpha,
864                             location,
865                             indices,
866                             verbose=False,
867                             use_cache=False,
868                             max_read_lines=None):
[6226]869        """Set quantity based on arbitrary points in a points file using
870        attribute_name selects name of attribute present in file.
[5897]871        If attribute_name is not specified, use first available attribute
[6226]872        as defined in geospatial_data.
[5897]873        """
874
875        from types import StringType
[6226]876
[5897]877        msg = 'Filename must be a text string'
878        assert type(filename) == StringType, msg
879
880        if location != 'vertices':
[6226]881            msg = "set_values_from_file is only defined for location='vertices'"
882            raise Exception, msg
[5897]883
[6244]884           
[6541]885        if True:
[5897]886            # Use mesh as defined by domain
[6226]887            # This used to cause problems for caching due to quantities
888            # changing, but it now works using the appropriate Mesh object.
[6541]889            # This addressed ticket:242 and was made to work when bug
890            # in ticket:314 was fixed 18 March 2009.
[5897]891            vertex_attributes = fit_to_mesh(filename,
[6226]892                                            mesh=self.domain.mesh,
[5897]893                                            alpha=alpha,
894                                            attribute_name=attribute_name,
895                                            use_cache=use_cache,
896                                            verbose=verbose,
897                                            max_read_lines=max_read_lines)
898        else:
899            # This variant will cause Mesh object to be recreated
[6226]900            # in fit_to_mesh thus doubling up on the neighbour structure
[6541]901            # FIXME(Ole): This is now obsolete 19 Jan 2009 except for bug
902            # (ticket:314) which was fixed 18 March 2009.
[5897]903            nodes = self.domain.get_nodes(absolute=True)
[6226]904            triangles = self.domain.get_triangles()
[5897]905            vertex_attributes = fit_to_mesh(filename,
[6226]906                                            nodes, triangles,
[5897]907                                            mesh=None,
908                                            alpha=alpha,
909                                            attribute_name=attribute_name,
910                                            use_cache=use_cache,
911                                            verbose=verbose,
912                                            max_read_lines=max_read_lines)
[6226]913
[5897]914        # Call underlying method using array values
[6228]915        if verbose:
[7317]916            log.critical('Applying fitted data to domain')
[6226]917        self.set_values_from_array(vertex_attributes, location,
[6228]918                                   indices, use_cache=use_cache,
919                                   verbose=verbose)
[5897]920
[6226]921    ##
922    # @brief Get index for maximum or minimum value of quantity.
923    # @param mode Either 'max' or 'min'.
924    # @param indices Set of IDs of elements to work on.
[5897]925    def get_extremum_index(self, mode=None, indices=None):
926        """Return index for maximum or minimum value of quantity (on centroids)
927
928        Optional arguments:
929            mode is either 'max'(default) or 'min'.
930            indices is the set of element ids that the operation applies to.
931
932        Usage:
933            i = get_extreme_index()
934
935        Notes:
936            We do not seek the extremum at vertices as each vertex can
937            have multiple values - one for each triangle sharing it.
938
939            If there are multiple cells with same maximum value, the
940            first cell encountered in the triangle array is returned.
941        """
942
943        V = self.get_values(location='centroids', indices=indices)
944
945        # Always return absolute indices
946        if mode is None or mode == 'max':
[6145]947            i = num.argmax(V)
[6226]948        elif mode == 'min':
[6145]949            i = num.argmin(V)
[6226]950        else:
951            raise ValueError, 'Bad mode value, got: %s' % str(mode)
[5897]952
953        if indices is None:
954            return i
955        else:
956            return indices[i]
957
[6226]958    ##
959    # @brief Get index for maximum value of quantity.
960    # @param indices Set of IDs of elements to work on.
[5897]961    def get_maximum_index(self, indices=None):
[6226]962        """See get extreme index for details"""
[5897]963
[6226]964        return self.get_extremum_index(mode='max', indices=indices)
[5897]965
[6226]966    ##
967    # @brief Return maximum value of quantity (on centroids).
968    # @param indices Set of IDs of elements to work on.
[5897]969    def get_maximum_value(self, indices=None):
970        """Return maximum value of quantity (on centroids)
971
972        Optional argument:
973            indices is the set of element ids that the operation applies to.
974
975        Usage:
976            v = get_maximum_value()
977
978        Note, we do not seek the maximum at vertices as each vertex can
[6226]979        have multiple values - one for each triangle sharing it
[5897]980        """
981
[6226]982        i = self.get_maximum_index(indices)
983        V = self.get_values(location='centroids')      #, indices=indices)
[5897]984
985        return V[i]
986
[6226]987    ##
988    # @brief Get location of maximum value of quantity (on centroids).
989    # @param indices Set of IDs of elements to work on.
[5897]990    def get_maximum_location(self, indices=None):
991        """Return location of maximum value of quantity (on centroids)
992
993        Optional argument:
994            indices is the set of element ids that the operation applies to.
995
996        Usage:
997            x, y = get_maximum_location()
998
999        Notes:
1000            We do not seek the maximum at vertices as each vertex can
1001            have multiple values - one for each triangle sharing it.
1002
1003            If there are multiple cells with same maximum value, the
[6226]1004            first cell encountered in the triangle array is returned.
[5897]1005        """
1006
1007        i = self.get_maximum_index(indices)
1008        x, y = self.domain.get_centroid_coordinates()[i]
1009
1010        return x, y
1011
[6226]1012    ##
1013    # @brief  Get index for minimum value of quantity.
1014    # @param indices Set of IDs of elements to work on.
[5897]1015    def get_minimum_index(self, indices=None):
[6226]1016        """See get extreme index for details"""
[5897]1017
[6226]1018        return self.get_extremum_index(mode='min', indices=indices)
[5897]1019
[6226]1020    ##
1021    # @brief Return minimum value of quantity (on centroids).
1022    # @param indices Set of IDs of elements to work on.
[5897]1023    def get_minimum_value(self, indices=None):
1024        """Return minimum value of quantity (on centroids)
1025
1026        Optional argument:
1027            indices is the set of element ids that the operation applies to.
1028
1029        Usage:
1030            v = get_minimum_value()
1031
[6226]1032        See get_maximum_value for more details.
[5897]1033        """
1034
1035        i = self.get_minimum_index(indices)
1036        V = self.get_values(location='centroids')
[6226]1037
[5897]1038        return V[i]
1039
[6226]1040
1041    ##
1042    # @brief Get location of minimum value of quantity (on centroids).
1043    # @param indices Set of IDs of elements to work on.
[5897]1044    def get_minimum_location(self, indices=None):
1045        """Return location of minimum value of quantity (on centroids)
1046
1047        Optional argument:
1048            indices is the set of element ids that the operation applies to.
1049
1050        Usage:
1051            x, y = get_minimum_location()
1052
1053        Notes:
1054            We do not seek the maximum at vertices as each vertex can
1055            have multiple values - one for each triangle sharing it.
1056
1057            If there are multiple cells with same maximum value, the
[6226]1058            first cell encountered in the triangle array is returned.
[5897]1059        """
1060
1061        i = self.get_minimum_index(indices)
1062        x, y = self.domain.get_centroid_coordinates()[i]
1063
1064        return x, y
1065
[6226]1066    ##
1067    # @brief Get values at interpolation points.
1068    # @param interpolation_points List of UTM coords or geospatial data object.
1069    # @param use_cache ??
1070    # @param verbose True if this method is to be verbose.
[7276]1071    def get_interpolated_values(self,
1072                                interpolation_points,
1073                                use_cache=False,
1074                                verbose=False):
[6226]1075        """Get values at interpolation points
[5897]1076
[6226]1077        The argument interpolation points must be given as either a
[5897]1078        list of absolute UTM coordinates or a geospatial data object.
1079        """
1080
1081        # FIXME (Ole): Points might be converted to coordinates relative to mesh origin
[6226]1082        # This could all be refactored using the
1083        # 'change_points_geo_ref' method of Class geo_reference.
[5897]1084        # The purpose is to make interpolation points relative
1085        # to the mesh origin.
1086        #
1087        # Speed is also a consideration here.
[6226]1088
1089        # Ensure that interpolation points is either a list of
[7276]1090        # points, Nx2 array, or geospatial and convert to numeric array
[6226]1091        if isinstance(interpolation_points, Geospatial_data):
[5897]1092            # Ensure interpolation points are in absolute UTM coordinates
[6226]1093            interpolation_points = \
1094                interpolation_points.get_data_points(absolute=True)
1095
[5897]1096        # Reconcile interpolation points with georeference of domain
[6226]1097        interpolation_points = \
1098            self.domain.geo_reference.get_relative(interpolation_points)
[5897]1099        interpolation_points = ensure_numeric(interpolation_points)
1100
[6226]1101
[5897]1102        # Get internal representation (disconnected) of vertex values
1103        vertex_values, triangles = self.get_vertex_values(xy=False,
[6226]1104                                                          smooth=False)
1105
[5897]1106        # Get possibly precomputed interpolation object
1107        I = self.domain.get_interpolation_object()
1108
[6226]1109        # Call interpolate method with interpolation points
[5897]1110        result = I.interpolate_block(vertex_values, interpolation_points,
[6226]1111                                     use_cache=use_cache, verbose=verbose)
1112
[5897]1113        return result
1114
[6226]1115    ##
1116    # @brief Get values as an array.
1117    # @param interpolation_points List of coords to get values at.
1118    # @param location Where to store results.
1119    # @param indices Set of IDs of elements to work on.
1120    # @param use_cache
1121    # @param verbose True if this method is to be verbose.
[7276]1122    def get_values(self,
[6654]1123                   interpolation_points=None,
1124                   location='vertices',
1125                   indices=None,
1126                   use_cache=False,
1127                   verbose=False):
[6226]1128        """Get values for quantity
[5897]1129
[7276]1130        Extract values for quantity as a numeric array.
[5897]1131
1132        Inputs:
1133           interpolation_points: List of x, y coordinates where value is
[6226]1134                                 sought (using interpolation). If points
1135                                 are given, values of location and indices
[5897]1136                                 are ignored. Assume either absolute UTM
1137                                 coordinates or geospatial data object.
[6226]1138
[5897]1139           location: Where values are to be stored.
1140                     Permissible options are: vertices, edges, centroids
1141                     and unique vertices. Default is 'vertices'
1142
1143
[5976]1144        The returned values will have the leading dimension equal to length of the indices list or
1145        N (all values) if indices is None.
[5897]1146
1147        In case of location == 'centroids' the dimension of returned
[7276]1148        values will be a list or a numerical array of length N, N being
[5897]1149        the number of elements.
[6226]1150
[5897]1151        In case of location == 'vertices' or 'edges' the dimension of
1152        returned values will be of dimension Nx3
1153
1154        In case of location == 'unique vertices' the average value at
1155        each vertex will be returned and the dimension of returned values
[6226]1156        will be a 1d array of length "number of vertices"
1157
[5897]1158        Indices is the set of element ids that the operation applies to.
1159
1160        The values will be stored in elements following their
1161        internal ordering.
1162        """
1163
1164        # FIXME (Ole): I reckon we should have the option of passing a
1165        #              polygon into get_values. The question becomes how
1166        #              resulting values should be ordered.
[6226]1167
[5897]1168        if verbose is True:
[7317]1169            log.critical('Getting values from %s' % location)
[5897]1170
1171        if interpolation_points is not None:
1172            return self.get_interpolated_values(interpolation_points,
1173                                                use_cache=use_cache,
1174                                                verbose=verbose)
[6226]1175
[5897]1176        # FIXME (Ole): Consider deprecating 'edges' - but not if it is used
[6226]1177        # elsewhere in ANUGA.
[5897]1178        # Edges have already been deprecated in set_values, see changeset:5521,
1179        # but *might* be useful in get_values. Any thoughts anyone?
[7044]1180        # YES (Ole): Edge values are necessary for volumetric balance check and inflow boundary. Keep them.
[6226]1181
1182        if location not in ['vertices', 'centroids',
1183                            'edges', 'unique vertices']:
1184            msg = 'Invalid location: %s' % location
[7276]1185            raise Exception, msg
[5897]1186
[6145]1187        import types
[6226]1188
[7276]1189        msg = "'indices' must be a list, array or None"
1190        assert isinstance(indices, (NoneType, list, num.ndarray)), msg
[5897]1191
1192        if location == 'centroids':
1193            if (indices ==  None):
1194                indices = range(len(self))
[7276]1195            return num.take(self.centroid_values, indices, axis=0)
[5897]1196        elif location == 'edges':
1197            if (indices ==  None):
1198                indices = range(len(self))
[7276]1199            return num.take(self.edge_values, indices, axis=0)
[5897]1200        elif location == 'unique vertices':
1201            if (indices ==  None):
[6191]1202                indices=range(self.domain.get_number_of_nodes())
[5897]1203            vert_values = []
1204
1205            # Go through list of unique vertices
1206            for unique_vert_id in indices:
1207                triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
[6226]1208
[5897]1209                # In case there are unused points
1210                if len(triangles) == 0:
1211                    msg = 'Unique vertex not associated with triangles'
[6226]1212                    raise Exception, msg
[5897]1213
1214                # Go through all triangle, vertex pairs
1215                # Average the values
1216                # FIXME (Ole): Should we merge this with get_vertex_values
1217                sum = 0
1218                for triangle_id, vertex_id in triangles:
1219                    sum += self.vertex_values[triangle_id, vertex_id]
[6226]1220                vert_values.append(sum / len(triangles))
[7276]1221            return num.array(vert_values, num.float)
[5897]1222        else:
1223            if (indices is None):
1224                indices = range(len(self))
[7276]1225            return num.take(self.vertex_values, indices, axis=0)
[5897]1226
[6226]1227    ##
1228    # @brief Set vertex values for all unique vertices based on array.
1229    # @param A Array to set values with.
1230    # @param indices Set of IDs of elements to work on.
[6228]1231    # @param use_cache ??
1232    # @param verbose??
[7276]1233    def set_vertex_values(self,
1234                          A,
1235                          indices=None,
1236                          use_cache=False,
1237                          verbose=False):
[5897]1238        """Set vertex values for all unique vertices based on input array A
[6226]1239        which has one entry per unique vertex, i.e. one value for each row in
1240        array self.domain.nodes.
[5897]1241
1242        indices is the list of vertex_id's that will be set.
1243
1244        This function is used by set_values_from_array
1245        """
1246
[6226]1247        # Check that A can be converted to array and is of appropriate dim
[7276]1248        A = ensure_numeric(A, num.float)
[5897]1249        assert len(A.shape) == 1
1250
1251        if indices is None:
1252            assert A.shape[0] == self.domain.get_nodes().shape[0]
1253            vertex_list = range(A.shape[0])
1254        else:
1255            assert A.shape[0] == len(indices)
1256            vertex_list = indices
1257
[6228]1258        #FIXME(Ole): This function ought to be faster.
1259        # We need to get the triangles_and_vertices list
1260        # from domain in one hit, then cache the computation of the
1261        # Nx3 array of vertex values that can then be assigned using
1262        # set_values_from_array.
1263        #
1264        # Alternatively, some C code would be handy
1265        #
1266        self._set_vertex_values(vertex_list, A)
1267           
1268    ##
1269    # @brief Go through list of unique vertices.
1270    # @param vertex_list ??
1271    # @param A ??
1272    def _set_vertex_values(self, vertex_list, A):
1273        """Go through list of unique vertices
1274        This is the common case e.g. when values
1275        are obtained from a pts file through fitting
1276        """
1277
[6226]1278        # Go through list of unique vertices
[5897]1279        for i_index, unique_vert_id in enumerate(vertex_list):
[6226]1280            triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
[5897]1281
1282            # In case there are unused points
[6226]1283            if len(triangles) == 0:
1284                continue
[5897]1285
1286            # Go through all triangle, vertex pairs
1287            # touching vertex unique_vert_id and set corresponding vertex value
1288            for triangle_id, vertex_id in triangles:
1289                self.vertex_values[triangle_id, vertex_id] = A[i_index]
1290
1291        # Intialise centroid and edge_values
1292        self.interpolate()
1293
[6226]1294    ##
1295    # @brief Smooth vertex values.
[6228]1296    def smooth_vertex_values(self, use_cache=False, verbose=False):
[6226]1297        """Smooths vertex values."""
[5897]1298
[6226]1299        A, V = self.get_vertex_values(xy=False, smooth=True)
[6228]1300        self.set_vertex_values(A, use_cache=use_cache, verbose=verbose)
[5897]1301
[6226]1302    ############################################################################
1303    # Methods for outputting model results
1304    ############################################################################
[5897]1305
[6226]1306    ##
1307    # @brief Get vertex values like an OBJ format i.e. one value per node.
1308    # @param xy True if we return X and Y as well as A and V.
1309    # @param smooth True if vertex values are to be smoothed.
1310    # @param precision The type of the result values (default float).
[6228]1311    def get_vertex_values(self, xy=True, smooth=None, precision=None):
[5897]1312        """Return vertex values like an OBJ format i.e. one value per node.
1313
1314        The vertex values are returned as one sequence in the 1D float array A.
1315        If requested the coordinates will be returned in 1D arrays X and Y.
1316
1317        The connectivity is represented as an integer array, V, of dimension
1318        Mx3, where M is the number of triangles. Each row has three indices
1319        defining the triangle and they correspond to elements in the arrays
[6226]1320        X, Y and A.
[5897]1321
[6226]1322        If smooth is True, vertex values corresponding to one common coordinate
1323        set will be smoothed by taking the average of vertex values for each
1324        node.  In this case vertex coordinates will be de-duplicated
1325        corresponding to the original nodes as obtained from the method
1326        general_mesh.get_nodes()
[5897]1327
[6226]1328        If no smoothings is required, vertex coordinates and values will be
1329        aggregated as a concatenation of values at vertices 0, vertices 1 and
1330        vertices 2.  This corresponds to the node coordinates obtained from the
1331        method general_mesh.get_vertex_coordinates()
[5897]1332
1333        Calling convention
1334        if xy is True:
[6226]1335           X, Y, A, V = get_vertex_values
[5897]1336        else:
[6226]1337           A, V = get_vertex_values
[5897]1338        """
1339
1340        if smooth is None:
1341            # Take default from domain
1342            try:
1343                smooth = self.domain.smooth
1344            except:
1345                smooth = False
1346
1347        if precision is None:
[7276]1348            precision = num.float
[5897]1349
1350        if smooth is True:
[6226]1351            # Ensure continuous vertex values by averaging values at each node
[5897]1352            V = self.domain.get_triangles()
1353            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
[7276]1354            A = num.zeros(N, num.float)
[6226]1355            points = self.domain.get_nodes()
1356
[5897]1357            if 1:
1358                # Fast C version
1359                average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1360                                      ensure_numeric(self.domain.number_of_triangles_per_node),
1361                                      ensure_numeric(self.vertex_values),
1362                                      A)
1363                A = A.astype(precision)
[6226]1364            else:
[5897]1365                # Slow Python version
1366                current_node = 0
1367                k = 0 # Track triangles touching on node
1368                total = 0.0
1369                for index in self.domain.vertex_value_indices:
1370                    if current_node == N:
[6226]1371                        msg = 'Current node exceeding number of nodes (%d) ' % N
[7276]1372                        raise Exception, msg
[5897]1373
1374                    k += 1
[6226]1375
[5897]1376                    volume_id = index / 3
1377                    vertex_id = index % 3
[6226]1378
[5897]1379                    v = self.vertex_values[volume_id, vertex_id]
1380                    total += v
1381
1382                    if self.domain.number_of_triangles_per_node[current_node] == k:
1383                        A[current_node] = total/k
[6226]1384
[5897]1385                        # Move on to next node
1386                        total = 0.0
1387                        k = 0
1388                        current_node += 1
1389        else:
[6226]1390            # Return disconnected internal vertex values
[5897]1391            V = self.domain.get_disconnected_triangles()
1392            points = self.domain.get_vertex_coordinates()
[7276]1393            A = self.vertex_values.flatten().astype(precision)
[5897]1394
[6226]1395        # Return
[5897]1396        if xy is True:
1397            X = points[:,0].astype(precision)
1398            Y = points[:,1].astype(precision)
[6226]1399
[5897]1400            return X, Y, A, V
1401        else:
[6226]1402            return A, V
[5897]1403
[6226]1404    ##
1405    # @brief Extrapolate conserved quantities from centroid.
[5897]1406    def extrapolate_first_order(self):
[6226]1407        """Extrapolate conserved quantities from centroid to vertices and edges
1408        for each volume using first order scheme.
[5897]1409        """
1410
1411        qc = self.centroid_values
1412        qv = self.vertex_values
1413        qe = self.edge_values
1414
1415        for i in range(3):
1416            qv[:,i] = qc
1417            qe[:,i] = qc
1418
1419        self.x_gradient *= 0.0
1420        self.y_gradient *= 0.0
1421
[6226]1422    ##
1423    # @brief Compute the integral of quantity across entire domain.
1424    # @return The integral.
1425    def get_integral(self):
1426        """Compute the integral of quantity across entire domain."""
[5897]1427
[6191]1428        areas = self.domain.get_areas()
[5897]1429        integral = 0
1430        for k in range(len(self.domain)):
[6191]1431            area = areas[k]
[5897]1432            qc = self.centroid_values[k]
1433            integral += qc*area
1434
1435        return integral
1436
[6226]1437    ##
1438    # @brief get the gradients.
[5897]1439    def get_gradients(self):
[6226]1440        """Provide gradients. Use compute_gradients first."""
[5897]1441
1442        return self.x_gradient, self.y_gradient
1443
[6226]1444    ##
1445    # @brief ??
1446    # @param timestep ??
[5897]1447    def update(self, timestep):
1448        # Call correct module function
1449        # (either from this module or C-extension)
1450        return update(self, timestep)
1451
[6226]1452    ##
1453    # @brief ??
[5897]1454    def compute_gradients(self):
1455        # Call correct module function
1456        # (either from this module or C-extension)
1457        return compute_gradients(self)
1458
[6226]1459    ##
1460    # @brief ??
[5897]1461    def limit(self):
1462        # Call correct module depending on whether
1463        # basing limit calculations on edges or vertices
1464        limit_old(self)
1465
[6226]1466    ##
1467    # @brief ??
[5897]1468    def limit_vertices_by_all_neighbours(self):
1469        # Call correct module function
1470        # (either from this module or C-extension)
1471        limit_vertices_by_all_neighbours(self)
1472
[6226]1473    ##
1474    # @brief ??
[5897]1475    def limit_edges_by_all_neighbours(self):
1476        # Call correct module function
1477        # (either from this module or C-extension)
1478        limit_edges_by_all_neighbours(self)
1479
[6226]1480    ##
1481    # @brief ??
[5897]1482    def limit_edges_by_neighbour(self):
1483        # Call correct module function
1484        # (either from this module or C-extension)
[6226]1485        limit_edges_by_neighbour(self)
[5897]1486
[6226]1487    ##
1488    # @brief ??
[5897]1489    def extrapolate_second_order(self):
1490        # Call correct module function
1491        # (either from this module or C-extension)
1492        compute_gradients(self)
1493        extrapolate_from_gradient(self)
[6226]1494
1495    ##
1496    # @brief ??
[5897]1497    def extrapolate_second_order_and_limit_by_edge(self):
1498        # Call correct module function
1499        # (either from this module or C-extension)
1500        extrapolate_second_order_and_limit_by_edge(self)
1501
[6226]1502    ##
1503    # @brief ??
[5897]1504    def extrapolate_second_order_and_limit_by_vertex(self):
1505        # Call correct module function
1506        # (either from this module or C-extension)
1507        extrapolate_second_order_and_limit_by_vertex(self)
1508
[6226]1509    ##
1510    # @brief ??
1511    # @param bound ??
[5897]1512    def bound_vertices_below_by_constant(self, bound):
1513        # Call correct module function
1514        # (either from this module or C-extension)
1515        bound_vertices_below_by_constant(self, bound)
1516
[6226]1517    ##
1518    # @brief ??
1519    # @param quantity ??
[5897]1520    def bound_vertices_below_by_quantity(self, quantity):
1521        # Call correct module function
1522        # (either from this module or C-extension)
1523
1524        # check consistency
1525        assert self.domain == quantity.domain
[6226]1526        bound_vertices_below_by_quantity(self, quantity)
[5897]1527
[6226]1528    ##
1529    # @brief ??
[5897]1530    def backup_centroid_values(self):
1531        # Call correct module function
1532        # (either from this module or C-extension)
1533        backup_centroid_values(self)
1534
[6226]1535    ##
1536    # @brief ??
1537    # @param a ??
1538    # @param b ??
1539    def saxpy_centroid_values(self, a, b):
[5897]1540        # Call correct module function
1541        # (either from this module or C-extension)
[6226]1542        saxpy_centroid_values(self, a, b)
[5897]1543
[6226]1544
1545##
1546# @brief OBSOLETE!
[5897]1547class Conserved_quantity(Quantity):
[6226]1548    """Class conserved quantity being removed, use Quantity."""
[5897]1549
1550    def __init__(self, domain, vertex_values=None):
1551        msg = 'ERROR: Use Quantity instead of Conserved_quantity'
1552        raise Exception, msg
1553
1554
[6226]1555######
1556# Prepare the C extensions.
1557######
[5897]1558
1559from anuga.utilities import compile
1560
[6226]1561if compile.can_use_C_extension('quantity_ext.c'):
1562    # Underlying C implementations can be accessed
1563
[5897]1564    from quantity_ext import \
1565         average_vertex_values,\
1566         backup_centroid_values,\
1567         saxpy_centroid_values,\
1568         compute_gradients,\
1569         limit_old,\
1570         limit_vertices_by_all_neighbours,\
1571         limit_edges_by_all_neighbours,\
1572         limit_edges_by_neighbour,\
1573         limit_gradient_by_neighbour,\
1574         extrapolate_from_gradient,\
1575         extrapolate_second_order_and_limit_by_edge,\
1576         extrapolate_second_order_and_limit_by_vertex,\
1577         bound_vertices_below_by_constant,\
1578         bound_vertices_below_by_quantity,\
1579         interpolate_from_vertices_to_edges,\
1580         interpolate_from_edges_to_vertices,\
[6226]1581         update
[5897]1582else:
[6226]1583    msg = 'C implementations could not be accessed by %s.\n ' % __file__
[5897]1584    msg += 'Make sure compile_all.py has been run as described in '
1585    msg += 'the ANUGA installation guide.'
1586    raise Exception, msg
Note: See TracBrowser for help on using the repository browser.