[374] | 1 | """Class pmesh2domain - Converting .tsh files to doamains |
---|
[371] | 2 | |
---|
| 3 | |
---|
| 4 | Copyright 2004 |
---|
| 5 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
| 6 | Geoscience Australia |
---|
| 7 | """ |
---|
| 8 | |
---|
[2056] | 9 | # FXME (DSG-DSG): Combine with function pmesh_to_domain_instance |
---|
| 10 | def pmesh_instance_to_domain_instance(mesh, |
---|
| 11 | DomainClass, |
---|
| 12 | setting_function=None): |
---|
| 13 | """ |
---|
| 14 | """ |
---|
| 15 | import sys |
---|
| 16 | from domain import Domain |
---|
[371] | 17 | |
---|
[2056] | 18 | vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \ |
---|
| 19 | ,tagged_elements_dict, geo_reference = \ |
---|
| 20 | pmesh_to_domain(mesh_instance=mesh, |
---|
| 21 | setting_function=setting_function) |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain." |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | domain = DomainClass(coordinates = vertex_coordinates, |
---|
| 28 | vertices = vertices, |
---|
| 29 | boundary = tag_dict, |
---|
| 30 | tagged_elements = tagged_elements_dict, |
---|
| 31 | geo_reference = geo_reference ) |
---|
| 32 | |
---|
| 33 | # set the water stage to be the elevation |
---|
| 34 | if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): |
---|
| 35 | vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] |
---|
| 36 | |
---|
| 37 | domain.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
| 38 | #print "vertex_quantity_dict",vertex_quantity_dict |
---|
| 39 | return domain |
---|
| 40 | |
---|
[415] | 41 | def pmesh_to_domain_instance(fileName, DomainClass, setting_function = None): |
---|
[371] | 42 | """ |
---|
| 43 | """ |
---|
| 44 | import sys |
---|
[389] | 45 | from domain import Domain |
---|
[371] | 46 | |
---|
[1556] | 47 | vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \ |
---|
[1178] | 48 | ,tagged_elements_dict, geo_reference = \ |
---|
[2056] | 49 | pmesh_to_domain(fileName=fileName, |
---|
| 50 | setting_function=setting_function) |
---|
[1556] | 51 | |
---|
| 52 | |
---|
[1178] | 53 | assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain." |
---|
[636] | 54 | |
---|
| 55 | |
---|
[1556] | 56 | domain = DomainClass(coordinates = vertex_coordinates, |
---|
| 57 | vertices = vertices, |
---|
| 58 | boundary = tag_dict, |
---|
[1178] | 59 | tagged_elements = tagged_elements_dict, |
---|
| 60 | geo_reference = geo_reference ) |
---|
[371] | 61 | |
---|
[636] | 62 | |
---|
[1556] | 63 | |
---|
[636] | 64 | #FIXME (Ole): Is this really the right place to apply the a default |
---|
| 65 | #value specific to the shallow water wave equation? |
---|
[1556] | 66 | #The 'assert' above indicates that any subclass of Domain is acceptable. |
---|
[636] | 67 | #Suggestion - module shallow_water.py will eventually take care of this |
---|
[1556] | 68 | #(when I get around to it) so it should be removed from here. |
---|
| 69 | |
---|
[905] | 70 | # This doesn't work on the domain instance. |
---|
| 71 | # This is still needed so -ve elevations don't cuase 'lakes' |
---|
[1178] | 72 | # The fixme we discussed was to only create a quantity when its values |
---|
| 73 | #are set. |
---|
[905] | 74 | # I think that's the way to go still |
---|
[1556] | 75 | |
---|
[773] | 76 | # set the water stage to be the elevation |
---|
| 77 | if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): |
---|
| 78 | vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] |
---|
[1556] | 79 | |
---|
[389] | 80 | domain.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
| 81 | #print "vertex_quantity_dict",vertex_quantity_dict |
---|
| 82 | return domain |
---|
[371] | 83 | |
---|
| 84 | |
---|
[389] | 85 | |
---|
[2056] | 86 | def pmesh_to_domain(fileName=None, |
---|
| 87 | mesh_instance=None, |
---|
| 88 | setting_function=None): |
---|
[371] | 89 | """ |
---|
| 90 | convert a pmesh dictionary to a list of Volumes. |
---|
| 91 | Also, return a list of triangles which have boundary tags |
---|
[389] | 92 | mesh_dict structure; |
---|
[1556] | 93 | generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) |
---|
[371] | 94 | generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...] |
---|
| 95 | generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) |
---|
[969] | 96 | generated segment tag list: [S1Tag, S2Tag, ...] (list of ints) |
---|
[371] | 97 | triangle list: [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers) |
---|
| 98 | triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor |
---|
[415] | 99 | triangle attribute list: [T1att, T2att, ...] (list of strings) |
---|
[371] | 100 | """ |
---|
[1556] | 101 | |
---|
[389] | 102 | from Numeric import transpose |
---|
[1423] | 103 | from load_mesh.loadASCII import import_mesh_file |
---|
[1556] | 104 | |
---|
[2056] | 105 | if fileName is None: |
---|
| 106 | mesh_dict = mesh_instance.Mesh2IODict() |
---|
| 107 | else: |
---|
| 108 | mesh_dict = import_mesh_file(fileName) |
---|
[389] | 109 | #print "mesh_dict",mesh_dict |
---|
[968] | 110 | vertex_coordinates = mesh_dict['vertices'] |
---|
| 111 | volumes = mesh_dict['triangles'] |
---|
[389] | 112 | vertex_quantity_dict = {} |
---|
[968] | 113 | point_atts = transpose(mesh_dict['vertex_attributes']) |
---|
| 114 | point_titles = mesh_dict['vertex_attribute_titles'] |
---|
[1178] | 115 | geo_reference = mesh_dict['geo_reference'] |
---|
[389] | 116 | for quantity, value_vector in map (None, point_titles, point_atts): |
---|
| 117 | vertex_quantity_dict[quantity] = value_vector |
---|
[969] | 118 | tag_dict = pmesh_dict_to_tag_dict(mesh_dict) |
---|
[415] | 119 | tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict) |
---|
[1556] | 120 | return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference |
---|
[371] | 121 | |
---|
| 122 | |
---|
| 123 | |
---|
[415] | 124 | def build_tagged_elements_dictionary(mesh_dict): |
---|
| 125 | """Build the dictionary of element tags. |
---|
| 126 | tagged_elements is a dictionary of element arrays, |
---|
| 127 | keyed by tag: |
---|
| 128 | { (tag): [e1, e2, e3..] } |
---|
| 129 | """ |
---|
[969] | 130 | tri_atts = mesh_dict['triangle_tags'] |
---|
[415] | 131 | #print "tri_atts", tri_atts |
---|
| 132 | tagged_elements = {} |
---|
| 133 | for tri_att_index in range(len(tri_atts)): |
---|
[1183] | 134 | tagged_elements.setdefault(tri_atts[tri_att_index],[]).append(tri_att_index) |
---|
[415] | 135 | #print "DSG pm2do tagged_elements", tagged_elements |
---|
| 136 | return tagged_elements |
---|
[371] | 137 | |
---|
[969] | 138 | def pmesh_dict_to_tag_dict(mesh_dict): |
---|
[389] | 139 | """ Convert the pmesh dictionary (mesh_dict) description of boundary tags |
---|
[969] | 140 | to a dictionary of tags, indexed with volume id and face number. |
---|
[371] | 141 | """ |
---|
[968] | 142 | triangles = mesh_dict['triangles'] |
---|
[371] | 143 | sides = calc_sides(triangles) |
---|
[969] | 144 | tag_dict = {} |
---|
| 145 | for seg, tag in map(None,mesh_dict['segments'], |
---|
| 146 | mesh_dict['segment_tags']): |
---|
[371] | 147 | v1 = seg[0] |
---|
| 148 | v2 = seg[1] |
---|
| 149 | for key in [(v1,v2),(v2,v1)]: |
---|
[969] | 150 | if sides.has_key(key) and tag <> "": |
---|
[686] | 151 | #"" represents null. Don't put these into the dictionary |
---|
[969] | 152 | #this creates a dict of lists of faces, indexed by tag |
---|
| 153 | #tagged_edges.setdefault(tag,[]).append(sides[key]) |
---|
| 154 | tag_dict[sides[key]] = tag |
---|
[686] | 155 | |
---|
[969] | 156 | return tag_dict |
---|
[371] | 157 | |
---|
[1556] | 158 | |
---|
[371] | 159 | def calc_sides(triangles): |
---|
| 160 | #Build dictionary mapping from sides (2-tuple of points) |
---|
[1556] | 161 | #to left hand side neighbouring triangle |
---|
[371] | 162 | sides = {} |
---|
| 163 | for id, triangle in enumerate(triangles): |
---|
| 164 | a = triangle[0] |
---|
| 165 | b = triangle[1] |
---|
[1556] | 166 | c = triangle[2] |
---|
[371] | 167 | sides[a,b] = (id, 2) #(id, face) |
---|
| 168 | sides[b,c] = (id, 0) #(id, face) |
---|
| 169 | sides[c,a] = (id, 1) #(id, face) |
---|
[1556] | 170 | return sides |
---|