Changeset 4592

Timestamp:

Jul 5, 2007, 9:43:02 AM (18 years ago)

Author:

ole

Message:

Allowed set_quantity restriction by polygon to observe geo_referencing.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• neighbour_mesh.py (modified) (1 diff)
• quantity.py (modified) (1 diff)
• test_quantity.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -769 +769 @@
         return self.lone_vertices
 
-    def get_centroid_coordinates(self):
+    def get_centroid_coordinates(self, absolute=False):
         """Return all centroid coordinates.
         Return all centroid coordinates for all triangles as an Nx2 array
         (ordered as x0, y0 for each triangle)
-        """
-        return self.centroid_coordinates
-        
-        
+
+        Boolean keyword argument absolute determines whether coordinates
+        are to be made absolute by taking georeference into account
+        Default is False as many parts of ANUGA expects relative coordinates.
+        """
+
+        V = self.centroid_coordinates
+        if absolute is True:
+            if not self.geo_reference.is_absolute():
+                V = self.geo_reference.get_absolute(V)
+            
+        return V
+
+        
     def get_radii(self):
         """Return all radii.

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -317 +317 @@
 
 
-            points = self.domain.get_centroid_coordinates()
+            points = self.domain.get_centroid_coordinates(absolute=True)
             indices = inside_polygon(points, polygon)
             

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -433 +433 @@
             raise msg
 
-            
+
+
+
+
+    def test_set_vertex_values_using_general_interface_subset_and_geo(self):
+        """test_set_vertex_values_using_general_interface_with_subset(self):
+        Test that indices and polygon works using georeferencing
+        """
+        
+        quantity = Quantity(self.mesh4)
+        G = Geo_reference(56, 10, 100)
+        quantity.domain.geo_reference = G
+
+        #print quantity.domain.get_nodes(absolute=True)
+
+
+        # 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]])        
+        
+
+
+        # Now try with polygon (pick points where y>2)
+        polygon = array([[0,2.1], [4,2.1], [4,7], [0,7]])
+        polygon += [G.xllcorner, G.yllcorner]
+        
+        quantity.set_values(0.0)
+        quantity.set_values(3.14, polygon=polygon)
+        
+        assert allclose(quantity.vertex_values,
+                        [[0,0,0], [0,0,0], [0,0,0],
+                         [3.14,3.14,3.14]])                
+
+
+        # Another polygon (pick triangle 1 and 2 (rightmost triangles)
+        polygon = array([[2.1, 0.0], [3.5,0.1], [2,2.2], [0.2,2]])
+        polygon += [G.xllcorner, G.yllcorner]
+        
+        quantity.set_values(0.0)
+        quantity.set_values(3.14, polygon=polygon)
+
+        assert allclose(quantity.vertex_values,
+                        [[0,0,0],
+                         [3.14,3.14,3.14],
+                         [3.14,3.14,3.14],                         
+                         [0,0,0]])                
+
+
 
     def test_set_values_using_fit(self):