Changeset 4599


Ignore:
Timestamp:
Jul 6, 2007, 12:28:53 PM (17 years ago)
Author:
ole
Message:

More work on search functions

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

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

    r4593 r4599  
    403403        Return list of triangle_ids, vertex_ids for specified node.
    404404        If node in None or absent, this information will be returned
    405         for all (full) nodes.
     405        for all (full) nodes in a list L where L[v] is the triangle
     406        list for node v.
    406407        """
    407408
     
    421422
    422423                triangle_list.append( (volume_id, vertex_id) )
     424
     425            triangle_list = array(triangle_list)   
    423426        else:
    424427            # Get info for all nodes recursively.
     
    426429            # working directly with the inverted triangle structure
    427430            for i in range(self.number_of_full_nodes):
     431               
    428432                L = self.get_triangles_and_vertices_per_node(node=i)
     433
    429434                triangle_list.append(L)
    430435
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r4593 r4599  
    166166        # One node
    167167        L = domain.get_triangles_and_vertices_per_node(node=2)
    168 
    169168        assert allclose(L[0], [0, 2])
    170169        assert allclose(L[1], [1, 1])
     
    173172        # All nodes
    174173        ALL = domain.get_triangles_and_vertices_per_node()
    175 
    176174        assert len(ALL) == 6
    177175        for i, Lref in enumerate(ALL):
     
    224222if __name__ == "__main__":
    225223    suite = unittest.makeSuite(Test_General_Mesh,'test')   
    226     #suite = unittest.makeSuite(Test_General_Mesh,'test_vertex_value_indices')
    227224    runner = unittest.TextTestRunner()
    228225    runner.run(suite)
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r4593 r4599  
    6060    return element_found, sigma0, sigma1, sigma2, k
    6161
     62
    6263def _search_triangles_of_vertices(mesh, candidate_vertices, x):
    6364    """Search for triangle containing x amongs candidate_vertices in mesh
     
    8283    #For all vertices in same cell as point x
    8384    for v in candidate_vertices:
     85       
    8486        #FIXME (DSG-DSG): this catches verts with no triangle.
    8587        #Currently pmesh is producing these.
     
    8991            continue
    9092       
    91         #for each triangle id (k) which has v as a vertex
    92         vertexlist = mesh.get_triangles_and_vertices_per_node(node=v)
    93         for k, _ in vertexlist:
    94             #Get the three vertex_points of candidate triangle k
    95             xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)           
     93        # Get all triangles which has v as a vertex
     94        # The list has elements (triangle, vertex), but only the
     95        # first component will be used here
     96        triangle_list = mesh.get_triangles_and_vertices_per_node(node=v)
    9697
    97             # Get the three normals (faster than using API)
    98             n0 = mesh.normals[k,0:2]
    99             n1 = mesh.normals[k,2:4]
    100             n2 = mesh.normals[k,4:6]           
     98        # Find triangle that contains x (if any) and interpolate
     99        element_found, sigma0, sigma1, sigma2, k =\
     100                       find_triangle_compute_interpolation(mesh,
     101                                                           triangle_list,
     102                                                           x)
    101103
    102            
    103             #Compute interpolation
    104             sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
    105             sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
    106             sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
     104        if element_found is True:
     105            # Don't look for any other triangle           
     106            break
    107107
    108             # Integrity check - machine precision is too hard
    109             # so we use hardwired single precision
    110             epsilon = 1.0e-6
    111             msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
    112                   %(abs(sigma0+sigma1+sigma2-1), epsilon)
    113             assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon, msg
    114 
    115 
    116             #Check that this triangle contains the data point
    117            
    118             # Sigmas are allowed to get negative within
    119             # machine precision on some machines (e.g. nautilus)
    120             epsilon = get_machine_precision() * 2
    121             if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    122                 element_found = True
    123                 break
    124            
    125         if element_found is True:
    126             # Don't look for any other triangle
    127             break
    128108       
    129109    return element_found, sigma0, sigma1, sigma2, k
    130110
    131111
     112           
     113def find_triangle_compute_interpolation(mesh, triangle_list, x):
     114    """Compute linear interpolation of point x and triangle k in mesh.
     115    It is assumed that x belongs to triangle k.
     116    """
     117
     118    element_found = False
     119    for k, _ in triangle_list:
     120        # Get the three vertex_points of candidate triangle k
     121        xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)           
     122       
     123        # Get the three normals
     124        n0 = mesh.get_normal(k, 0)
     125        n1 = mesh.get_normal(k, 1)
     126        n2 = mesh.get_normal(k, 2)           
     127       
     128       
     129        # Compute interpolation
     130        sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
     131        sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
     132        sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
     133
     134        # Integrity check - machine precision is too hard
     135        # so we use hardwired single precision
     136        epsilon = 1.0e-6
     137        delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
     138        msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
     139              %(delta, epsilon)
     140        assert delta < epsilon, msg
     141
     142
     143        # Check that this triangle contains the data point
     144        # Sigmas are allowed to get negative within
     145        # machine precision on some machines (e.g. nautilus)
     146        epsilon = get_machine_precision() * 2
     147        if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
     148            element_found = True
     149            break
     150           
     151    return element_found, sigma0, sigma1, sigma2, k
Note: See TracChangeset for help on using the changeset viewer.