1 | ## Automatically adapted for numpy.oldnumeric Oct 28, 2008 by alter_code1.py |
---|
2 | |
---|
3 | """Class pmesh2domain - Converting .tsh files to doamains |
---|
4 | |
---|
5 | |
---|
6 | Copyright 2004 |
---|
7 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
---|
8 | Geoscience Australia |
---|
9 | """ |
---|
10 | |
---|
11 | import sys |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | def pmesh_instance_to_domain_instance(mesh, |
---|
16 | DomainClass): |
---|
17 | """ |
---|
18 | Convert a pmesh instance/object into a domain instance. |
---|
19 | |
---|
20 | Use pmesh_to_domain_instance to convert a mesh file to a domain instance. |
---|
21 | """ |
---|
22 | |
---|
23 | vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \ |
---|
24 | ,tagged_elements_dict, geo_reference = \ |
---|
25 | pmesh_to_domain(mesh_instance=mesh) |
---|
26 | |
---|
27 | # NOTE(Ole): This import cannot be at the module level due to mutual |
---|
28 | # dependency with domain.py |
---|
29 | from anuga.abstract_2d_finite_volumes.domain import Domain |
---|
30 | |
---|
31 | |
---|
32 | |
---|
33 | |
---|
34 | msg = 'The class %s is not a subclass of the generic domain class %s'\ |
---|
35 | %(DomainClass, Domain) |
---|
36 | assert issubclass(DomainClass, Domain), msg |
---|
37 | |
---|
38 | |
---|
39 | domain = DomainClass(coordinates = vertex_coordinates, |
---|
40 | vertices = vertices, |
---|
41 | boundary = tag_dict, |
---|
42 | tagged_elements = tagged_elements_dict, |
---|
43 | geo_reference = geo_reference ) |
---|
44 | |
---|
45 | # set the water stage to be the elevation |
---|
46 | if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): |
---|
47 | vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] |
---|
48 | |
---|
49 | domain.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
50 | #print "vertex_quantity_dict",vertex_quantity_dict |
---|
51 | return domain |
---|
52 | |
---|
53 | |
---|
54 | |
---|
55 | def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False): |
---|
56 | """ |
---|
57 | Converts a mesh file(.tsh or .msh), to a Domain instance. |
---|
58 | |
---|
59 | file_name is the name of the mesh file to convert, including the extension |
---|
60 | |
---|
61 | DomainClass is the Class that will be returned. |
---|
62 | It must be a subclass of Domain, with the same interface as domain. |
---|
63 | |
---|
64 | use_cache: True means that caching is attempted for the computed domain. |
---|
65 | """ |
---|
66 | |
---|
67 | if use_cache is True: |
---|
68 | from caching import cache |
---|
69 | result = cache(_pmesh_to_domain_instance, (file_name, DomainClass), |
---|
70 | dependencies = [file_name], |
---|
71 | verbose = verbose) |
---|
72 | |
---|
73 | else: |
---|
74 | result = apply(_pmesh_to_domain_instance, (file_name, DomainClass)) |
---|
75 | |
---|
76 | return result |
---|
77 | |
---|
78 | |
---|
79 | |
---|
80 | |
---|
81 | def _pmesh_to_domain_instance(file_name, DomainClass): |
---|
82 | """ |
---|
83 | Converts a mesh file(.tsh or .msh), to a Domain instance. |
---|
84 | |
---|
85 | Internal function. See public interface pmesh_to_domain_instance for details |
---|
86 | """ |
---|
87 | |
---|
88 | vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, \ |
---|
89 | tagged_elements_dict, geo_reference = \ |
---|
90 | pmesh_to_domain(file_name=file_name) |
---|
91 | |
---|
92 | |
---|
93 | # NOTE(Ole): This import cannot be at the module level due to mutual |
---|
94 | # dependency with domain.py |
---|
95 | from anuga.abstract_2d_finite_volumes.domain import Domain |
---|
96 | |
---|
97 | |
---|
98 | msg = 'The class %s is not a subclass of the generic domain class %s'\ |
---|
99 | %(DomainClass, Domain) |
---|
100 | assert issubclass(DomainClass, Domain), msg |
---|
101 | |
---|
102 | |
---|
103 | |
---|
104 | domain = DomainClass(coordinates = vertex_coordinates, |
---|
105 | vertices = vertices, |
---|
106 | boundary = tag_dict, |
---|
107 | tagged_elements = tagged_elements_dict, |
---|
108 | geo_reference = geo_reference ) |
---|
109 | |
---|
110 | |
---|
111 | |
---|
112 | #FIXME (Ole): Is this really the right place to apply the a default |
---|
113 | #value specific to the shallow water wave equation? |
---|
114 | #The 'assert' above indicates that any subclass of Domain is acceptable. |
---|
115 | #Suggestion - module shallow_water.py will eventually take care of this |
---|
116 | #(when I get around to it) so it should be removed from here. |
---|
117 | |
---|
118 | # This doesn't work on the domain instance. |
---|
119 | # This is still needed so -ve elevations don't cuase 'lakes' |
---|
120 | # The fixme we discussed was to only create a quantity when its values |
---|
121 | #are set. |
---|
122 | # I think that's the way to go still |
---|
123 | |
---|
124 | # set the water stage to be the elevation |
---|
125 | if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'): |
---|
126 | vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation'] |
---|
127 | |
---|
128 | domain.set_quantity_vertices_dict(vertex_quantity_dict) |
---|
129 | #print "vertex_quantity_dict",vertex_quantity_dict |
---|
130 | return domain |
---|
131 | |
---|
132 | |
---|
133 | def pmesh_to_domain(file_name=None, |
---|
134 | mesh_instance=None, |
---|
135 | use_cache=False, |
---|
136 | verbose=False): |
---|
137 | """ |
---|
138 | Convert a pmesh file or a pmesh mesh instance to a bunch of lists |
---|
139 | that can be used to instanciate a domain object. |
---|
140 | |
---|
141 | use_cache: True means that caching is attempted for the computed domain. |
---|
142 | """ |
---|
143 | |
---|
144 | if use_cache is True: |
---|
145 | from caching import cache |
---|
146 | result = cache(_pmesh_to_domain, (file_name, mesh_instance), |
---|
147 | dependencies = [file_name], |
---|
148 | verbose = verbose) |
---|
149 | |
---|
150 | else: |
---|
151 | result = apply(_pmesh_to_domain, (file_name, mesh_instance)) |
---|
152 | |
---|
153 | return result |
---|
154 | |
---|
155 | |
---|
156 | def _pmesh_to_domain(file_name=None, |
---|
157 | mesh_instance=None, |
---|
158 | use_cache=False, |
---|
159 | verbose=False): |
---|
160 | """ |
---|
161 | Convert a pmesh file or a pmesh mesh instance to a bunch of lists |
---|
162 | that can be used to instanciate a domain object. |
---|
163 | """ |
---|
164 | |
---|
165 | # from numpy.oldnumeric import transpose |
---|
166 | from numpy import transpose |
---|
167 | from load_mesh.loadASCII import import_mesh_file |
---|
168 | |
---|
169 | if file_name is None: |
---|
170 | mesh_dict = mesh_instance.Mesh2IODict() |
---|
171 | else: |
---|
172 | mesh_dict = import_mesh_file(file_name) |
---|
173 | #print "mesh_dict",mesh_dict |
---|
174 | vertex_coordinates = mesh_dict['vertices'] |
---|
175 | volumes = mesh_dict['triangles'] |
---|
176 | vertex_quantity_dict = {} |
---|
177 | point_atts = transpose(mesh_dict['vertex_attributes']) |
---|
178 | point_titles = mesh_dict['vertex_attribute_titles'] |
---|
179 | geo_reference = mesh_dict['geo_reference'] |
---|
180 | if point_atts != None: |
---|
181 | print 'type(point_atts)=%s' % type(point_atts) |
---|
182 | for quantity, value_vector in map(None, point_titles, point_atts): |
---|
183 | vertex_quantity_dict[quantity] = value_vector |
---|
184 | tag_dict = pmesh_dict_to_tag_dict(mesh_dict) |
---|
185 | tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict) |
---|
186 | return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference |
---|
187 | |
---|
188 | |
---|
189 | |
---|
190 | def build_tagged_elements_dictionary(mesh_dict): |
---|
191 | """Build the dictionary of element tags. |
---|
192 | tagged_elements is a dictionary of element arrays, |
---|
193 | keyed by tag: |
---|
194 | { (tag): [e1, e2, e3..] } |
---|
195 | """ |
---|
196 | tri_atts = mesh_dict['triangle_tags'] |
---|
197 | tagged_elements = {} |
---|
198 | if tri_atts is None: |
---|
199 | tagged_elements[''] = range(len(mesh_dict['triangles'])) |
---|
200 | else: |
---|
201 | for tri_att_index in range(len(tri_atts)): |
---|
202 | tagged_elements.setdefault(tri_atts[tri_att_index], |
---|
203 | []).append(tri_att_index) |
---|
204 | return tagged_elements |
---|
205 | |
---|
206 | def pmesh_dict_to_tag_dict(mesh_dict): |
---|
207 | """ Convert the pmesh dictionary (mesh_dict) description of boundary tags |
---|
208 | to a dictionary of tags, indexed with volume id and face number. |
---|
209 | """ |
---|
210 | triangles = mesh_dict['triangles'] |
---|
211 | sides = calc_sides(triangles) |
---|
212 | tag_dict = {} |
---|
213 | for seg, tag in map(None, mesh_dict['segments'], |
---|
214 | mesh_dict['segment_tags']): |
---|
215 | v1 = int(seg[0]) |
---|
216 | v2 = int(seg[1]) |
---|
217 | for key in [(v1,v2),(v2,v1)]: |
---|
218 | if sides.has_key(key) and tag <> "": |
---|
219 | #"" represents null. Don't put these into the dictionary |
---|
220 | #this creates a dict of lists of faces, indexed by tag |
---|
221 | #tagged_edges.setdefault(tag,[]).append(sides[key]) |
---|
222 | tag_dict[sides[key]] = tag |
---|
223 | |
---|
224 | return tag_dict |
---|
225 | |
---|
226 | |
---|
227 | def calc_sides(triangles): |
---|
228 | #Build dictionary mapping from sides (2-tuple of points) |
---|
229 | #to left hand side neighbouring triangle |
---|
230 | sides = {} |
---|
231 | for id, triangle in enumerate(triangles): |
---|
232 | a = int(triangle[0]) |
---|
233 | b = int(triangle[1]) |
---|
234 | c = int(triangle[2]) |
---|
235 | sides[a,b] = (id, 2) #(id, face) |
---|
236 | sides[b,c] = (id, 0) #(id, face) |
---|
237 | sides[c,a] = (id, 1) #(id, face) |
---|
238 | return sides |
---|