Changeset 3688


Ignore:
Timestamp:
Oct 4, 2006, 9:32:03 AM (17 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:
2 edited

Legend:

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

    r3684 r3688  
    517517
    518518
     519                # Check that previous are not in candidate list
     520                #for p in candidate_list:
     521                #    assert not allclose(p0, p)
     522
     523
    519524                # Choose vector against which all angles will be measured
    520525                if len(polygon) > 1:   
    521                     v_prev = p0 - polygon[-2] # Vector that leads to p0
    522 
    523                     # Normalise
    524                     #v_prev /= sqrt(sum(v_prev**2)) 
     526                    v_prev = p0 - polygon[-2] # Vector that leads to p0 from previous point
    525527                else:
    526528                    # FIXME (Ole): What do we do if the first point has multiple
     
    531533                    v_prev = [1.0, 0.0]
    532534                   
    533                    
     535
    534536                # Choose candidate with minimum angle   
    535537                minimum_angle = 2*pi
    536538                for pc in candidate_list:
    537                     vc = pc-p0  # Candidate vector
    538                     # Normalise
    539                     #vc /= sqrt(sum(vc**2))                     
     539
    540540                   
     541                    vc = pc-p0  # Candidate vector (from p0 to candidate pt)
     542                   
    541543                    # Angle between each candidate and the previous vector
    542544                    # in [-pi, pi]
    543 
    544545                    ac = angle(vc, v_prev)
    545                     if ac > pi: ac = pi-ac
     546                    if ac > pi:
     547                        # Give preference to angles on the right hand side of v_prev
     548                        #print 'pc = %s, changing angle from %f to %f' %(pc, ac*180/pi, (ac-2*pi)*180/pi)
     549                        ac = ac-2*pi
    546550
    547551                    # take the minimal angle corresponding to the rightmost vector
     
    553557                if verbose is True:
    554558                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
    555                
     559                   
    556560            else:
    557561                p1 = candidate_list[0]
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r3684 r3688  
    910910
    911911
    912     def xtest_boundary_polygon_VI(self):
     912    def test_boundary_polygon_VI(self):
    913913        """test_boundary_polygon_VI(self)
    914914
     
    942942                  [  75735.4765625 ,  23762.00585938],
    943943                  [  52341.70703125,  38563.39453125]]
     944
     945        ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
    944946
    945947        triangles = [[19, 0,15],
     
    969971        Pref = mesh.get_boundary_polygon()
    970972
    971         plot_polygons([ensure_numeric(Pref)], 'goodP')
     973        #plot_polygons([ensure_numeric(Pref)], 'goodP')
    972974
    973975        for p in points:
     
    10441046                  [  35406.3359375 ,  79332.9140625 ]]
    10451047
    1046         #points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
     1048        scaled_points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation
    10471049
    10481050        triangles = [[ 0, 1, 2],
     
    10681070                     [60,61,62]]
    10691071
     1072
     1073        # First use scaled points for ease of debugging
     1074        mesh = Mesh(scaled_points, triangles)
     1075        mesh.check_integrity()
     1076        P = mesh.get_boundary_polygon()
     1077
     1078        for p in scaled_points:
     1079            assert is_inside_polygon(p, P)           
     1080
     1081        # Then use original points and test       
    10701082        mesh = Mesh(points, triangles)
    10711083        mesh.check_integrity()
    10721084        P = mesh.get_boundary_polygon()
    10731085
    1074         plot_polygons([ensure_numeric(P)], 'badP')
    1075         print P
    1076 
    10771086        for p in points:
    1078             #assert is_inside_polygon(p, P)           
    1079             if not is_inside_polygon(p, P):
    1080                 print 'Point %s is not in P' %(p)
    1081             else:
    1082                 print 'Point %s OK' %(p)
    1083 
    1084         #for i in range(len(Pref)):
    1085         #    print P[i], Pref[i]
    1086         #
    1087         #print ensure_numeric(P) - ensure_numeric(Pref)
    1088 
    1089 
     1087            assert is_inside_polygon(p, P)           
     1088
     1089        assert allclose(P, Pref)   
    10901090
    10911091    def test_lone_vertices(self):
Note: See TracChangeset for help on using the changeset viewer.