source: branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py @ 6304

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

Initial commit of numpy changes. Still a long way to go.

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