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

Last change on this file since 4127 was 4127, checked in by ole, 17 years ago

Fixed ticket:121 - it was a simple case of commenting out an old diagnostics output that relied on variables that were in existence before the change to geospatial_data objects. The okushiri validation now runs past the point where it failed in changeset:4107

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