Changeset 4658


Ignore:
Timestamp:
Aug 6, 2007, 5:00:03 PM (10 years ago)
Author:
duncan
Message:

minor optimising changes - minimal speed up to fitting

Location:
anuga_core/source/anuga/fit_interpolate
Files:
2 edited

Legend:

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

    r4651 r4658  
    136136                pfile = open(profile_file, "w")
    137137                sys.stdout = pfile
    138                 s = S.sort_stats('cumulative').print_stats(30)
     138                s = S.sort_stats('cumulative').print_stats(60)
    139139                sys.stdout = saveout
    140140                pfile.close()
     
    164164                pfile = open(profile_file, "w")
    165165                sys.stdout = pfile
    166                 s = S.sort_stats('cumulative').print_stats(30)
     166                s = S.sort_stats('cumulative').print_stats(60)
    167167                sys.stdout = saveout
    168168                pfile.close()
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r4655 r4658  
    7171    fit and interpolate.
    7272    """
    73    
     73
     74    # these statments are needed if triangles is empty
    7475    #Find triangle containing x:
    7576    element_found = False
     
    111112    epsilon = 1.0e-6
    112113   
    113     xmax = max(xi0[0], xi1[0], xi2[0])
    114     xmin = min(xi0[0], xi1[0], xi2[0])
    115     ymax = max(xi0[1], xi1[1], xi2[1])
    116     ymin = min(xi0[1], xi1[1], xi2[1])
    117 
    118    
    119     if  x[0] > xmax + epsilon:
     114    if  x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:
    120115        return False,0,0,0
    121     if  x[0] < xmin - epsilon:
     116    if  x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:
    122117        return False,0,0,0
    123     if  x[1] > ymax + epsilon:
     118    if  x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:
    124119        return False,0,0,0
    125     if  x[1] < ymin - epsilon:
     120    if  x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:
    126121        return False,0,0,0
    127122   
    128     # Get the three normals
    129     #n0 = norms[0] 
    130     #n1 = norms[1]
    131     #n2 = norms[2]
    132        
    133     # Compute interpolation
     123    # machine precision on some machines (e.g. nautilus)
     124    epsilon = get_machine_precision() * 2
     125   
     126    # Compute interpolation - return as soon as possible
     127   
     128    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
     129    if sigma0 < -epsilon:
     130        return False,0,0,0
     131    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
     132    if sigma1 < -epsilon:
     133        return False,0,0,0
    134134    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
    135     sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
    136     sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
    137 
    138     delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
    139     msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
    140           %(delta, epsilon)
    141     assert delta < epsilon, msg
     135    if sigma2 < -epsilon:
     136        return False,0,0,0
     137   
     138    # epsilon = 1.0e-6
     139    # we want to speed this up, so don't do assertions
     140    #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
     141    #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
     142    #      %(delta, epsilon)
     143    #assert delta < epsilon, msg
    142144
    143145
     
    145147    # Sigmas are allowed to get negative within
    146148    # machine precision on some machines (e.g. nautilus)
    147     epsilon = get_machine_precision() * 2
    148     if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    149         element_found = True
    150     else:
    151         element_found = False
    152     return element_found, sigma0, sigma1, sigma2
     149    #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
     150    #    element_found = True
     151    #else:
     152    #    element_found = False
     153    return True, sigma0, sigma1, sigma2
Note: See TracChangeset for help on using the changeset viewer.