source: branches/inundation-numpy-branch/pyvolution/pmesh2domain.py @ 7248

"""Class pmesh2domain - Converting .tsh files to doamains


   Copyright 2004
   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
   Geoscience Australia
"""

def pmesh_instance_to_domain_instance(mesh,
                                      DomainClass):
    """
    """
    import sys
    from domain import Domain

    vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
                        ,tagged_elements_dict, geo_reference = \
                        _pmesh_to_domain(mesh_instance=mesh)


    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."


    domain = DomainClass(coordinates = vertex_coordinates,
                         vertices = vertices,
                         boundary = tag_dict,
                         tagged_elements = tagged_elements_dict,
                         geo_reference = geo_reference )

    # set the water stage to be the elevation
    if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']

    domain.set_quantity_vertices_dict(vertex_quantity_dict)
    #print "vertex_quantity_dict",vertex_quantity_dict
    return domain



def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False):
    """
    Converts a mesh file(.tsh or .msh), to a Domain instance.

    file_name is the name of the mesh file to convert, including the extension

    DomainClass is the Class that will be returned.
    It must be a subclass of Domain, with the same interface as domain.

    use_cache: True means that caching is attempted for the computed domain.    
    """

    if use_cache is True:
        from caching import cache
        result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
                       dependencies = [file_name],                     
                       verbose = verbose)

    else:
        result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))        
        
    return result




def _pmesh_to_domain_instance(file_name, DomainClass):
    """
    Converts a mesh file(.tsh or .msh), to a Domain instance.

    Internal function. See public interface pmesh_to_domain_instance for details
    """
    
    import sys
    from domain import Domain

    vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
                        ,tagged_elements_dict, geo_reference = \
                        _pmesh_to_domain(file_name=file_name)


    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."


    domain = DomainClass(coordinates = vertex_coordinates,
                         vertices = vertices,
                         boundary = tag_dict,
                         tagged_elements = tagged_elements_dict,
                         geo_reference = geo_reference )



    #FIXME (Ole): Is this really the right place to apply the a default
    #value specific to the shallow water wave equation?
    #The 'assert' above indicates that any subclass of Domain is acceptable.
    #Suggestion - module shallow_water.py will eventually take care of this
    #(when I get around to it) so it should be removed from here.

    # This doesn't work on the domain instance.
    # This is still needed so -ve elevations don't cuase 'lakes'
    # The fixme we discussed was to only create a quantity when its values
    #are set.
    # I think that's the way to go still

    # set the water stage to be the elevation
    if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']

    domain.set_quantity_vertices_dict(vertex_quantity_dict)
    #print "vertex_quantity_dict",vertex_quantity_dict
    return domain



def _pmesh_to_domain(file_name=None,
                    mesh_instance=None):
    """
    convert a pmesh dictionary to a list of Volumes.
    Also, return a list of triangles which have boundary tags
    mesh_dict structure;
    generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
    generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
    generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
    generated segment tag list: [S1Tag, S2Tag, ...] (list of ints)
    triangle list:  [(point1, point2, point3),(p5,p4, p1),...] (Tuples of integers)
    triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
    triangle attribute list: [T1att, T2att, ...] (list of strings)
    """

    from Numeric import transpose
    from load_mesh.loadASCII import import_mesh_file

    if file_name is None:
        mesh_dict = mesh_instance.Mesh2IODict()
    else:
        mesh_dict = import_mesh_file(file_name)
    #print "mesh_dict",mesh_dict
    vertex_coordinates = mesh_dict['vertices']
    volumes = mesh_dict['triangles']
    vertex_quantity_dict = {}
    point_atts = transpose(mesh_dict['vertex_attributes'])
    point_titles  = mesh_dict['vertex_attribute_titles']
    geo_reference  = mesh_dict['geo_reference']
    for quantity, value_vector in map (None, point_titles, point_atts):
        vertex_quantity_dict[quantity] = value_vector
    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference



def build_tagged_elements_dictionary(mesh_dict):
    """Build the dictionary of element tags.
    tagged_elements is a dictionary of element arrays,
    keyed by tag:
    { (tag): [e1, e2, e3..] }
    """
    tri_atts = mesh_dict['triangle_tags']
    #print "tri_atts", tri_atts
    tagged_elements = {}
    for tri_att_index in range(len(tri_atts)):
        tagged_elements.setdefault(tri_atts[tri_att_index],[]).append(tri_att_index)
    #print "DSG pm2do tagged_elements", tagged_elements
    return tagged_elements

def pmesh_dict_to_tag_dict(mesh_dict):
    """ Convert the pmesh dictionary (mesh_dict) description of boundary tags
    to a dictionary of tags, indexed with volume id and face number.
    """
    triangles = mesh_dict['triangles']
    sides = calc_sides(triangles)
    tag_dict = {}
    for seg, tag in map(None, mesh_dict['segments'],
                        mesh_dict['segment_tags']):
        v1 = seg[0]
        v2 = seg[1]
        for key in [(v1,v2),(v2,v1)]:
            if sides.has_key(key) and tag <> "":
                #"" represents null.  Don't put these into the dictionary
                #this creates a dict of lists of faces, indexed by tag
                #tagged_edges.setdefault(tag,[]).append(sides[key])
                tag_dict[sides[key]] = tag

    return tag_dict


def calc_sides(triangles):
    #Build dictionary mapping from sides (2-tuple of points)
    #to left hand side neighbouring triangle
    sides = {}
    for id, triangle in enumerate(triangles):
        a = triangle[0]
        b = triangle[1]
        c = triangle[2]
        sides[a,b] = (id, 2) #(id, face)
        sides[b,c] = (id, 0) #(id, face)
        sides[c,a] = (id, 1) #(id, face)
    return sides