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

"""Class mesh2domain - 2D triangular domains for finite-volume computations of
   the shallow water wave equation


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


def pmesh_to_domain(fileName, tags = None, setting_function = None):
    """
    """
    #FIXME:  The current design doesn't seem to accomadate tags and
    # setting_functions being passed into the domain at this point.

    #FIXME: Plus the names of the functions are no longer appropriate,
    #since domain objects aren't being returned.
    
    import sys
    from config import pmesh_filename
    sys.path.append(pmesh_filename)
    from load_mesh.loadASCII import import_trianglulation

    try:
        meshdic = import_trianglulation(fileName)
    except IOError, e:        
        msg = 'Could not open file %s ' %fileName
        raise IOError, msg 

        
    return pmesh_dictionary_to_domain(meshdic,
                                      setting_function = setting_function)



def pmesh_dictionary_to_domain(meshdic, setting_function = None):
    """
    convert a pmesh dictionary to a list of Volumes.
    Also, return a list of triangles which have boundary tags
    meshdic 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 floats)
    """

    vertex_coordinates = meshdic['generatedpointlist']
    volumes = meshdic['generatedtrianglelist']

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

    return vertex_coordinates, volumes, marker_dict




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

#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.

   
def pmesh_dict_to_marker_dict(triangles,meshdic,Vertices):
    """ Convert the pmesh dictionary (meshdic) description of boundary tags
    to a dictionary of markers, indexed with volume id and face number.
    """
    sides = calc_sides(triangles)
    marker_dict = {}
    for seg, marker in map(None,meshdic['generatedsegmentlist'],
                           meshdic['generatedsegmentmarkerlist']):
        v1 = seg[0]
        v2 = seg[1]
        for key in [(v1,v2),(v2,v1)]:
            if sides.has_key(key):
                #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