Changeset 705 for inundation

Timestamp:

Dec 14, 2004, 12:59:59 PM (20 years ago)

Author:

duncan

Message:

fixed bug in least_squares algorithm.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (3 diffs)
• quad.py (modified) (1 diff)
• test_least_squares.py (modified) (3 diffs)
• test_quad.py (modified) (1 diff)
• wiki/issues.txt (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -245 +245 @@
             #Find vertices near x
             candidate_vertices = root.search(x[0], x[1])
-
             is_more_elements = True
-            if candidate_vertices == []:
-                # The point isn't even within the root cell!
-                is_more_elements = False
-                element_found = False
-            else:
-                element_found, sigma0, sigma1, sigma2, k = \
-                    self.search_triangles_of_vertices(candidate_vertices, x)
-
+            
+            element_found, sigma0, sigma1, sigma2, k = \
+                self.search_triangles_of_vertices(candidate_vertices, x)
             while not element_found and is_more_elements: 
-                candidate_vertices = root.expand_search()
-                if candidate_vertices == []:
-                    # All the triangles have been searched.
+                candidate_vertices, branch = root.expand_search()
+                if branch == []:
+                    # Searching all the verts from the root cell that haven't
+                    # been searched.  This is the last try
+                    element_found, sigma0, sigma1, sigma2, k = \
+                      self.search_triangles_of_vertices(candidate_vertices, x)
                     is_more_elements = False
                 else:
                     element_found, sigma0, sigma1, sigma2, k = \
                       self.search_triangles_of_vertices(candidate_vertices, x)
-                    
-                
+                     
             
             #Update interpolation matrix A if necessary     
@@ -301 +297 @@
             sigma0 = -10.0
             sigma1 = -10.0
+            k = -10.0
 
             #For all vertices in same cell as point x
@@ -317 +314 @@
                     #print "PDSG - xi1", xi1
                     #print "PDSG - xi2", xi2
-                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))" \
-                    #      % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
+                   #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))" \
+                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
                     
                     #Get the three normals

inundation/ga/storm_surge/pyvolution/quad.py

@@ -122 +122 @@
                 #print "cell ", cell.show() 
                 points += cell.retrieve()
-        return points
+        return points, self.branch
 
 

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -83 +83 @@
                       0., 0. , 0., 0., 0., 0.]]
         assert allclose(interp.get_A(), answer)
-        
+     
+    def test_quad_treeII(self):
+        p0 = [-66.0, 14.0]
+        p1 = [14.0, -66.0]
+        p2 = [14.0, 14.0]
+        p3 = [60.0, 20.0]
+        p4 = [10.0, 60.0]
+        p5 = [60.0, 60.0]
+        
+        points = [p0, p1, p2, p3, p4, p5]
+        triangles = [ [0, 1, 2],
+                      [3, 2, 1],
+                      [0, 2, 4],
+                      [0, 2, 4],
+                      [4, 2, 5],
+                      [5, 2, 3]]   
+
+        data = [ [-26.0,-26.0] ] 
+        
+        interp = Interpolation(points, triangles, data, alpha = 0.0)
+        #print "PDSG - interp.get_A()", interp.get_A()
+        answer =  [ [ 0.5,  0.5,  0.0,  0.,
+                      0., 0.]]
+        assert allclose(interp.get_A(), answer)
+        interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
+        #print "PDSG -30,-30 - interp.get_A()", interp.get_A()
+        answer =  [ [ 0.0,  0.0,  0.0,  0.,
+                      0., 0.]]
+        assert allclose(interp.get_A(), answer)
+
+        
+        #point outside of quad tree root cell
+        interp.set_point_coordinates([[-70, -70]])
+        #print "PDSG -70,-70 interp.get_A()", interp.get_A()
+        answer =  [ [ 0.0,  0.0,  0.0,  0.,
+                      0., 0. ]]
+        assert allclose(interp.get_A(), answer)
+           
 
     def test_datapoints_at_vertices(self):
@@ -158 +195 @@
         
         interp = Interpolation(points, vertices, data)
-
+        #print "interp.get_A()", interp.get_A()
         assert allclose(sum(interp.get_A(), axis=1), 1.0)
         
@@ -793 +830 @@
 
     suite = unittest.makeSuite(TestCase,'test')
+    #suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

inundation/ga/storm_surge/pyvolution/test_quad.py

@@ -127 +127 @@
         assert type(result) in [types.ListType,types.TupleType],\
                                 'should be a list'
-        
+        #print "result",result
         self.assertEqual(result, [0,1,2,3])
 

inundation/ga/storm_surge/pyvolution/wiki/issues.txt

@@ -2 +2 @@
 OPEN ISSUES:
 ------------
+Issue (least_squares): The current is-a-point-in-a-triangle algorithm
+is slow if a point is outside of the mesh.  It will check each
+triangle 3 times to see if the point is in that triangle.
+Importance: low
+Suggested Action: keep a list of triangles searched./ form a polygon
+of the mesh boundary and throw out points not in the polygon. 
+
 Issue (interpolate_sww): Give a warning if points are outside the mesh.
 Importance: Mid