| 1 | """Class pmesh2domain - Converting .tsh files to doamains |
|---|
| 2 | |
|---|
| 3 | |
|---|
| 4 | Copyright 2004 |
|---|
| 5 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
|---|
| 6 | Geoscience Australia |
|---|
| 7 | """ |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | def pmesh_to_domain_instance(fileName, DomainClass,tags = None, setting_function = None): |
|---|
| 11 | """ |
|---|
| 12 | """ |
|---|
| 13 | #FIXME: The current design doesn't seem to accomadate tags |
|---|
| 14 | # being passed into the domain at this point. |
|---|
| 15 | |
|---|
| 16 | import sys |
|---|
| 17 | from domain import Domain |
|---|
| 18 | |
|---|
| 19 | vertex_coordinates, volumes, marker_dict, vertex_quantity_dict = \ |
|---|
| 20 | pmesh_to_domain(fileName, |
|---|
| 21 | setting_function=setting_function) |
|---|
| 22 | |
|---|
| 23 | assert issubclass(DomainClass, Domain), "DomainClass is not a subclass of Domain." |
|---|
| 24 | domain = DomainClass(vertex_coordinates, volumes, marker_dict) |
|---|
| 25 | |
|---|
| 26 | # set the water level to be the elevation |
|---|
| 27 | if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('level'): |
|---|
| 28 | vertex_quantity_dict['level'] = vertex_quantity_dict['elevation'] |
|---|
| 29 | domain.set_quantity_vertices_dict(vertex_quantity_dict) |
|---|
| 30 | #print "vertex_quantity_dict",vertex_quantity_dict |
|---|
| 31 | return domain |
|---|
| 32 | |
|---|
| 33 | |
|---|
| 34 | |
|---|
| 35 | def pmesh_to_domain(fileName, setting_function = None): |
|---|
| 36 | """ |
|---|
| 37 | convert a pmesh dictionary to a list of Volumes. |
|---|
| 38 | Also, return a list of triangles which have boundary tags |
|---|
| 39 | mesh_dict structure; |
|---|
| 40 | generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) |
|---|
| 41 | generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...] |
|---|
| 42 | generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) |
|---|
| 43 | generated segment marker list: [S1Marker, S2Marker, ...] (list of ints) |
|---|
| 44 | triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers) |
|---|
| 45 | triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor |
|---|
| 46 | triangle attribute list: [T1att, T2att, ...] (list of floats) |
|---|
| 47 | """ |
|---|
| 48 | |
|---|
| 49 | from Numeric import transpose |
|---|
| 50 | from load_mesh.loadASCII import mesh_file_to_mesh_dictionary |
|---|
| 51 | |
|---|
| 52 | mesh_dict = mesh_file_to_mesh_dictionary(fileName) |
|---|
| 53 | #print "mesh_dict",mesh_dict |
|---|
| 54 | vertex_coordinates = mesh_dict['generatedpointlist'] |
|---|
| 55 | volumes = mesh_dict['generatedtrianglelist'] |
|---|
| 56 | |
|---|
| 57 | #if setting_function: |
|---|
| 58 | # if not type(setting_function) is ListType: |
|---|
| 59 | # setting_function = [setting_function] |
|---|
| 60 | # for funct in setting_function: |
|---|
| 61 | # mesh_dict = funct(mesh_dict, vertices = mesh_vertices, |
|---|
| 62 | # volumes = volumes) |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | vertex_quantity_dict = {} |
|---|
| 66 | point_atts = transpose(mesh_dict['generatedpointattributelist']) |
|---|
| 67 | point_titles = mesh_dict['generatedpointattributetitlelist'] |
|---|
| 68 | #print "point_titles",point_titles |
|---|
| 69 | for quantity, value_vector in map (None, point_titles, point_atts): |
|---|
| 70 | vertex_quantity_dict[quantity] = value_vector |
|---|
| 71 | marker_dict = pmesh_dict_to_marker_dict(volumes, mesh_dict, |
|---|
| 72 | vertex_coordinates) |
|---|
| 73 | |
|---|
| 74 | return vertex_coordinates, volumes, marker_dict, vertex_quantity_dict |
|---|
| 75 | |
|---|
| 76 | |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | #FIXME: The issue is whether this format should be stored in the tsh file |
|---|
| 80 | #instead of having to be created here? |
|---|
| 81 | |
|---|
| 82 | #FIXME: Another issue is that the tsh file stores consecutive |
|---|
| 83 | #indices explicitly. This is really redundant. |
|---|
| 84 | #Suggest looking at obj and our own sww format and also consider |
|---|
| 85 | #using netCDF. |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | def pmesh_dict_to_marker_dict(triangles,mesh_dict,Vertices): |
|---|
| 89 | """ Convert the pmesh dictionary (mesh_dict) description of boundary tags |
|---|
| 90 | to a dictionary of markers, indexed with volume id and face number. |
|---|
| 91 | """ |
|---|
| 92 | sides = calc_sides(triangles) |
|---|
| 93 | marker_dict = {} |
|---|
| 94 | for seg, marker in map(None,mesh_dict['generatedsegmentlist'], |
|---|
| 95 | mesh_dict['generatedsegmentmarkerlist']): |
|---|
| 96 | v1 = seg[0] |
|---|
| 97 | v2 = seg[1] |
|---|
| 98 | for key in [(v1,v2),(v2,v1)]: |
|---|
| 99 | if sides.has_key(key): |
|---|
| 100 | #this creates a dict of lists of faces, indexed by marker |
|---|
| 101 | #tagged_edges.setdefault(marker,[]).append(sides[key]) |
|---|
| 102 | marker_dict[sides[key]] = marker |
|---|
| 103 | return marker_dict |
|---|
| 104 | |
|---|
| 105 | |
|---|
| 106 | def calc_sides(triangles): |
|---|
| 107 | #Build dictionary mapping from sides (2-tuple of points) |
|---|
| 108 | #to left hand side neighbouring triangle |
|---|
| 109 | sides = {} |
|---|
| 110 | for id, triangle in enumerate(triangles): |
|---|
| 111 | a = triangle[0] |
|---|
| 112 | b = triangle[1] |
|---|
| 113 | c = triangle[2] |
|---|
| 114 | sides[a,b] = (id, 2) #(id, face) |
|---|
| 115 | sides[b,c] = (id, 0) #(id, face) |
|---|
| 116 | sides[c,a] = (id, 1) #(id, face) |
|---|
| 117 | return sides |
|---|
| 118 | |
|---|
