Changeset 7793 for anuga_work/development/anuga_1d/quantity.py

Timestamp:

Jun 6, 2010, 10:48:41 PM (13 years ago)

Author:

steve

Message:

Setting up test function for Quantity

File:

• anuga_work/development/anuga_1d/quantity.py (modified) (5 diffs)

anuga_work/development/anuga_1d/quantity.py

@@ -1 @@
-"""Class Quantity - Implements values at each 1d element
+"""
+Class Quantity - Implements values at each 1d element
 
 To create:
@@ -85 +86 @@
         """
         return self.centroid_values.shape[0]
-    
+
+    def __neg__(self):
+        """Negate all values in this quantity giving meaning to the
+        expression -Q where Q is an instance of class Quantity
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(-self.vertex_values)
+        return Q
+
+
+
+    def __add__(self, other):
+        """Add to self anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        Q = Quantity(self.domain)
+        Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
+        return result
+
+    def __radd__(self, other):
+        """Handle cases like 7+Q, where Q is an instance of class Quantity
+        """
+
+        return self + other
+
+    def __sub__(self, other):
+        return self + -other            # Invoke self.__neg__()
+
+    def __mul__(self, other):
+        """Multiply self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The product of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be products from q1 and q2
+        result.vertex_values = self.vertex_values * Q.vertex_values
+        result.centroid_values = self.centroid_values * Q.centroid_values
+
+        return result
+
+    def __rmul__(self, other):
+        """Handle cases like 3*Q, where Q is an instance of class Quantity
+        """
+
+        return self * other
+
+    def __div__(self, other):
+        """Divide self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values_from_numeric for details
+
+        Zero division is dealt with by adding an epsilon to the divisore
+        FIXME (Ole): Replace this with native INF once we migrate to NumPy
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The quotient of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be quotient of q1 and q2
+        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
+        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
+
+        return result
+
+    def __rdiv__(self, other):
+        """Handle cases like 3/Q, where Q is an instance of class Quantity
+        """
+
+        return self / other
+
+    def __pow__(self, other):
+        """Raise quantity to (numerical) power
+
+        As with __mul__ vertex values are processed entry by entry
+        while centroid and edge values are re-interpolated.
+
+        Example using __pow__:
+          Q = (Q1**2 + Q2**2)**0.5
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:
+            Q = Quantity(self.domain)
+            Q.set_values_from_numeric(other)
+
+        result = Quantity(self.domain)
+
+        # The power of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values_from_numeric (which calls interpolate). Otherwise
+        # centroid values wouldn't be correct
+        result.vertex_values = self.vertex_values ** other
+        result.centroid_values = self.centroid_values ** other
+
+        return result
+
+    def set_values_from_numeric(self, numeric):
+        import Numeric
+        x = Numeric.array([1.0,2.0])
+        y = [1.0,2.0]
+
+        if type(numeric) == type(y):
+            self.set_values_from_array(numeric)
+        elif type(numeric) == type(x):
+            self.set_values_from_array(numeric)
+        elif callable(numeric):
+            self.set_values_from_function(numeric)
+        elif isinstance(numeric, Quantity):
+            self.set_values_from_quantity(numeric)
+        else:   # see if it's coercible to a float (float, int or long, etc)
+            try:
+                numeric = float(numeric)
+            except ValueError:
+                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
+                raise Exception, msg
+            self.set_values_from_constant(numeric)
+
+    def set_values_from_constant(self,numeric):
+
+        self.vertex_values[:]   = numeric
+        self.centroid_values[:] = numeric
+
+
+    def set_values_from_array(self,numeric):
+
+        self.vertex_values[:]   = numeric
+        self.interpolate()
+
+
+    def set_values_from_quantity(self,quantity):
+
+        self.vertex_values[:]   = quantity.vertex_values
+        self.centroid_values[:] = quantity.centroid_values
+
+    def set_values_from_function(self,function):
+
+        self.vertex_values[:]   = map(function, self.domain.vertices)
+        self.centroid_values[:] = map(function, self.domain.centroids)
+
+
     def interpolate(self):
         """
@@ -218 +390 @@
             assert values.shape[1] == 2, msg
 
-            self.vertex_values = values
+            self.vertex_values[:] = values
 
 
@@ -428 +600 @@
         """
 
-        #print 'compute_gradient'
 
         from Numeric import array, zeros, Float
@@ -908 +1079 @@
         self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
         
-class Conserved_quantity(Quantity):
-    """Class conserved quantity adds to Quantity:
-
-    storage and method for updating, and
-    methods for extrapolation from centropid to vertices inluding
-    gradients and limiters
-    """
-
-    def __init__(self, domain, vertex_values=None):
-        Quantity.__init__(self, domain, vertex_values)
-
-        print "Use Quantity instead of Conserved_quantity"
-
-
-##
-##def newLinePlot(title='Simple Plot'):
-##    import Gnuplot
-##    g = Gnuplot.Gnuplot()
-##    g.title(title)
-##    g('set data style linespoints') 
-##    g.xlabel('x')
-##    g.ylabel('y')
-##    return g
-##
-##def linePlot(g,x,y):
-##    import Gnuplot
-##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
+
 
 def newLinePlot(title='Simple Plot'):