Changeset 459

Timestamp:

Oct 28, 2004, 4:20:00 PM (20 years ago)

Author:

duncan

Message:

adding region concept with quantities. Redundant functions.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• domain.py (modified) (1 diff)
• quantity.py (modified) (5 diffs)
• test_domain.py (modified) (2 diffs)
• test_quantity.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/domain.py

@@ -209 +209 @@
                 raise msg
 
- 
+    def set_region(self, functions):
+        for tag in self.tagged_elements.keys():
+            for function in functions:
+                print "tag",tag
+                print "function",function
+                function(tag, self.tagged_elements[tag], self)
+
+                #Do we need to do this sort of thing?
+                #self = function(tag, self.tagged_elements[tag], self)
+            
     #MISC        
     def check_integrity(self):

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -154 +154 @@
 
             
-    def set_array_values(self, values, location='vertices'):
+    def set_array_values(self, values, location='vertices', indexes = None):
         """Set values for quantity
 
@@ -175 +175 @@
         """
 
-        from Numeric import array, Float
+        from Numeric import array, Float, Int
 
         values = array(values).astype(Float)
 
+        if (indexes <> None):
+            indexes = array(indexes).astype(Int)
+            msg = 'Number of values must match number of indexes'
+            assert values.shape[0] == indexes.shape[0], msg 
+            
         N = self.centroid_values.shape[0]
         
@@ -184 +189 @@
             assert len(values.shape) == 1, 'Values array must be 1d'
 
-            msg = 'Number of values must match number of elements'
-            assert values.shape[0] == N, msg
-            
-            self.centroid_values = values                        
+            if indexes == None:
+                msg = 'Number of values must match number of elements'
+                assert values.shape[0] == N, msg
+            
+                self.centroid_values = values
+            else:
+                msg = 'Number of values must match number of indexes'
+                assert values.shape[0] == indexes.shape[0], msg
+                
+                #Brute force
+                for i in range(len(indexes)):
+                    self.centroid_values[indexes[i]] = values[i]
+                    
         elif location == 'edges':
             assert len(values.shape) == 2, 'Values array must be 2d'
@@ -200 +214 @@
         else:
             if len(values.shape) == 1:
-                #Values are being specified once for each unique vertex
-                msg = 'Number of values must match number of vertices'
-                assert values.shape[0] == self.domain.coordinates.shape[0], msg
-
-                self.set_vertex_values(values)
+
+                if indexes == None:
+                    #Values are being specified once for each unique vertex
+                    msg = 'Number of values must match number of vertices'
+                    assert values.shape[0] == self.domain.coordinates.shape[0], msg
+                    self.set_vertex_values(values)
+                else:
+                    for element_index, value in map(None, indexes, values): 
+                        self.vertex_values[element_index, 0] = value
+                        self.vertex_values[element_index, 1] = value
+                        self.vertex_values[element_index, 2] = value
             elif len(values.shape) == 2:
                 #Vertex values are given as a triplet for each triangle
-                msg = 'Number of values must match number of elements'
-                assert values.shape[0] == N, msg
                 
                 msg = 'Array must be N x 3'            
                 assert values.shape[1] == 3, msg
-                self.vertex_values = values            
+                
+                if indexes == None:
+                    self.vertex_values = values
+                else:
+                    for element_index, value in map(None, indexes, values): 
+                        self.vertex_values[element_index] = value
             else:    
                 msg = 'Values array must be 1d or 2d'
                 raise msg
+            
+ #FIXME - roll this back into set_array_values, if possible             
+    def set_array_values_by_index(self, values, location='vertices', indexes = None):
+        from Numeric import array, Float, Int
+
+        values = array(values).astype(Float)
+        
+        if not (indexes == None):
+            indexes = array(indexes).astype(Int)
+
+        N = self.centroid_values.shape[0]
+        
+        if location == 'centroids':
+            assert len(values.shape) == 1, 'Values array must be 1d'
+
+            if indexes == None:
+                msg = 'Number of values must match number of elements'
+                assert values.shape[0] == N, msg
+            
+                self.centroid_values = values
+            else:
+                msg = 'Number of values must match number of indexes'
+                assert values.shape[0] == indexes.shape[0], msg
+                
+                #Brute force
+                for i in range(len(indexes)):
+                    self.centroid_values[indexes[i]] = values[i]
                 
 
@@ -246 +296 @@
             for triangle_id, vertex_id in L:
                 self.vertex_values[triangle_id, vertex_id] = A[k]
+                
+        #Intialise centroid and edge_values
+        self.interpolate()
+ 
+    def set_some_vertex_values(self, A, indexes):
+        """Set vertex values for some triangles based on input array A
+        which is assumed to have one entry per (unique) vertex, i.e.
+        one value for each row in array self.domain.coordinates.
+        """
+
+        from Numeric import array, Float
+        
+        #Assert that A can be converted to a Numeric array of appropriate dim
+        A = array(A, Float)
+        msg = 'Number of values must match number of indexes'
+        assert A.shape[0] == len(indexes), msg 
+        if (len(A.shape) == 2):
+            assert A.shape[0] == len(indexes), msg  
+            msg = 'Each vertex must have a list of three values, or one value'
+            assert A.shape[1] == 3, msg  
+        elif (len(A.shape) <> 1):  
+            msg = 'vertex values must be 1D or 2D.'
+            raise msg 
+        
+        #Go through list of triangles
+        #FIXME This is brute force.  Can it be optimised?
+        for element_index, value in map(None, indexes, A):            
+            if len(A.shape) == 2:
+                self.vertex_values[element_index] = value
+            else:
+                self.vertex_values[element_index, 0] = value
+                self.vertex_values[element_index, 1] = value
+                self.vertex_values[element_index, 2] = value
                 
         #Intialise centroid and edge_values

inundation/ga/storm_surge/pyvolution/test_domain.py

@@ -7 +7 @@
 from config import epsilon
 from Numeric import allclose, array
+
+
+def add_to_verts(tag, elements, domain):
+    if tag == "mound":
+        print "Selected"
+        
+
 
 class TestCase(unittest.TestCase):
@@ -267 +274 @@
         #    assert allclose(domain.quantities[name].centroid_values, x)   
 
-        
+
+    def testz_set_region(self):
+        """Domain implements a default first order gradient limiter
+        """
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0,0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+        boundary = { (0, 0): 'Third',
+                     (0, 2): 'First',
+                     (2, 0): 'Second',
+                     (2, 1): 'Second',
+                     (3, 1): 'Second',
+                     (3, 2): 'Third'}             
+
+        domain = Domain(points, vertices, boundary,
+                        conserved_quantities =\
+                        ['level', 'xmomentum', 'ymomentum'])                
+        domain.check_integrity()
+
+        domain.set_quantity('level', [[1,2,3], [5,5,5],
+                                      [0,0,9], [-6, 3, 3]])
+
+        assert allclose( domain.quantities['level'].centroid_values,
+                         [2,5,3,0] )
+                         
+        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
+                                          [3,3,3], [4, 4, 4]])
+
+        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
+                                          [30,30,30], [40, 40, 40]])        
+
+
+        domain.distribute_to_vertices_and_edges()
+
+        #First order extrapolation
+        assert allclose( domain.quantities['level'].vertex_values,
+                         [[ 2.,  2.,  2.],
+                          [ 5.,  5.,  5.],
+                          [ 3.,  3.,  3.],
+                          [ 0.,  0.,  0.]])
+
+        domain.build_tagged_elements_dictionary({'mound':[0,1]})
+        domain.set_region([add_to_verts])
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(TestCase,'testz')
     runner = unittest.TextTestRunner()
     runner.run(suite)

inundation/ga/storm_surge/pyvolution/test_quantity.py

@@ -611 +611 @@
         
         
-                                 
-
-        
+ 
+    def set_array_values_by_index(self):
+
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        from Numeric import zeros, Float
+        
+        #Create basic mesh
+        points, vertices, boundary = rectangular(1, 1)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        #print "domain.number_of_elements ",domain.number_of_elements 
+        quantity = Quantity(domain,[[1,1,1],[2,2,2]])
+        value = [7]
+        indexes = [1]
+        quantity.set_array_values_by_index(value,
+                                           location = 'centroids',
+                                           indexes = indexes)
+        #print "quantity.centroid_values",quantity.centroid_values
+        
+        assert allclose(quantity.centroid_values, [1,7])
+
+        quantity.set_array_values([15,20,25], indexes = indexes)
+        assert allclose(quantity.centroid_values, [1,20])
+
+        quantity.set_array_values([15,20,25], indexes = indexes)
+        assert allclose(quantity.centroid_values, [1,20])
+        
+    def test_setting_some_vertex_values(self):
+
+        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)
+        #print "domain.number_of_elements ",domain.number_of_elements 
+        quantity = Quantity(domain,[[1,1,1],[2,2,2],[3,3,3],
+                                    [4,4,4],[5,5,5],[6,6,6]])
+        value = [7]
+        indexes = [1]
+        quantity.set_array_values(value,
+                                  location = 'centroids',
+                                  indexes = indexes)
+        #print "quantity.centroid_values",quantity.centroid_values
+        assert allclose(quantity.centroid_values, [1,7,3,4,5,6])
+        
+        value = [[15,20,25]]
+        quantity.set_array_values(value, indexes = indexes)
+        #print "1 quantity.vertex_values",quantity.vertex_values
+        assert allclose(quantity.vertex_values[1], value[0])
+
+        values = [10,100,50]
+        quantity.set_array_values(values, indexes = [0,1,5])
+        #print "2 quantity.vertex_values",quantity.vertex_values
+        assert allclose(quantity.vertex_values[0], [10,10,10])
+        assert allclose(quantity.vertex_values[5], [50,50,50])
+        quantity.interpolate()
+        #print "quantity.centroid_values",quantity.centroid_values
+        assert allclose(quantity.centroid_values, [10,100,3,4,5,50])
+       
+        values = [[31,30,29],[400,400,400],[1000,999,998]]
+        quantity.set_array_values(values, indexes = [3,3,5])
+        quantity.interpolate()
+        assert allclose(quantity.centroid_values, [10,100,3,400,5,999])
 #-------------------------------------------------------------
 if __name__ == "__main__":