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

Last change on this file since 4712 was 4098, checked in by jakeman, 17 years ago

Wrapped elements extracted from Numeric arrays in int() as they appeared
as arrays on the ANU 64 bit machine. This is similar to what was done in
changeset:2778

Committed by Ole after joint debugging with John

o-This line, and those below, will be ignored--

M pmesh2domain.py

File size: 8.1 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    for quantity, value_vector in map (None, point_titles, point_atts):
178        vertex_quantity_dict[quantity] = value_vector
179    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
180    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
181    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
182
183
184
185def build_tagged_elements_dictionary(mesh_dict):
186    """Build the dictionary of element tags.
187    tagged_elements is a dictionary of element arrays,
188    keyed by tag:
189    { (tag): [e1, e2, e3..] }
190    """
191    tri_atts = mesh_dict['triangle_tags']
192    #print "tri_atts", tri_atts
193    tagged_elements = {}
194    for tri_att_index in range(len(tri_atts)):
195        tagged_elements.setdefault(tri_atts[tri_att_index],[]).append(tri_att_index)
196    #print "DSG pm2do tagged_elements", tagged_elements
197    return tagged_elements
198
199def pmesh_dict_to_tag_dict(mesh_dict):
200    """ Convert the pmesh dictionary (mesh_dict) description of boundary tags
201    to a dictionary of tags, indexed with volume id and face number.
202    """
203    triangles = mesh_dict['triangles']
204    sides = calc_sides(triangles)
205    tag_dict = {}
206    for seg, tag in map(None, mesh_dict['segments'],
207                        mesh_dict['segment_tags']):
208        v1 = int(seg[0])
209        v2 = int(seg[1])
210        for key in [(v1,v2),(v2,v1)]:
211            if sides.has_key(key) and tag <> "":
212                #"" represents null.  Don't put these into the dictionary
213                #this creates a dict of lists of faces, indexed by tag
214                #tagged_edges.setdefault(tag,[]).append(sides[key])
215                tag_dict[sides[key]] = tag
216
217    return tag_dict
218
219
220def calc_sides(triangles):
221    #Build dictionary mapping from sides (2-tuple of points)
222    #to left hand side neighbouring triangle
223    sides = {}
224    for id, triangle in enumerate(triangles):
225        a = int(triangle[0])
226        b = int(triangle[1])
227        c = int(triangle[2])
228        sides[a,b] = (id, 2) #(id, face)
229        sides[b,c] = (id, 0) #(id, face)
230        sides[c,a] = (id, 1) #(id, face)
231    return sides
Note: See TracBrowser for help on using the repository browser.