source: anuga_core/source/anuga/fit_interpolate/search_functions.py @ 3941

Last change on this file since 3941 was 3865, checked in by ole, 17 years ago

Ooops - one of them should be machine precision. Sorry.

File size: 3.9 KB
Line 
1"""
2General functions used in fit and interpolate.
3
4   Ole Nielsen, Stephen Roberts, Duncan Gray
5   Geoscience Australia, 2006.
6
7"""
8from Numeric import dot
9
10from numerical_tools import get_machine_precision
11
12def search_tree_of_vertices(root, mesh, x):
13    """
14    Find the triangle (element) that the point x is in.
15   
16    root: A quad tree of the vertices
17    Return the associated sigma and k values
18    (and if the element was found) .
19    """
20    #Find triangle containing x:
21    element_found = False
22
23    # This will be returned if element_found = False
24    sigma2 = -10.0
25    sigma0 = -10.0
26    sigma1 = -10.0
27    k = -10.0
28           
29    #Find vertices near x
30    candidate_vertices = root.search(x[0], x[1])
31    is_more_elements = True
32
33    element_found, sigma0, sigma1, sigma2, k = \
34                   _search_triangles_of_vertices(mesh,
35                                                 candidate_vertices, x)
36    while not element_found and is_more_elements:
37        candidate_vertices, branch = root.expand_search()
38        if branch == []:
39            # Searching all the verts from the root cell that haven't
40            # been searched.  This is the last try
41            element_found, sigma0, sigma1, sigma2, k = \
42                           _search_triangles_of_vertices(mesh,
43                                                         candidate_vertices, x)
44            is_more_elements = False
45        else:
46            element_found, sigma0, sigma1, sigma2, k = \
47                       _search_triangles_of_vertices(mesh,
48                                                     candidate_vertices, x)
49
50    return element_found, sigma0, sigma1, sigma2, k
51
52def _search_triangles_of_vertices(mesh, candidate_vertices, x):
53    #Find triangle containing x:
54    element_found = False
55
56    # This will be returned if element_found = False
57    sigma2 = -10.0
58    sigma0 = -10.0
59    sigma1 = -10.0
60    k = -10.0
61    #print "*$* candidate_vertices", candidate_vertices
62    #For all vertices in same cell as point x
63    for v in candidate_vertices:
64        #FIXME (DSG-DSG): this catches verts with no triangle.
65        #Currently pmesh is producing these.
66        #this should be stopped,
67        if mesh.vertexlist[v] is None:
68            continue
69        #for each triangle id (k) which has v as a vertex
70        for k, _ in mesh.vertexlist[v]:
71            #Get the three vertex_points of candidate triangle
72            xi0 = mesh.get_vertex_coordinate(k, 0)
73            xi1 = mesh.get_vertex_coordinate(k, 1)
74            xi2 = mesh.get_vertex_coordinate(k, 2)
75
76            #Get the three normals
77            n0 = mesh.get_normal(k, 0)
78            n1 = mesh.get_normal(k, 1)
79            n2 = mesh.get_normal(k, 2)
80           
81            #Compute interpolation
82            sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
83            sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
84            sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
85
86            # Integrity check - machine precision is too hard.
87            #epsilon = get_machine_precision()           
88            # Use single precision as we used to here
89            epsilon = 1.0e-6
90            msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
91                  %(abs(sigma0+sigma1+sigma2-1), epsilon)
92            assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon, msg
93
94            #Check that this triangle contains the data point
95           
96            # Sigmas are allowed to get negative within
97            # machine precision on some machines (e.g nautilus)
98            epsilon = get_machine_precision()           
99            if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
100                element_found = True
101                break
102
103        if element_found is True:
104            #Don't look for any other triangle
105            break
106    return element_found, sigma0, sigma1, sigma2, k
107
108
Note: See TracBrowser for help on using the repository browser.