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

Last change on this file since 6541 was 6541, checked in by ole, 15 years ago

Implemented fix to ticket:314.

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