Changeset 4579

Timestamp:

Jul 3, 2007, 1:30:49 PM (18 years ago)

Author:

ole

Message:

Thoughts towards allowing 'polygon' to be specified in set_quantity.
This is far from done.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (1 diff)
• utilities/polygon.py (modified) (4 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -179 +179 @@
                    alpha = None,
                    location = 'vertices',
+                   polygon = None,
                    indices = None,
                    verbose = False,
@@ -262 +263 @@
                   will be left undefined.
 
+
+        polygon: Restrict update of quantity to locations that fall
+                 inside polygon. Polygon works by selecting indices
+                 and calling set_values recursively.
+
+        indices: Restrict update of quantity to locations that are 
+                 identified by indices (e.g. node ids if location
+                 is 'vertices')        
+        
         verbose: True means that output to stdout is generated
 
@@ -278 +288 @@
         from types import FloatType, IntType, LongType, ListType, NoneType
         from Numeric import ArrayType
+
+
+        # Polygon situation
+        #if polygon is not None:
+        #    if indices is not None:
+        #        msg = 'Only one of polygon and indices can be specified'
+        #        raise Exception, msg
+        #
+        #    Need to get candidate points. I think we should
+        #    simplify this whole thing a bit. Do we really need location?
+        #    point_indices = inside_polygon(points, polygon)
+
+
 
         #General input checks
@@ -388 +411 @@
         """
 
+        # FIXME (Ole): Somehow indices refer to centroids
+        # rather than vertices as default. See unit test
+        # test_set_vertex_values_using_general_interface_with_subset(self):
+        
 
         if location == 'centroids':

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -377 +377 @@
 
 
+    def test_set_vertex_values_using_general_interface_with_subset(self):
+        """test_set_vertex_values_using_general_interface_with_subset(self):
+        Test that indices and polygon works
+        """
+        
+        quantity = Quantity(self.mesh4)
+
+
+        quantity.set_values([0,2,3,5], indices=[0,2,3,5])
+        assert allclose(quantity.vertex_values,
+                        [[0,0,2], [0,2,0], [0,2,5], [3,0,0]])
+
+
+        # Constant
+        quantity.set_values(0.0)
+        quantity.set_values(3.14, indices=[0,2], location='vertices')
+
+        # Indices refer to triangle numbers here - not vertices (why?)
+        assert allclose(quantity.vertex_values,
+                        [[3.14,3.14,3.14], [0,0,0],
+                         [3.14,3.14,3.14], [0,0,0]])        
+        
+
+
+        # FIXME: Not done yet
+        # Now try with polygon (pick points where y>2)
+        #polygon = [[0,2.1], [4,2.1], [4,7], [0,7]]
+        #quantity.set_values(0.0)
+        #quantity.set_values(3.14, polygon=polygon, location='vertices')
+        #
+        #print quantity.vertex_values
+        #
+        ## Indices refer to triangle numbers here - not vertices (why?)
+        #assert allclose(quantity.vertex_values,
+        #                [[0,0,0], [0,0,0], [0,0,0],
+        #                 [3.14,3.14,3.14]])                
+
+        
 
 

anuga_core/source/anuga/utilities/polygon.py

@@ -622 +622 @@
         return z
 
-def read_polygon(filename,split=','):
+
+def read_polygon(filename, split=','):
     """Read points assumed to form a polygon.
        There must be exactly two numbers in each line separated by a comma.
@@ -639 +640 @@
     return polygon
 
-def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
+
+def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
     """Populate given polygon with uniformly distributed points.
 
@@ -694 +696 @@
 
     return points
+
 
 def point_in_polygon(polygon, delta=1e-8):
@@ -739 +742 @@
     return point
 
+
 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
-    """Calcalutes the approximate number of triangles inside the bounding polygon
-    and the other interior regions
+    """Calcalutes the approximate number of triangles inside the
+    bounding polygon and the other interior regions
 
     Polygon areas are converted to square Kms