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

Last change on this file since 4652 was 4652, checked in by duncan, 17 years ago

changing max points per cell

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