"""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, tag_dict, vertex_quantity_dict \ ,tagged_elements_dict, geo_reference = \ pmesh_to_domain(fileName, setting_function=setting_function) assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain." domain = DomainClass(vertex_coordinates, volumes, 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(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 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 mesh_file_to_mesh_dictionary mesh_dict = mesh_file_to_mesh_dictionary(fileName) #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 #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_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