source: anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/pmesh2domain.py @ 5899

Last change on this file since 5899 was 5899, checked in by rwilson, 15 years ago

Initial NumPy? changes (again!).

File size: 8.4 KB
Line 
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
11import sys
12
13
14
15def 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
55def 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
81def _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
133def 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
156def _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
190def 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
206def 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
227def 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
Note: See TracBrowser for help on using the repository browser.