Changeset 3688

Timestamp:

Oct 4, 2006, 9:32:03 AM (16 years ago)

Author:

ole

Message:

Fixed get_boundary polygon based on information from new unit test.
It now works for all examples of meshes provided.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• neighbour_mesh.py (modified) (3 diffs)
• test_neighbour_mesh.py (modified) (5 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -517 +517 @@
 
 
+                # Check that previous are not in candidate list
+                #for p in candidate_list:
+                #    assert not allclose(p0, p)
+
+
                 # Choose vector against which all angles will be measured
                 if len(polygon) > 1:    
-                    v_prev = p0 - polygon[-2] # Vector that leads to p0
-
-                    # Normalise
-                    #v_prev /= sqrt(sum(v_prev**2))  
+                    v_prev = p0 - polygon[-2] # Vector that leads to p0 from previous point
                 else:
                     # FIXME (Ole): What do we do if the first point has multiple
@@ -531 +533 @@
                     v_prev = [1.0, 0.0]
                     
-                    
+
                 # Choose candidate with minimum angle    
                 minimum_angle = 2*pi
                 for pc in candidate_list:
-                    vc = pc-p0  # Candidate vector
-                    # Normalise
-                    #vc /= sqrt(sum(vc**2))                      
+
                     
+                    vc = pc-p0  # Candidate vector (from p0 to candidate pt)
+                   
                     # Angle between each candidate and the previous vector
                     # in [-pi, pi]
-
                     ac = angle(vc, v_prev)
-                    if ac > pi: ac = pi-ac
+                    if ac > pi:
+                        # Give preference to angles on the right hand side of v_prev
+                        #print 'pc = %s, changing angle from %f to %f' %(pc, ac*180/pi, (ac-2*pi)*180/pi)
+                        ac = ac-2*pi
 
                     # take the minimal angle corresponding to the rightmost vector
@@ -553 +557 @@
                 if verbose is True:
                     print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
-                
+                    
             else:
                 p1 = candidate_list[0]

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -910 +910 @@
 
 
-    def xtest_boundary_polygon_VI(self):
+    def test_boundary_polygon_VI(self):
         """test_boundary_polygon_VI(self)
 
@@ -942 +942 @@
                   [  75735.4765625 ,  23762.00585938],
                   [  52341.70703125,  38563.39453125]]
+
+        ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation        
 
         triangles = [[19, 0,15],
@@ -969 +971 @@
         Pref = mesh.get_boundary_polygon()
 
-        plot_polygons([ensure_numeric(Pref)], 'goodP')
+        #plot_polygons([ensure_numeric(Pref)], 'goodP')
 
         for p in points:
@@ -1044 +1046 @@
                   [  35406.3359375 ,  79332.9140625 ]]
 
-        #points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
+        scaled_points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
 
         triangles = [[ 0, 1, 2],
@@ -1068 +1070 @@
                      [60,61,62]]
 
+
+        # First use scaled points for ease of debugging
+        mesh = Mesh(scaled_points, triangles)
+        mesh.check_integrity()
+        P = mesh.get_boundary_polygon()
+
+        for p in scaled_points:
+            assert is_inside_polygon(p, P)            
+
+        # Then use original points and test        
         mesh = Mesh(points, triangles)
         mesh.check_integrity()
         P = mesh.get_boundary_polygon()
 
-        plot_polygons([ensure_numeric(P)], 'badP')
-        print P
-
         for p in points:
-            #assert is_inside_polygon(p, P)            
-            if not is_inside_polygon(p, P):
-                print 'Point %s is not in P' %(p)
-            else:
-                print 'Point %s OK' %(p)
-
-        #for i in range(len(Pref)):
-        #    print P[i], Pref[i]
-        #
-        #print ensure_numeric(P) - ensure_numeric(Pref)
-
-
+            assert is_inside_polygon(p, P)            
+
+        assert allclose(P, Pref)    
 
     def test_lone_vertices(self):