source: inundation/ga/storm_surge/pyvolution/pmesh2domain.py @ 905

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


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


def pmesh_to_domain_instance(fileName, DomainClass, setting_function = None):
    """
    """
    import sys
    from domain import Domain

    vertex_coordinates, volumes, marker_dict, vertex_quantity_dict ,tagged_elements_dict = \
                        pmesh_to_domain(fileName,
                                                   setting_function=setting_function)
        
    assert issubclass(DomainClass, Domain), "DomainClass is not a subclass of Domain."


    
    domain = DomainClass(vertex_coordinates, volumes, marker_dict,
                         tagged_elements = tagged_elements_dict )


    
    #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(fileName, setting_function = 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 marker list: [S1Marker, S2Marker, ...] (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 mesh_file_to_mesh_dictionary
    
    mesh_dict = mesh_file_to_mesh_dictionary(fileName)
    #print "mesh_dict",mesh_dict
    vertex_coordinates = mesh_dict['generatedpointlist']
    volumes = mesh_dict['generatedtrianglelist']

    #if setting_function:
    #    if not type(setting_function) is ListType:
    #        setting_function = [setting_function]
    #    for funct in setting_function:
    #        mesh_dict = funct(mesh_dict, vertices = mesh_vertices,
    #                        volumes = volumes)


    vertex_quantity_dict = {}
    point_atts = transpose(mesh_dict['generatedpointattributelist'])
    point_titles  = mesh_dict['generatedpointattributetitlelist']
    #print "point_titles",point_titles
    for quantity, value_vector in map (None, point_titles, point_atts):
        vertex_quantity_dict[quantity] = value_vector
    marker_dict = pmesh_dict_to_marker_dict(mesh_dict)
    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    return vertex_coordinates, volumes, marker_dict, vertex_quantity_dict,tagged_elements_dict



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['generatedtriangleattributelist']
    #print "tri_atts", tri_atts
    tagged_elements = {}
    for tri_att_index in range(len(tri_atts)):
        tagged_elements.setdefault(tri_atts[tri_att_index][0],[]).append(tri_att_index)
    #print "DSG pm2do tagged_elements", tagged_elements
    return tagged_elements

#FIXME: The issue is whether this format should be stored in the tsh file
#instead of having to be created here?

#This information is pyvolution focused, not mesh generation focused.
#This is an appropriate place for the info to be created. -DSG

#FIXME: Another issue is that the tsh file stores consecutive
#indices explicitly. This is really redundant.
#Suggest looking at obj and our own sww format and also consider
#using netCDF.

# It is redundant.  It was the format originally decided on.
# I'm happy for it to change. -DSG

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

    return marker_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