source: anuga_core/source/anuga/abstract_2d_finite_volumes/pmesh2domain.py @ 5521

Last change on this file since 5521 was 4902, checked in by duncan, 16 years ago

Removing some of the old pmesh mesh triangulation data structure

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