Changeset 5208

Timestamp:

Apr 11, 2008, 4:20:45 PM (16 years ago)

Author:

duncan

Message:

Add averaging to set region.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• region.py (modified) (6 diffs)
• test_region.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/region.py

@@ -1 @@
-"""boundary.py - Classes for implementing region conditions
+"""region.py - Classes for implementing region conditions
 
 NOTE: I've only set it up for testing X values of constants and functions.
@@ -8 +8 @@
 # FIXME (DSG-DSG) add better comments
 
+from Numeric import average
 class Region:
     """Base class for modifying quantities based on a region.
@@ -61 +62 @@
     """
     Will add a value to the current quantity value.
+    
+    quantity = initial_quantity + X.
+    
+    This method does not work with functions
+    
+    Use average to add X to a mean average of the quantity values.
+    eg if applying this to elevation this will give a flat roof.
     """
     
-    def __init__(self, tag, quantity, X, location='vertices', initial_quantity=None):
+    def __init__(self, tag, quantity, X, location='vertices', initial_quantity=None, average=False):
         #I have to get this going!
         #Region.__init__(self)
@@ -70 +78 @@
         self.location = location
         self.X = X
+        self.average = average
         if initial_quantity is None:
             self.quantity_initial_value = quantity
@@ -87 +96 @@
             #              indices=self.build_indices(elements, domain),
             #              location=self.location) + self.X
-
             Q = domain.get_quantity(self.quantity_initial_value)
-            new_values = Q.get_values(indices=self.build_indices(elements, domain),
-                                      location=self.location) + self.X            
+            if self.average is True:
+                # Average the points, and then the verts in the averaged point.  
+                values = Q.get_values(indices=self.build_indices(elements, domain),
+                                      location=self.location)
+                av = average(values)
+                if self.location == "vertices":
+                    av = average(av)
+                new_values = av + self.X    
+            else:
+                new_values = Q.get_values(indices=self.build_indices(elements, domain),
+                                      location=self.location) + self.X    
             domain.set_quantity(self.quantity_answer, new_values,
                                 indices=self.build_indices(elements, domain),
@@ -97 +114 @@
 class Add_quantities(Region):
     """
-    Will add a value to the current quantity value.
+    Will add a quantity to the current quantity value.
     """
     

anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py

@@ -7 +7 @@
 from region import *
 #from anuga.config import epsilon
-from Numeric import allclose #, array, ones, Float
+from Numeric import allclose, average #, array, ones, Float
 
 
@@ -174 +174 @@
                          [ 0.07,  0.07,  0.07],
                          [ 0.07,  0.07,  0.07]])
+                         
+    def test_unique_vertices_average_loc_vert(self):
+        """
+        get values based on triangle lists.
+        """
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.build_tagged_elements_dictionary({'bottom':[0,1],
+                                                 'top':[4,5],
+                                                 'not_bottom':[2,3,4,5]})
+
+        #Set friction
+        domain.set_quantity('friction', add_x_y)
+        av_bottom = 2.0/3.0
+        add = 60.0
+        calc_frict = av_bottom + add
+        domain.set_region(Add_value_to_region('bottom', 'friction', add,
+                          initial_quantity='friction',
+                           #location='unique vertices',
+                           location='vertices',
+                           average=True
+                          ))
+
+        #print domain.quantities['friction'].get_values()
+        frict_points = domain.quantities['friction'].get_values()
+        assert allclose(frict_points[0],\
+                        [ calc_frict, calc_frict, calc_frict])
+        assert allclose(frict_points[1],\
+                        [ calc_frict, calc_frict, calc_frict])
+  
+    def test_unique_vertices_average_loc_unique_vert(self):
+        """
+        get values based on triangle lists.
+        """
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 3)
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.build_tagged_elements_dictionary({'bottom':[0,1],
+                                                 'top':[4,5],
+                                                 'not_bottom':[2,3,4,5]})
+
+        #Set friction
+        domain.set_quantity('friction', add_x_y)
+        av_bottom = 2.0/3.0
+        add = 60.0
+        calc_frict = av_bottom + add
+        domain.set_region(Add_value_to_region('bottom', 'friction', add,
+                          initial_quantity='friction',
+                           location='unique vertices',
+                           average=True
+                          ))
+
+        #print domain.quantities['friction'].get_values()
+        frict_points = domain.quantities['friction'].get_values()
+        assert allclose(frict_points[0],\
+                        [ calc_frict, calc_frict, calc_frict])
+        assert allclose(frict_points[1],\
+                        [ calc_frict, calc_frict, calc_frict])
+        assert allclose(frict_points[2],\
+                        [ calc_frict, 1.0 + 2.0/3.0, calc_frict])
+        assert allclose(frict_points[3],\
+                        [ 2.0/3.0,calc_frict, 1.0 + 2.0/3.0])
+                                                
+                         
 #-------------------------------------------------------------
 if __name__ == "__main__":