Changeset 6544


Ignore:
Timestamp:
Mar 18, 2009, 2:34:46 PM (15 years ago)
Author:
ole
Message:

Optimised fitting a bit more and cleaned up in search_functions.

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r6541 r6544  
    1919search_more_cells_time = initial_search_value
    2020
    21 #FIXME test what happens if a
    22 LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]),
    23                         num.array([max_float,max_float]),
    24                         num.array([max_float,max_float])),
    25                        (num.array([1,1], num.Int),      #array default#
    26                         num.array([0,0], num.Int),      #array default#
    27                         num.array([-1.1,-1.1]))]]]
     21LAST_TRIANGLE = [[-10,
     22                   (num.array([[max_float,max_float],
     23                               [max_float,max_float],
     24                               [max_float,max_float]]),
     25                    (num.array([1,1], num.Int),      #array default#
     26                     num.array([0,0], num.Int),      #array default#
     27                     num.array([-1.1,-1.1])))]]
    2828
    2929def search_tree_of_vertices(root, mesh, x):
     
    4848    global search_more_cells_time
    4949
    50     #Find triangle containing x:
    51     element_found = False
    52 
    5350    # This will be returned if element_found = False
     51    element_found = False   
    5452    sigma2 = -10.0
    5553    sigma0 = -10.0
     
    6159        element_found, sigma0, sigma1, sigma2, k = \
    6260            _search_triangles_of_vertices(mesh, last_triangle, x)
     61           
    6362    except:
     63        print 'This should never happen:', last_triangle
    6464        element_found = False
    6565                   
    66     #print "last_triangle", last_triangle
    6766    if element_found is True:
    68         #print "last_triangle", last_triangle
    6967        return element_found, sigma0, sigma1, sigma2, k
    7068
    71     # This was only slightly faster than just checking the
    72     # last triangle and it significantly slowed down
    73     # non-gridded fitting
    74     # If the last element was a dud, search its neighbours
    75     #print "last_triangle[0][0]", last_triangle[0][0]
    76     #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0])
    77     #print "neighbours", neighbours
    78     #neighbours = []
    79   #   for k in neighbours:
    80 #         if k >= 0:
    81 #             tri = mesh.get_vertex_coordinates(k,
    82 #                                                    absolute=True)
    83 #             n0 = mesh.get_normal(k, 0)
    84 #             n1 = mesh.get_normal(k, 1)
    85 #             n2 = mesh.get_normal(k, 2)
    86 #             triangle =[[k,(tri, (n0, n1, n2))]]
    87 #             element_found, sigma0, sigma1, sigma2, k = \
    88 #                            _search_triangles_of_vertices(mesh,
    89 #                                                          triangle, x)
    90 #             if element_found is True:
    91 #                 return element_found, sigma0, sigma1, sigma2, k
    92            
    93     #t0 = time.time()
     69   
    9470    # Get triangles in the cell that the point is in.
    9571    # Triangle is a list, first element triangle_id,
     
    9773    triangles = root.search(x[0], x[1])
    9874    is_more_elements = True
    99    
    10075    element_found, sigma0, sigma1, sigma2, k = \
    10176                   _search_triangles_of_vertices(mesh,
     
    133108    global last_triangle
    134109   
    135     # these statments are needed if triangles is empty
    136     #Find triangle containing x:
     110    # These statments are needed if triangles is empty
    137111    element_found = False
    138 
    139     # This will be returned if element_found = False
    140112    sigma2 = -10.0
    141113    sigma0 = -10.0
    142114    sigma1 = -10.0
    143115    k = -10
    144    
    145     #For all vertices in same cell as point x
     116
     117    # For all vertices in same cell as point x
    146118    for k, tri_verts_norms in triangles:
    147119        tri = tri_verts_norms[0]
     
    150122        # tri is a list of verts (x, y), representing a tringle
    151123        # Find triangle that contains x (if any) and interpolate
    152         if is_inside_triangle(x, tri, closed=True):
     124       
     125        # Input check disabled to speed things up.
     126        if is_inside_triangle(x, tri,
     127                              closed=True,
     128                              check_inputs=False):
     129           
    153130            n0, n1, n2 = tri_verts_norms[1]       
    154             element_found, sigma0, sigma1, sigma2 =\
    155                 find_triangle_compute_interpolation(tri, n0, n1, n2, x)
    156         else:
    157             element_found = False
     131            sigma0, sigma1, sigma2 =\
     132                compute_interpolation_values(tri, n0, n1, n2, x)
     133               
     134            element_found = True   
     135
     136            # Don't look for any other triangles in the triangle list
     137            last_triangle = [[k, tri_verts_norms]]
     138            break           
    158139
    159140           
    160         if element_found is True:
    161             # Don't look for any other triangles in the triangle list
    162             last_triangle = [[k,tri_verts_norms]]
    163             break
    164141    return element_found, sigma0, sigma1, sigma2, k
    165142
    166 
    167143           
    168 def find_triangle_compute_interpolation(triangle, n0, n1, n2, x):
    169     """Compute linear interpolation of point x and triangle k in mesh.
    170     It is assumed that x belongs to triangle k.max_float
     144def compute_interpolation_values(triangle, n0, n1, n2, x):
     145    """Compute linear interpolation of point x and triangle.
     146   
     147    n0, n1, n2 are normal to the tree edges.
    171148    """
    172149
     
    174151    xi0, xi1, xi2 = triangle
    175152
    176     # this is where we can call some fast c code.
    177      
    178     # Integrity check - machine precision is too hard
    179     # so we use hardwired single precision
    180     #epsilon = 1.0e-6
     153    sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
     154    sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
     155    sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
     156
     157    return sigma0, sigma1, sigma2
     158
    181159   
    182     #if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
    183     #    return False,0,0,0
    184     #if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:
    185     #    return False,0,0,0
    186     #if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:
    187     #    return False,0,0,0
    188     #if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:
    189     #    return False,0,0,0
    190    
    191     # machine precision on some machines (e.g. nautilus)
    192     epsilon = get_machine_precision() * 2
    193    
    194     # Compute interpolation - return as soon as possible
    195     #  print "(xi0-xi1)", (xi0-xi1)
    196     # print "n0", n0
    197     # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)
    198    
    199     sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)
    200     #if sigma0 < -epsilon:
    201         #print 'sigma0', sigma0
    202         #sigma0 = 0.0
    203         #return False,0,0,0
    204     sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)
    205     #if sigma1 < -epsilon:
    206         #print 'sigma1', sigma1   
    207         #sigma1 = 0.0       
    208         #return False,0,0,0
    209     sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)
    210     #if sigma2 < -epsilon:
    211         #print 'sigma2', sigma2   
    212         #sigma2 = 0.0       
    213         #return False,0,0,0
    214    
    215     # epsilon = 1.0e-6
    216     # we want to speed this up, so don't do assertions
    217     #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
    218     #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
    219     #      %(delta, epsilon)
    220     #assert delta < epsilon, msg
    221 
    222 
    223     # Check that this triangle contains the data point
    224     # Sigmas are allowed to get negative within
    225     # machine precision on some machines (e.g. nautilus)
    226     #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    227     #    element_found = True
    228     #else:
    229     #    element_found = False
    230     return True, sigma0, sigma1, sigma2
    231 
    232160def set_last_triangle():
    233161    global last_triangle
  • anuga_core/source/anuga/fit_interpolate/test_search_functions.py

    r6541 r6544  
    55from search_functions import search_tree_of_vertices, set_last_triangle
    66from search_functions import _search_triangles_of_vertices
    7 from search_functions import find_triangle_compute_interpolation
     7from search_functions import compute_interpolation_values
    88
    99from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
     
    1111from anuga.utilities.polygon import is_inside_polygon
    1212from anuga.utilities.quad import build_quadtree, Cell
     13from anuga.utilities.numerical_tools import ensure_numeric
    1314from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle   
    1415
     
    5253
    5354        x = [0.2, 0.7]
    54         found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     55        found, s0, s1, s2, k = search_tree_of_vertices(root,
     56                                                       mesh,
     57                                                       ensure_numeric(x))
    5558        assert k == 1 # Triangle one
    5659        assert found is True
    5760
    5861    def test_bigger(self):
    59         """test_larger mesh
     62        """test_bigger
     63       
     64        test larger mesh
    6065        """
    6166
     
    7277                  [0.1,0.9], [0.4,0.6], [0.9,0.1],
    7378                  [10, 3]]:
    74                
    75             found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     79
     80           
     81            found, s0, s1, s2, k = search_tree_of_vertices(root,
     82                                                           mesh,
     83                                                           ensure_numeric(x))
    7684
    7785            if k >= 0:
     
    95103        mesh.check_integrity()
    96104
    97        
    98105        for m in range(8):
    99106            root = build_quadtree(mesh, max_points_per_cell = m)
     
    105112                      [10, 3]]:
    106113               
    107                 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x)
     114                found, s0, s1, s2, k = search_tree_of_vertices(root,
     115                                                               mesh,
     116                                                               x)
    108117
    109118                if k >= 0:
    110119                    V = mesh.get_vertex_coordinates(k) # nodes for triangle k
     120
     121                    assert is_inside_triangle(x, V, closed=True)
    111122                    assert is_inside_polygon(x, V)
     123                   
    112124                    assert found is True
    113125                else:
    114126                    assert found is False               
    115127
    116                
     128           
     129                if m == k == 0: return   
    117130
    118131    def test_underlying_function(self):
     
    127140
    128141        # One point
    129         x = [0.5, 0.5]
     142        x = ensure_numeric([0.5, 0.5])
    130143        candidate_vertices = root.search(x[0], x[1])
    131144
    132         #print x, candidate_vertices
     145        # print x, candidate_vertices
    133146        found, sigma0, sigma1, sigma2, k = \
    134147               _search_triangles_of_vertices(mesh,
     
    156169                   _search_triangles_of_vertices(mesh,
    157170                                                 candidate_vertices,
    158                                                  x)
     171                                                 ensure_numeric(x))
    159172            if k >= 0:
    160173                V = mesh.get_vertex_coordinates(k) # nodes for triangle k
     
    211224       
    212225       
    213     def test_triangle_compute_interpolation(self):
    214         """test_triangle_compute_interpolation
    215        
    216         Test that triangle can be found if point is inside it
     226    def test_compute_interpolation_values(self):
     227        """test_compute_interpolation_values
     228       
     229        Test that interpolation values are correc
     230       
     231        This test used to check element_found as output from
     232        find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
     233        and that failed before 18th March 2009.
     234       
     235        Now this function no longer returns this flag, so the test
     236        is merely checknig the sigmas.
     237       
    217238        """
    218239       
     
    232253                                 closed=True, verbose=False)
    233254        assert is_inside_triangle(x, triangle,
    234                                   closed=True, verbose=False)                                 
    235        
    236         element_found, sigma0, sigma1, sigma2 = \
    237             find_triangle_compute_interpolation(triangle, n0, n1, n2, x)
     255                                  closed=True, verbose=False)
     256       
     257        sigma0, sigma1, sigma2 = \
     258            compute_interpolation_values(triangle, n0, n1, n2, x)
    238259           
    239260        msg = 'Point which is clearly inside triangle was not found'
    240         assert element_found is True, msg
    241        
    242         #print sigma0, sigma1, sigma2
    243         #print sigma0 + sigma1 + sigma2
    244                
    245         assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6
     261        assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg
    246262
    247263#-------------------------------------------------------------
  • anuga_core/source/anuga/utilities/polygon.py

    r6541 r6544  
    232232    """
    233233
     234    triangle = ensure_numeric(triangle)       
     235    point = ensure_numeric(point, num.Float)   
     236   
    234237    if check_inputs is True:
    235         triangle = ensure_numeric(triangle)   
    236 
    237         point = ensure_numeric(point, num.Float) # Convert point to array of points
    238238        msg = 'is_inside_triangle must be invoked with one point only'
    239239        assert num.allclose(point.shape, [2]), msg
    240240   
    241        
     241   
     242    # Use C-implementation
     243    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
     244   
     245   
     246
     247    # FIXME (Ole): The rest of this function has been made
     248    # obsolete by the C extension.
     249   
    242250    # Quickly reject points that are clearly outside
    243251    if point[0] < min(triangle[:,0]): return False
     
    245253    if point[1] < min(triangle[:,1]): return False
    246254    if point[1] > max(triangle[:,1]): return False       
    247 
    248     # Use C-implementation
    249     return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
    250    
    251    
    252 
    253     # FIXME (Ole): The rest of this function has been made
    254     # obsolete by the C extension.
    255255       
    256256    # Start search   
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r6535 r6544  
    254254  double a00, a10, a01, a11, b0, b1;
    255255  double denom, alpha, beta;
     256 
     257  double x, y; // Point coordinates
    256258  int i, j, res;
    257259
     260  x = point[0];
     261  y = point[1];
     262 
     263  // Quickly reject points that are clearly outside
     264  if ((x < triangle[0]) &&
     265      (x < triangle[2]) &&
     266      (x < triangle[4])) return 0;       
     267     
     268  if ((x > triangle[0]) &&
     269      (x > triangle[2]) &&
     270      (x > triangle[4])) return 0;             
     271 
     272  if ((y < triangle[1]) &&
     273      (y < triangle[3]) &&
     274      (y < triangle[5])) return 0;       
     275     
     276  if ((y > triangle[1]) &&
     277      (y > triangle[3]) &&
     278      (y > triangle[5])) return 0;             
     279 
     280 
    258281  // v0 = C-A
    259282  v0x = triangle[4]-triangle[0];
     
    274297  if (fabs(denom) > 0.0) {
    275298    // v = point-A 
    276     vx = point[0] - triangle[0];
    277     vy = point[1] - triangle[1];     
     299    vx = x - triangle[0];
     300    vy = y - triangle[1];     
    278301   
    279302    b0 = v0x*vx + v0y*vy; // innerproduct(v0, v)       
     
    291314    for (i=0; i<3; i++) {
    292315      j = (i+1) % 3; // Circular index into triangle array
    293       res = __point_on_line(point[0], point[1],
     316      res = __point_on_line(x, y,
    294317                            triangle[2*i], triangle[2*i+1],
    295318                            triangle[2*j], triangle[2*j+1],                         
  • anuga_core/source/anuga/utilities/test_polygon.py

    r6535 r6544  
    18961896        assert res is False
    18971897       
    1898        
     1898
     1899        res = is_inside_triangle([10, 3],
     1900                                 [[0.1, 0.],
     1901                                  [0.1, 0.08333333],
     1902                                  [0., 0.]])
     1903        assert res is False
     1904               
    18991905               
    19001906#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.