source: branches/numpy/anuga/abstract_2d_finite_volumes/pmesh2domain.py @ 6883

Last change on this file since 6883 was 6689, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk to numpy branch.

File size: 9.2 KB
RevLine 
[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
9import sys
[6304]10import 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).
18def 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.
59def 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]86def _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 ??
138def 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 ??
164def _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
199def 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]218def 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
239def 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
Note: See TracBrowser for help on using the repository browser.