source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py @ 7937

Last change on this file since 7937 was 7937, checked in by steve, 14 years ago

Fixed an error in pmesh2domain

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