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

Last change on this file since 4815 was 4815, checked in by ole, 16 years ago

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