Changeset 2709


Ignore:
Timestamp:
Apr 12, 2006, 2:28:05 PM (18 years ago)
Author:
ole
Message:

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

Location:
inundation
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/domain.py

    r2648 r2709  
    11"""Class Domain - 2D triangular domains for finite-volume computations of
    2    the shallow water wave equation
     2   conservation laws.
    33
    44
  • inundation/pyvolution/mesh.py

    r2683 r2709  
    66
    77from general_mesh import General_mesh
     8from math import pi, sqrt
     9       
    810
    911class Mesh(General_mesh):
     
    125127                #of inscribed circle is computed
    126128
    127                 from math import sqrt
    128129                a = sqrt((x0-x1)**2+(y0-y1)**2)
    129130                b = sqrt((x1-x2)**2+(y1-y2)**2)
     
    165166    def set_to_inscribed_circle(self,safety_factor = 1):
    166167        #FIXME phase out eventually
    167         from math import sqrt
    168168        N = self.number_of_elements
    169169        V = self.vertex_coordinates
     
    414414
    415415
    416     def get_boundary_polygon(self):
    417         """Return bounding polygon as a list of points
    418 
    419         FIXME: If triangles are listed as discontinuous
    420         (e.g vertex coordinates listed multiple times),
    421         this may not work as expected.
    422         """
     416    def get_boundary_polygon(self, verbose = False):
     417        """Return bounding polygon for mesh (counter clockwise)
     418
     419        Using the mesh boundary, derive a bounding polygon for this mesh.
     420        If multiple vertex values are present, the algorithm will select the
     421        path that contains the mesh.
     422        """
     423       
    423424        from Numeric import allclose, sqrt, array, minimum, maximum
    424 
    425 
    426 
    427         #V = self.get_vertex_coordinates()
     425        from utilities.numerical_tools import angle, ensure_numeric       
     426
     427
     428        # Get mesh extent
     429        xmin, xmax, ymin, ymax = self.get_extent()
     430        pmin = ensure_numeric([xmin, ymin])
     431        pmax = ensure_numeric([xmax, ymax])       
     432
     433
     434        # Assemble dictionary of boundary segments and choose starting point
    428435        segments = {}
    429 
    430         #pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
    431         #pmax = (max(self.coordinates[:,0]), max(self.coordinates[:,1]))
    432 
    433         #FIXME:Can this be written more compactly, e.g.
    434         #using minimum and maximium?
    435         pmin = array( [min(self.coordinates[:,0]),
    436                        min(self.coordinates[:,1]) ] )
    437 
    438         pmax = array( [max(self.coordinates[:,0]),
    439                        max(self.coordinates[:,1]) ] )
    440 
    441         mindist = sqrt(sum( (pmax-pmin)**2 ))
     436        inverse_segments = {}
     437        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
    442438        for i, edge_id in self.boundary.keys():
    443             #Find vertex ids for boundary segment
     439            # Find vertex ids for boundary segment
    444440            if edge_id == 0: a = 1; b = 2
    445441            if edge_id == 1: a = 2; b = 0
    446442            if edge_id == 2: a = 0; b = 1
    447443
    448             A = tuple(self.get_vertex_coordinate(i, a))
    449             B = tuple(self.get_vertex_coordinate(i, b))
    450 
    451             #Take the point closest to pmin as starting point
    452             #Note: Could be arbitrary, but nice to have
    453             #a unique way of selecting
    454             dist_A = sqrt(sum( (A-pmin)**2 ))
    455             dist_B = sqrt(sum( (B-pmin)**2 ))
    456 
    457             #Find minimal point
     444            A = self.get_vertex_coordinate(i, a) # Start
     445            B = self.get_vertex_coordinate(i, b) # End
     446
     447            # Take the point closest to pmin as starting point
     448            # Note: Could be arbitrary, but nice to have
     449            # a unique way of selecting
     450            dist_A = sqrt(sum((A-pmin)**2))
     451            dist_B = sqrt(sum((B-pmin)**2))
     452
     453            #Find lower leftmost point
    458454            if dist_A < mindist:
    459455                mindist = dist_A
     
    464460
    465461
     462            # Sanity check
    466463            if p0 is None:
    467                 raise 'Weird'
    468                 p0 = A #We need a starting point (FIXME)
    469 
    470                 print 'A', A
    471                 print 'B', B
    472                 print 'pmin', pmin
    473                 print
    474 
    475             segments[A] = B
     464                raise Exception('Impossible')
     465
     466
     467            # Register potential paths from A to B
     468            if not segments.has_key(tuple(A)):
     469                segments[tuple(A)] = [] # Empty list for candidate points               
     470               
     471            segments[tuple(A)].append(B)               
     472
     473
    476474
    477475
    478476        #Start with smallest point and follow boundary (counter clock wise)
    479         polygon = [p0]
    480         while len(polygon) < len(self.boundary):
    481             p1 = segments[p0]
     477        polygon = [p0]      # Storage for final boundary polygon
     478        point_registry = {} # Keep track of storage to avoid multiple runs around boundary
     479                            # This will only be the case if there are more than one candidate
     480                            # FIXME (Ole): Perhaps we can do away with polygon and use
     481                            # only point_registry to save space.
     482
     483        point_registry[tuple(p0)] = 0                           
     484                           
     485        #while len(polygon) < len(self.boundary):
     486        while len(point_registry) < len(self.boundary):
     487
     488            candidate_list = segments[tuple(p0)]
     489            if len(candidate_list) > 1:
     490                # Multiple points detected
     491                # Take the candidate that is furthest to the clockwise direction,
     492                # as that will follow the boundary.
     493
     494
     495                if verbose:
     496                    print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list)
     497
     498
     499                # Choose vector against which all angles will be measured
     500                if len(polygon) > 1:   
     501                    v_prev = p0 - polygon[-2] # Vector that leads to p0
     502                else:
     503                    # FIXME (Ole): What do we do if the first point has multiple candidates?
     504                    # Being the lower left corner, perhaps we can use the
     505                    # vector [1, 0], but I really don't know if this is completely
     506                    # watertight.
     507                    # Another option might be v_prev = [1.0, 0.0]
     508                    v_prev = [1.0, 0.0]
     509                   
     510
     511                   
     512                # Choose candidate with minimum angle   
     513                minimum_angle = 2*pi
     514                for pc in candidate_list:
     515                    vc = pc-p0  # Candidate vector
     516                   
     517                    # Angle between each candidate and the previous vector in [-pi, pi]
     518                    ac = angle(vc, v_prev)
     519                    if ac > pi: ac = pi-ac
     520
     521                    # take the minimal angle corresponding to the rightmost vector
     522                    if ac < minimum_angle:
     523                        minimum_angle = ac
     524                        p1 = pc             # Best candidate
     525                       
     526
     527                if verbose is True:
     528                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
     529               
     530            else:
     531                p1 = candidate_list[0]
     532
     533            if point_registry.has_key(tuple(p1)):
     534                # We have completed the boundary polygon - yeehaa
     535                break
     536            else:
     537                point_registry[tuple(p1)] = len(point_registry)
     538           
    482539            polygon.append(p1)
    483540            p0 = p1
     541
    484542
    485543        return polygon
     
    494552
    495553        from config import epsilon
    496         from math import pi
    497554        from utilities.numerical_tools import anglediff
    498555
     
    537594
    538595                #Normalise
    539                 from math import sqrt
    540596                l_u = sqrt(u[0]*u[0] + u[1]*u[1])
    541597                l_v = sqrt(v[0]*v[0] + v[1]*v[1])               
  • inundation/pyvolution/test_mesh.py

    r2532 r2709  
    849849
    850850
     851    def test_boundary_polygon_V(self):
     852        """Create a discontinuous mesh (duplicate vertices)
     853        and check that boundary is as expected
     854       
     855        """
     856        from mesh import Mesh
     857        from Numeric import zeros, Float
     858       
     859
     860        #Points
     861        a = [0.0, 0.0] #0
     862        b = [0.0, 0.5] #1
     863        c = [0.0, 1.0] #2
     864        d = [0.5, 0.0] #3
     865        e = [0.5, 0.5] #4
     866        f = [1.0, 0.0] #5
     867        g = [1.0, 0.5] #6
     868        h = [1.0, 1.0] #7
     869        i = [1.5, 0.5] #8
     870
     871        #Duplicate points for triangles edg [4,3,6] (central) and
     872        #gid [6,8,7] (top right boundary) to them disconnected
     873        #from the others
     874
     875        e0 = [0.5, 0.5] #9
     876        d0 = [0.5, 0.0] #10
     877        g0 = [1.0, 0.5] #11
     878        i0 = [1.5, 0.5] #12       
     879       
     880
     881        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
     882
     883
     884
     885        #dea, bae, bec, fgd,
     886        #edg, ghe, gfi, gih
     887        #vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
     888        #             [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
     889
     890
     891        #dea, bae, bec, fgd,
     892        #e0d0g0, ghe, gfi, g0i0h
     893        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
     894                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]       
     895
     896        mesh = Mesh(points, vertices)
     897
     898        mesh.check_integrity()
     899
     900        P = mesh.get_boundary_polygon()
     901
     902        #print P
     903       
     904        assert len(P) == 8
     905        assert allclose(P, [a, d, f, i, h, e, c, b])
     906        assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
     907       
     908
     909        for p in points:
     910            #print p, P
     911            assert inside_polygon(p, P)
    851912
    852913
    853914#-------------------------------------------------------------
    854915if __name__ == "__main__":
    855     suite = unittest.makeSuite(Test_Mesh,'test')
     916    suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
     917    #suite = unittest.makeSuite(Test_Mesh,'test')
    856918    runner = unittest.TextTestRunner()
    857919    runner.run(suite)
  • inundation/pyvolution/test_shallow_water.py

    r2620 r2709  
    30403040
    30413041
    3042         #Create an triangle shaped domain (reusing coordinates from domain 1),
     3042        #Create a triangle shaped domain (reusing coordinates from domain 1),
    30433043        #formed from the lower and right hand  boundaries and
    30443044        #the sw-ne diagonal
     
    33763376if __name__ == "__main__":
    33773377    suite = unittest.makeSuite(Test_Shallow_Water,'test')
    3378     runner = unittest.TextTestRunner()
     3378    runner = unittest.TextTestRunner(verbosity=1)   
    33793379    runner.run(suite)
  • inundation/utilities/numerical_tools.py

    r2704 r2709  
    5757   
    5858def anglediff(v0, v1):
    59     """Compute difference between angle of vector x0, y0 and x1, y1.
     59    """Compute difference between angle of vector v0 (x0, y0) and v1 (x1, y1).
    6060    This is used for determining the ordering of vertices,
    6161    e.g. for checking if they are counter clockwise.
  • inundation/utilities/test_numerical_tools.py

    r2704 r2709  
    4747       
    4848        assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
    49         assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)   
     49        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
     50
     51        #From test_get_boundary_polygon_V
     52        v_prev = [-0.5, -0.5]
     53        vc = [ 0.0,  -0.5]
     54        assert allclose(angle(vc, v_prev)/pi*180, 45.0)
     55
     56        vc = [ 0.5,  0.0]
     57        assert allclose(angle(vc, v_prev)/pi*180, 135.0)
     58
     59        vc = [ -0.5,  0.5]
     60        assert allclose(angle(vc, v_prev)/pi*180, 270.0)               
     61
     62
    5063               
    5164       
Note: See TracChangeset for help on using the changeset viewer.