Changeset 4836

Timestamp:

Nov 20, 2007, 5:38:57 PM (17 years ago)

Author:

ole

Message:

Arranged for timestepping statistic for chosen triangle, e.g. one of the
gauges in the Okushiri example.
The underlying function is currently brute force, but OK for now.

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py (modified) (3 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py (modified) (2 diffs)
• anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• anuga_validation/okushiri_2005/run_okushiri.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -637 +637 @@
 
 
-    def timestepping_statistics(self, track_speeds=False):
+    def timestepping_statistics(self,
+                                track_speeds=False,
+                                triangle_id=None):
         """Return string with time stepping statistics for printing or logging
 
         Optional boolean keyword track_speeds decides whether to report
         location of smallest timestep as well as a histogram and percentile
-        report. 
+        report.
+
+        Optional keyword triangle_id can be used to specify a particular
+        triangle rather than the one with the largest speed. 
         """
 
@@ -716 +721 @@
             
             # Find index of largest computed flux speed
-            k = self.k = argmax(self.max_speed)
+            if triangle_id is None:
+                k = self.k = argmax(self.max_speed)
+            else:
+                errmsg = 'Triangle_id %d does not exist in mesh: %s' %(triangle_id,
+                                                                    str(self))
+                assert 0 <= triangle_id < len(self), errmsg
+                k = self.k = triangle_id
+            
 
             x, y = self.get_centroid_coordinates()[k]
@@ -725 +737 @@
             msg += '  Triangle #%d with centroid (%.4f, %.4f), ' %(k, x, y)
             msg += 'area = %.4f and radius = %.4f ' %(area, radius)
-            msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
-            msg += 'during last time interval. Quantities below '
-            msg += 'are reported at their present value, and not what '
-            msg += 'they were at the time the maximal speed was attained.'
-            msg += 'To see this, rerun the model with yieldsteps smaller '
-            msg += 'than the smallest internal timestep reported.'
+            if triangle_id is None:
+                msg += 'had the largest computed speed: %.6f m/s ' %(max_speed)
+            else:
+                msg += 'had computed speed: %.6f m/s ' %(max_speed)
+                
             if max_speed > 0.0:
                 msg += '(timestep=%.6f)\n' %(radius/max_speed)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -262 +262 @@
     def get_node(self, i,
                  absolute=False):
-        """Return node coordinates.
-
-        The nodes are ordered in an Nx2 array where N is the number of nodes.
-        This is the same format they were provided in the constructor
-        i.e. without any duplication.
+        """Return node coordinates for triangle i.
 
         Boolean keyword argument absolute determines whether coordinates
@@ -597 +593 @@
 
         
+        

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -860 +860 @@
         return str
 
+    def get_triangle_containing_point(self, point):
+        """Return triangle id for triangle containing specifiend point (x,y)
+
+        If point isn't within mesh, raise exception
+
+        """
+        
+        # FIXME(Ole): This function is currently brute force
+        # because I needed it for diagnostics.
+        # We should make it fast - probably based on the
+        # quad tree structure.
+        from anuga.utilities.polygon import is_outside_polygon,\
+             is_inside_polygon
+
+        polygon = self.get_boundary_polygon()
+        #print polygon, point
+        
+        if is_outside_polygon(point, polygon):
+            msg = 'Point %s is outside mesh' %str(point)
+            raise Exception, msg
+
+
+        V = self.get_vertex_coordinates(absolute=True)
+
+        # FIXME: Horrible brute force
+        for i, triangle in enumerate(self.triangles):
+            poly = V[3*i:3*i+3]
+            #print i, poly
+
+            if is_inside_polygon(point, poly, closed=True):
+                return i
+                
+        return

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1195 +1195 @@
 
         assert allclose(absolute_points, boundary_polygon)
-        
+
+    def test_get_triangle_containing_point(self):
+
+        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]]
+        mesh = Mesh(points, vertices)
+        
+        mesh.check_integrity()
+
+
+        try:
+            id = mesh.get_triangle_containing_point([3.0, 5.0])
+        except:
+            pass
+        else:
+            msg = 'Should have caught point outside polygon (Non)'            
+            raise Exception, msg
+            
+        id = mesh.get_triangle_containing_point([0.5, 1.0])
+        assert id == 0
+
+        id = mesh.get_triangle_containing_point([1.0, 3.0])
+        assert id == 3        
+
+        for i, point in enumerate(mesh.get_centroid_coordinates()):
+            id = mesh.get_triangle_containing_point(point)
+            assert id == i        
+            
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon_IIIa')
+    #suite = unittest.makeSuite(Test_Mesh,'test_get_triangle_containing_point')
     suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -462 +462 @@
 
         
-    def timestepping_statistics(self, track_speeds=False):
+    def timestepping_statistics(self,
+                                track_speeds=False,
+                                triangle_id=None):        
         """Return string with time stepping statistics for printing or logging
 
@@ -475 +477 @@
 
         # Call basic machinery from parent class
-        msg = Generic_Domain.timestepping_statistics(self, track_speeds)
+        msg = Generic_Domain.timestepping_statistics(self,
+                                                     track_speeds,
+                                                     triangle_id)
 
         if track_speeds is True:
@@ -482 +486 @@
             qwidth = self.qwidth
         
-            # Triangle with maximum speed
+            # Selected triangle
             k = self.k
 

anuga_validation/okushiri_2005/run_okushiri.py

@@ -82 +82 @@
 
 
+# Select triangle containing ch5 for diagnostic output
+# around known gauge
+triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
+# This should get triangle id 32833 with centroid (4.5244, 1.1972)
+
+
 #-------------------------
 # Evolve through time
@@ -89 +95 @@
 
 for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
-    domain.write_time(track_speeds=False)
+    print domain.timestepping_statistics(track_speeds=True,
+                                         triangle_id=triangle_id)
 
 print 'That took %.2f seconds' %(time.time()-t0)