Changeset 4515


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r4475 r4515  
    444444
    445445        Using the mesh boundary, derive a bounding polygon for this mesh.
    446         If multiple vertex values are present, the algorithm will select the
    447         path that contains the mesh.
     446        If multiple vertex values are present (vertices stored uniquely),
     447        the algorithm will select the path that contains the entire mesh.
    448448
    449449        All points are in absolute UTM coordinates
     
    451451       
    452452        from Numeric import allclose, sqrt, array, minimum, maximum
    453         from anuga.utilities.numerical_tools import angle, ensure_numeric        
     453        from anuga.utilities.numerical_tools import angle, ensure_numeric     
    454454
    455455
     
    464464        inverse_segments = {}
    465465        p0 = None
    466         mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
     466        mindist = sqrt(sum((pmax-pmin)**2)) # Start value across entire mesh
    467467        for i, edge_id in self.boundary.keys():
    468468            # Find vertex ids for boundary segment
     
    481481            dist_B = sqrt(sum((B-pmin)**2))
    482482
    483             #Find lower leftmost point
     483            # Find lower leftmost point
    484484            if dist_A < mindist:
    485485                mindist = dist_A
     
    497497            # Register potential paths from A to B
    498498            if not segments.has_key(tuple(A)):
    499                 segments[tuple(A)] = [] # Empty list for candidate points               
    500                
     499                segments[tuple(A)] = [] # Empty list for candidate points
     500                         
    501501            segments[tuple(A)].append(B)               
    502502
    503503
    504         #Start with smallest point and follow boundary (counter clock wise)
     504        # Start with smallest point and follow boundary (counter clock wise)
    505505        polygon = [p0]      # Storage for final boundary polygon
    506         point_registry = {} # Keep track of storage to avoid multiple runs around
    507                             # boundary. This will only be the case if there are
    508                             # more than one candidate.
     506        point_registry = {} # Keep track of storage to avoid multiple runs
     507                            # around boundary. This will only be the case if
     508                            # there are more than one candidate.
    509509                            # FIXME (Ole): Perhaps we can do away with polygon
    510510                            # and use only point_registry to save space.
     
    512512        point_registry[tuple(p0)] = 0                           
    513513                           
    514         #while len(polygon) < len(self.boundary):
    515514        while len(point_registry) < len(self.boundary):
    516515
    517516            candidate_list = segments[tuple(p0)]
    518517            if len(candidate_list) > 1:
    519                 # Multiple points detected (this will be the case for meshes with
    520                 # duplicate points as those used for discontinuous triangles).
    521                 # Take the candidate that is furthest to the clockwise direction,
    522                 # as that will follow the boundary.
    523 
     518                # Multiple points detected (this will be the case for meshes
     519                # with duplicate points as those used for discontinuous
     520                # triangles with vertices stored uniquely).
     521                # Take the candidate that is furthest to the clockwise
     522                # direction, as that will follow the boundary.
     523                #
     524                # This will also be the case for pathological triangles
     525                # that have no neighbours.
    524526
    525527                if verbose:
     
    527529                          %(str(p0), candidate_list)
    528530
    529 
    530531                # Check that previous are not in candidate list
    531532                #for p in candidate_list:
    532533                #    assert not allclose(p0, p)
    533534
    534 
    535535                # Choose vector against which all angles will be measured
    536536                if len(polygon) > 1:   
    537                     v_prev = p0 - polygon[-2] # Vector that leads to p0 from previous point
     537                    v_prev = p0 - polygon[-2] # Vector that leads to p0
     538                                              # from previous point
    538539                else:
    539                     # FIXME (Ole): What do we do if the first point has multiple
    540                     # candidates?
     540                    # FIXME (Ole): What do we do if the first point has
     541                    # multiple candidates?
    541542                    # Being the lower left corner, perhaps we can use the
    542                     # vector [1, 0], but I really don't know if this is completely
    543                     # watertight.
     543                    # vector [1, 0], but I really don't know if this is
     544                    # completely watertight.
    544545                    v_prev = [1.0, 0.0]
    545546                   
     
    549550                for pc in candidate_list:
    550551
    551                    
    552552                    vc = pc-p0  # Candidate vector (from p0 to candidate pt)
    553553                   
     
    556556                    ac = angle(vc, v_prev)
    557557                    if ac > pi:
    558                         # Give preference to angles on the right hand side of v_prev
    559                         #print 'pc = %s, changing angle from %f to %f' %(pc, ac*180/pi, (ac-2*pi)*180/pi)
     558                        # Give preference to angles on the right hand side
     559                        # of v_prev
     560                        # print 'pc = %s, changing angle from %f to %f'\
     561                        # %(pc, ac*180/pi, (ac-2*pi)*180/pi)
    560562                        ac = ac-2*pi
    561563
    562                     # take the minimal angle corresponding to the rightmost vector
     564                    # Take the minimal angle corresponding to the
     565                    # rightmost vector
    563566                    if ac < minimum_angle:
    564567                        minimum_angle = ac
     
    567570
    568571                if verbose is True:
    569                     print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
     572                    print '  Best candidate %s, angle %f'\
     573                          %(p1, minimum_angle*180/pi)
    570574                   
    571575            else:
    572576                p1 = candidate_list[0]
    573577
     578               
    574579            if point_registry.has_key(tuple(p1)):
    575                 # We have completed the boundary polygon - yeehaa
    576                 break
     580                # We have reached a point already visited.
     581               
     582                if allclose(p1, polygon[0]):
     583                    # If it is the initial point, the polygon is complete.
     584                   
     585                    if verbose is True:
     586                        print '  Stop criterion fulfilled at point %s' %p1
     587                        print polygon               
     588                       
     589                    # We have completed the boundary polygon - yeehaa
     590                    break
     591                else:   
     592                    # The point already visited is not the initial point
     593                    # This would be a pathological triangle, but the
     594                    # algorithm must be able to deal with this
     595                    pass
     596   
    577597            else:
     598                # We are still finding new points on the boundary
    578599                point_registry[tuple(p1)] = len(point_registry)
    579600           
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r4171 r4515  
    809809            assert is_inside_polygon(p, P)
    810810
     811           
     812    def test_boundary_polygon_IIIa(self):
     813        """test_boundary_polygon_IIIa - Check pathological situation where
     814        one triangle has no neighbours. This may be the case if a mesh
     815        is partitioned using pymetis.
     816        """
     817
     818        from Numeric import zeros, Float
     819
     820
     821        #Points
     822        a = [0.0, 0.0] #0
     823        b = [0.0, 0.5] #1
     824        c = [0.0, 1.0] #2
     825        d = [0.5, 0.0] #3
     826        e = [0.5, 0.5] #4
     827        f = [1.0, 0.0] #5
     828        g = [1.0, 0.5] #6
     829        h = [1.0, 1.0] #7
     830       
     831        # Add pathological triangle with no neighbours to an otherwise
     832        # trivial mesh
     833       
     834        points = [a, b, c, d, e, f, g, h]
     835
     836        #cbe, aeb, dea, fed, ghe (pathological triangle)
     837        vertices = [[2,1,4], [0,4,1], [3,4,0], [5,4,3],
     838                    [6,7,4]]   
     839                   
     840        mesh = Mesh(points, vertices)
     841        mesh.check_integrity()
     842
     843        P = mesh.get_boundary_polygon(verbose=False)
     844
     845       
     846        assert len(P) == 9
     847       
     848        # Note that point e appears twice!
     849        assert allclose(P, [a, d, f, e, g, h, e, c, b])
     850
     851        for p in points:
     852            msg = 'Point %s is not inside polygon %s'\
     853            %(p, P)     
     854            assert is_inside_polygon(p, P), msg             
     855
     856
     857                                       
     858           
     859           
    811860
    812861    def test_boundary_polygon_IV(self):
     
    11491198#-------------------------------------------------------------
    11501199if __name__ == "__main__":
    1151     #suite = unittest.makeSuite(Test_Mesh,'test_mesh_get_boundary_polygon_with_georeferencing')
     1200    #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon_IIIa')
    11521201    suite = unittest.makeSuite(Test_Mesh,'test')
    11531202    runner = unittest.TextTestRunner()
Note: See TracChangeset for help on using the changeset viewer.