source: inundation/pyvolution/quantity.py @ 2267

Last change on this file since 2267 was 2262, checked in by ole, 19 years ago

Fixed missing geo reference in set_quantity and added tests

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