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