Changeset 4515

Timestamp:

May 30, 2007, 5:25:51 AM (17 years ago)

Author:

ole

Message:

Modified get_bounding_polygon to deal cleanly with pathological (neighbourless) triangles. This should address ticket:160

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• neighbour_mesh.py (modified) (10 diffs)
• test_neighbour_mesh.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -444 +444 @@
 
         Using the mesh boundary, derive a bounding polygon for this mesh.
-        If multiple vertex values are present, the algorithm will select the
-        path that contains the mesh.
+        If multiple vertex values are present (vertices stored uniquely), 
+        the algorithm will select the path that contains the entire mesh.
 
         All points are in absolute UTM coordinates
@@ -451 +451 @@
         
         from Numeric import allclose, sqrt, array, minimum, maximum
-        from anuga.utilities.numerical_tools import angle, ensure_numeric        
+        from anuga.utilities.numerical_tools import angle, ensure_numeric     
 
 
@@ -464 +464 @@
         inverse_segments = {}
         p0 = None
-        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
+        mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh
         for i, edge_id in self.boundary.keys():
             # Find vertex ids for boundary segment
@@ -481 +481 @@
             dist_B = sqrt(sum((B-pmin)**2))
 
-            #Find lower leftmost point
+            # Find lower leftmost point
             if dist_A < mindist:
                 mindist = dist_A
@@ -497 +497 @@
             # Register potential paths from A to B
             if not segments.has_key(tuple(A)):
-                segments[tuple(A)] = [] # Empty list for candidate points                
-                
+                segments[tuple(A)] = [] # Empty list for candidate points
+                          
             segments[tuple(A)].append(B)                
 
 
-        #Start with smallest point and follow boundary (counter clock wise)
+        # Start with smallest point and follow boundary (counter clock wise)
         polygon = [p0]      # Storage for final boundary polygon
-        point_registry = {} # Keep track of storage to avoid multiple runs around
-                            # boundary. This will only be the case if there are
-                            # more than one candidate.
+        point_registry = {} # Keep track of storage to avoid multiple runs 
+                            # around boundary. This will only be the case if 
+                            # there are more than one candidate.
                             # FIXME (Ole): Perhaps we can do away with polygon
                             # and use only point_registry to save space.
@@ -512 +512 @@
         point_registry[tuple(p0)] = 0                            
                             
-        #while len(polygon) < len(self.boundary):
         while len(point_registry) < len(self.boundary):
 
             candidate_list = segments[tuple(p0)]
             if len(candidate_list) > 1:
-                # Multiple points detected (this will be the case for meshes with
-                # duplicate points as those used for discontinuous triangles).
-                # Take the candidate that is furthest to the clockwise direction,
-                # as that will follow the boundary.
-
+                # Multiple points detected (this will be the case for meshes 
+                # with duplicate points as those used for discontinuous 
+                # triangles with vertices stored uniquely).
+                # Take the candidate that is furthest to the clockwise 
+                # direction, as that will follow the boundary.
+                #
+                # This will also be the case for pathological triangles
+                # that have no neighbours.
 
                 if verbose:
@@ -527 +529 @@
                           %(str(p0), candidate_list)
 
-
                 # 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 from previous point
+                    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
-                    # candidates?
+                    # FIXME (Ole): What do we do if the first point has 
+                    # multiple candidates?
                     # Being the lower left corner, perhaps we can use the
-                    # vector [1, 0], but I really don't know if this is completely
-                    # watertight.
+                    # vector [1, 0], but I really don't know if this is 
+                    # completely watertight.
                     v_prev = [1.0, 0.0]
                     
@@ -549 +550 @@
                 for pc in candidate_list:
 
-                    
                     vc = pc-p0  # Candidate vector (from p0 to candidate pt)
                    
@@ -556 +556 @@
                     ac = angle(vc, v_prev)
                     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)
+                        # 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
+                    # Take the minimal angle corresponding to the 
+                    # rightmost vector
                     if ac < minimum_angle:
                         minimum_angle = ac
@@ -567 +570 @@
 
                 if verbose is True:
-                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
+                    print '  Best candidate %s, angle %f'\
+                          %(p1, minimum_angle*180/pi)
                     
             else:
                 p1 = candidate_list[0]
 
+                
             if point_registry.has_key(tuple(p1)):
-                # We have completed the boundary polygon - yeehaa
-                break 
+                # We have reached a point already visited. 
+                
+                if allclose(p1, polygon[0]):
+                    # If it is the initial point, the polygon is complete. 
+                    
+                    if verbose is True:
+                        print '  Stop criterion fulfilled at point %s' %p1
+                        print polygon               
+                        
+                    # We have completed the boundary polygon - yeehaa
+                    break
+                else:    
+                    # The point already visited is not the initial point
+                    # This would be a pathological triangle, but the 
+                    # algorithm must be able to deal with this
+                    pass
+   
             else:
+                # We are still finding new points on the boundary
                 point_registry[tuple(p1)] = len(point_registry)
             

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -809 +809 @@
             assert is_inside_polygon(p, P)
 
+            
+    def test_boundary_polygon_IIIa(self):
+        """test_boundary_polygon_IIIa - Check pathological situation where
+        one triangle has no neighbours. This may be the case if a mesh
+        is partitioned using pymetis.
+        """
+
+        from Numeric import zeros, Float
+
+
+        #Points
+        a = [0.0, 0.0] #0
+        b = [0.0, 0.5] #1
+        c = [0.0, 1.0] #2
+        d = [0.5, 0.0] #3
+        e = [0.5, 0.5] #4
+        f = [1.0, 0.0] #5
+        g = [1.0, 0.5] #6
+        h = [1.0, 1.0] #7
+       
+        # Add pathological triangle with no neighbours to an otherwise
+        # trivial mesh
+        
+        points = [a, b, c, d, e, f, g, h]
+
+        #cbe, aeb, dea, fed, ghe (pathological triangle)
+        vertices = [[2,1,4], [0,4,1], [3,4,0], [5,4,3],
+                    [6,7,4]]    
+                    
+        mesh = Mesh(points, vertices)
+        mesh.check_integrity()
+
+        P = mesh.get_boundary_polygon(verbose=False)
+
+        
+        assert len(P) == 9
+        
+        # Note that point e appears twice!
+        assert allclose(P, [a, d, f, e, g, h, e, c, b])
+
+        for p in points:
+            msg = 'Point %s is not inside polygon %s'\
+            %(p, P)     
+            assert is_inside_polygon(p, P), msg             
+
+
+                                        
+            
+            
 
     def test_boundary_polygon_IV(self):
@@ -1149 +1198 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
+    #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon_IIIa')
     suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()